diff --git a/src/lib/trace.js b/src/lib/trace.js index 6a01341..ca459de 100644 --- a/src/lib/trace.js +++ b/src/lib/trace.js @@ -269,6 +269,1018 @@ function simplify(pts, eps) { return pts.filter((_, i) => keep[i]); } +// --- fitting: straight where the artwork is straight, curved where it curves -- +// +// RDP alone cannot do this, and the reason is a measurement, not a matter of +// taste. On these masks a MATHEMATICALLY STRAIGHT edge still wanders through the +// alpha field -- print grain, JPEG ringing, and the shadow un-blend being a +// two-colour model of a soft gradient. Measured over the acceptance set, that +// wander peaks around 0.15-0.3px across a 40px run on a clean source. RDP's +// tolerance is 0.2px, i.e. BELOW the noise floor, so it faithfully reproduces +// the grain: the straight bottom of a "C" becomes a chain of short segments that +// reads as a ripple. Tightening eps makes it worse, and loosening it far enough +// to clear the noise would start cutting real curvature. +// +// Approximating each run with a cubic under a LOOSER tolerance is not enough on +// its own, and it is worth recording why, because it looks like it should be. +// A cubic fit pins its endpoints to actual sample points, so every join inherits +// that one point's share of the noise -- measured, that leaves the straight edges +// no better than RDP left them (ripple 0.56px -> 0.64px, i.e. slightly worse). +// Raising the tolerance does not help: the endpoints are pinned regardless. +// +// What removes the noise is FITTING A LINE to the whole run and throwing the +// sample points away. A total-least-squares line over a few hundred samples has +// an error of noise/sqrt(n) -- hundredths of a pixel, not tenths -- so the edge +// comes out as straight as the artwork meant it. Vertices are then the +// INTERSECTIONS of neighbouring primitives rather than sample points, which has +// a second payoff: the trace-time blur rounds every real corner into a small arc, +// and intersecting the two adjacent lines puts the sharp corner back. +// +// So the pass is: find the corners, cut the loop at them, split each piece into +// maximal straight runs and the curved gaps between them, fit lines to the +// former and cubics to the latter, then resolve the joints. Corner detection +// stays conservative -- a false corner is a visible kink in a smooth bowl, which +// is worse than the facets this exists to remove. + +/** Cumulative arc length; s[i] is the distance from pts[0] to pts[i]. */ +function arcLengths(pts) { + const n = pts.length; + const s = new Float64Array(n); + for (let i = 1; i < n; i++) + s[i] = s[i - 1] + Math.hypot(pts[i][0] - pts[i - 1][0], pts[i][1] - pts[i - 1][1]); + return s; +} + +/** + * Unit tangent at `i` over an arc-length window, by total-least-squares. + * + * A two-point difference is useless here: neighbouring contour points are ~1px + * apart and carry ~0.3px of jitter, so its direction is wrong by up to 20 + * degrees. Regressing the principal direction over a window several px long + * averages that down to a couple of degrees, which is what corner detection and + * the Bezier endpoints both need. + * + * `side` picks the window: -1 backward, +1 forward, 0 centred. `wrap` treats the + * index space as cyclic (closed loops); otherwise the window is clamped. + */ +function tangentAt(pts, s, i, win, side, wrap) { + const n = pts.length; + const total = wrap ? s[n - 1] + Math.hypot(pts[0][0] - pts[n - 1][0], pts[0][1] - pts[n - 1][1]) : 0; + const idx = (k) => (wrap ? ((k % n) + n) % n : k < 0 ? 0 : k > n - 1 ? n - 1 : k); + // Arc distance from i to k, signed, cyclic-aware. + const arc = (k) => { + let d = s[idx(k)] - s[i]; + if (wrap) { + const kk = Math.floor(k / n); + d += kk * total; + } + return d; + }; + const lo = side > 0 ? 0 : -win; + const hi = side < 0 ? 0 : win; + let k0 = i, k1 = i; + while (arc(k0 - 1) >= lo && i - k0 < n - 1) k0--; + while (arc(k1 + 1) <= hi && k1 - i < n - 1) k1++; + if (k1 - k0 < 1) { k0 = i - 1; k1 = i + 1; } + + let mx = 0, my = 0, cnt = 0; + for (let k = k0; k <= k1; k++) { + const p = pts[idx(k)]; + mx += p[0]; my += p[1]; cnt++; + } + mx /= cnt; my /= cnt; + let sxx = 0, sxy = 0, syy = 0; + for (let k = k0; k <= k1; k++) { + const p = pts[idx(k)]; + const dx = p[0] - mx, dy = p[1] - my; + sxx += dx * dx; sxy += dx * dy; syy += dy * dy; + } + // Principal eigenvector of the 2x2 covariance, in closed form. + const tr = sxx + syy; + const det = sxx * syy - sxy * sxy; + const disc = Math.max(0, (tr * tr) / 4 - det); + const lam = tr / 2 + Math.sqrt(disc); + let vx = lam - syy, vy = sxy; + if (Math.abs(vx) + Math.abs(vy) < 1e-12) { vx = sxy; vy = lam - sxx; } + let len = Math.hypot(vx, vy); + if (len < 1e-12) { + const a = pts[idx(k0)], b = pts[idx(k1)]; + vx = b[0] - a[0]; vy = b[1] - a[1]; + len = Math.hypot(vx, vy) || 1; + } + vx /= len; vy /= len; + // TLS gives an undirected line; orient it along increasing index. + const a = pts[idx(k0)], b = pts[idx(k1)]; + if ((b[0] - a[0]) * vx + (b[1] - a[1]) * vy < 0) { vx = -vx; vy = -vy; } + return [vx, vy]; +} + +/** + * Indices where the contour genuinely turns a corner. + * + * Conservative on purpose. A false corner puts a visible kink in what should be + * a smooth bowl, which is worse than the facets this whole pass exists to + * remove; a missed corner only costs a slightly soft terminal. So the turn is + * measured over a window well wider than the blur (which rounds every real + * corner into a small arc), and only local maxima above the threshold survive. + */ +function findCorners(pts, s, win, minTurn) { + const n = pts.length; + const turn = new Float64Array(n); + for (let i = 0; i < n; i++) { + const tb = tangentAt(pts, s, i, win, -1, true); + const tf = tangentAt(pts, s, i, win, +1, true); + turn[i] = Math.atan2(tb[0] * tf[1] - tb[1] * tf[0], tb[0] * tf[0] + tb[1] * tf[1]); + } + // Non-maximum suppression over the same window, so one corner reports once + // instead of as a cluster of neighbours that all clear the threshold. + const out = []; + const total = s[n - 1] + Math.hypot(pts[0][0] - pts[n - 1][0], pts[0][1] - pts[n - 1][1]); + const arcGap = (a, b) => { + let d = Math.abs(s[a] - s[b]); + return Math.min(d, total - d); + }; + for (let i = 0; i < n; i++) { + const m = Math.abs(turn[i]); + if (m < minTurn) continue; + let best = true; + for (let k = 1; k < n; k++) { + const j = (i + k) % n; + if (arcGap(i, j) > win) break; + if (Math.abs(turn[j]) > m) { best = false; break; } + } + if (!best) continue; + for (let k = 1; k < n; k++) { + const j = ((i - k) % n + n) % n; + if (arcGap(i, j) > win) break; + // Ties go to the earlier index so a plateau reports exactly once. + if (Math.abs(turn[j]) >= m) { best = false; break; } + } + if (best) out.push(i); + } + return out; +} + +/** + * Prefix moments, so a total-least-squares line over ANY index range costs O(1). + * + * Growing a straight run means testing many candidate ranges; doing that with an + * O(n) fit each time is quadratic on exactly the long straight edges this pass + * exists for. Coordinates are shifted to the loop's centroid first -- the raw + * ones run to a couple of thousand, and squaring them before differencing throws + * away precision the variance of a near-straight run cannot spare. + */ +function momentTable(pts) { + const n = pts.length; + let cx = 0, cy = 0; + for (const p of pts) { cx += p[0]; cy += p[1]; } + cx /= n; cy /= n; + const sx = new Float64Array(n + 1), sy = new Float64Array(n + 1); + const sxx = new Float64Array(n + 1), sxy = new Float64Array(n + 1), syy = new Float64Array(n + 1); + for (let i = 0; i < n; i++) { + const x = pts[i][0] - cx, y = pts[i][1] - cy; + sx[i + 1] = sx[i] + x; + sy[i + 1] = sy[i] + y; + sxx[i + 1] = sxx[i] + x * x; + sxy[i + 1] = sxy[i] + x * y; + syy[i + 1] = syy[i] + y * y; + } + return { sx, sy, sxx, sxy, syy, cx, cy }; +} + +/** + * Algebraic (Kasa) circle through a list of points: centre, radius, and the + * GEOMETRIC rms residual. Minimising the algebraic residual is linear in the + * three unknowns, so this is one 3x3 solve; the geometric rms is then measured + * directly, because the algebraic one is scaled by the radius and would call a + * huge-radius fit good no matter how badly it matched. + */ +function fitCircle(win) { + const n = win.length; + if (n < 6) return null; + let mx = 0, my = 0; + for (const p of win) { mx += p[0]; my += p[1]; } + mx /= n; my /= n; + let Sxx = 0, Sxy = 0, Syy = 0, Sxz = 0, Syz = 0, Sz = 0; + for (const p of win) { + const x = p[0] - mx, y = p[1] - my, z = x * x + y * y; + Sxx += x * x; Sxy += x * y; Syy += y * y; + Sxz += x * z; Syz += y * z; Sz += z; + } + // Centred, so the linear terms drop out and the solve is 2x2. + const det = Sxx * Syy - Sxy * Sxy; + if (Math.abs(det) < 1e-12) return null; + const cx = (Sxz * Syy - Syz * Sxy) / (2 * det) + mx; + const cy = (Syz * Sxx - Sxz * Sxy) / (2 * det) + my; + let sr = 0; + for (const p of win) sr += Math.hypot(p[0] - cx, p[1] - cy); + const r = sr / n; + if (!(r > 0) || !Number.isFinite(r)) return null; + let se = 0; + for (const p of win) { + const d = Math.hypot(p[0] - cx, p[1] - cy) - r; + se += d * d; + } + return { cx, cy, r, rms: Math.sqrt(se / n) }; +} + +/** + * Is pts[a..b] really a CHORD of an arc that carries on through it? + * + * A straight-run test built on a residual cannot answer this on its own. On a + * circle of radius R a chord of arc length L departs from its own best-fit line + * by only about 0.3 * L^2 / (8R), so at any tolerance loose enough to absorb + * print grain, a big enough bowl admits chords -- which is exactly how the + * counter of an "O" comes out as a polygon. + * + * The evidence that settles it lies OUTSIDE the run. Fit a circle to a window + * extended past both ends: if that window is a clean arc, and the arc's sagitta + * across the run is more than the eye would forgive, the run is a chord and not + * a straight edge. A genuinely straight edge either has no room to extend + * (corners bound it) or extends into something no circle fits. + * + * On a loop with no corners the window WRAPS. It has to: with nothing to cut at, + * the seam goes at the flattest point of the outline, which is precisely where a + * straight run wants to form -- clipping there would exempt the one place that + * most needs the test, and leave a flat spot at the seam of every "O". + */ +function isChordOfArc(q, s, a, b, from, to, wrap, N, o) { + const L = s[b] - s[a]; + if (L <= 0) return false; + const reach = L * 0.6; + const at = (i) => q[((i % N) + N) % N]; + + // Collect the window, walking outward by arc length from each end. + const back = [], fwd = []; + let got = 0; + for (let i = a - 1; ; i--) { + if (!wrap && i < from) break; + if (a - i >= N) break; + const p = at(i), pv = at(i + 1); + got += Math.hypot(pv[0] - p[0], pv[1] - p[1]); + back.push(p); + if (got >= reach) break; + } + const backArc = got; + got = 0; + for (let i = b + 1; ; i++) { + if (!wrap && i > to) break; + if (i - b >= N) break; + const p = at(i), pv = at(i - 1); + got += Math.hypot(p[0] - pv[0], p[1] - pv[1]); + fwd.push(p); + if (got >= reach) break; + } + // Without context on both sides there is nothing to infer from. + if (backArc < L * 0.25 || got < L * 0.25) return false; + + const win = back.reverse(); + for (let i = a; i <= b; i++) win.push(q[i]); + win.push(...fwd); + + const c = fitCircle(win); + if (!c) return false; + // The window has to actually BE an arc; "straight, then it bends away" is + // not, and a circle fitted to that would invent a radius meaning nothing. + if (c.rms > o.lineRms * 1.3) return false; + + // Sagitta of that arc across the run itself. + const chord = Math.hypot(q[b][0] - q[a][0], q[b][1] - q[a][1]); + const half = Math.min(c.r, chord / 2); + return c.r - Math.sqrt(Math.max(0, c.r * c.r - half * half)) > o.arcSagitta; +} + +/** + * TLS line over pts[a..b], with the RMS perpendicular residual. + * + * The residual is the smaller eigenvalue of the covariance, which is why this + * doubles as the straightness test: zero-mean grain contributes its own variance + * and nothing more, while a genuine bend contributes a systematic term that + * grows with the square of the run length. One threshold separates them. + */ +function lineOver(M, a, b) { + const n = b - a + 1; + if (n < 3) return null; + const mx = (M.sx[b + 1] - M.sx[a]) / n; + const my = (M.sy[b + 1] - M.sy[a]) / n; + const cxx = (M.sxx[b + 1] - M.sxx[a]) / n - mx * mx; + const cxy = (M.sxy[b + 1] - M.sxy[a]) / n - mx * my; + const cyy = (M.syy[b + 1] - M.syy[a]) / n - my * my; + const tr = cxx + cyy; + const disc = Math.sqrt(Math.max(0, (tr * tr) / 4 - (cxx * cyy - cxy * cxy))); + const lmax = tr / 2 + disc; + const lmin = Math.max(0, tr / 2 - disc); + let vx = lmax - cyy, vy = cxy; + if (Math.abs(vx) + Math.abs(vy) < 1e-12) { vx = cxy; vy = lmax - cxx; } + const len = Math.hypot(vx, vy); + if (!(len > 1e-12)) return null; + return { + ox: mx + M.cx, oy: my + M.cy, + dx: vx / len, dy: vy / len, + rms: Math.sqrt(lmin), + }; +} + +/** + * Point the line along increasing index. The principal eigenvector defines an + * undirected axis, and the sign it happens to come out with is arbitrary -- but + * the line doubles as the tangent handed to the neighbouring curve fit, and a + * reversed tangent there sends the cubic off backwards. + */ +function orient(ln, a, b) { + if (!ln) return null; + if ((b[0] - a[0]) * ln.dx + (b[1] - a[1]) * ln.dy < 0) { ln.dx = -ln.dx; ln.dy = -ln.dy; } + return ln; +} + +/** + * Undo greedy over-segmentation of a curve into chords. + * + * `straightRuns` scans forward and starts a new run wherever the last one stopped + * fitting, so a gently curving arc comes back as a chain of straight pieces + * meeting at small angles -- a polygon, which is the very defect this pass exists + * to remove. Any two adjacent lines whose junction is not a corner belong to one + * curve and are merged back. + * + * The angle threshold here can be far lower than the one `findCorners` uses, and + * that is the point: a corner between two FITTED LINES is measured from hundreds + * of averaged samples, not from two noisy windowed tangents, so a break this + * shallow is still real when both sides are confidently straight. + */ +function unpolygonise(pieces, isCorner, o) { + const out = []; + let i = 0; + while (i < pieces.length) { + if (pieces[i].kind !== "L") { out.push(pieces[i]); i++; continue; } + let j = i; + while ( + j + 1 < pieces.length && + pieces[j + 1].kind === "L" && + pieces[j + 1].a === pieces[j].b && + !isCorner.has(pieces[j].b) && + pieces[j].line && pieces[j + 1].line && + Math.abs(Math.atan2( + pieces[j].line.dx * pieces[j + 1].line.dy - pieces[j].line.dy * pieces[j + 1].line.dx, + pieces[j].line.dx * pieces[j + 1].line.dx + pieces[j].line.dy * pieces[j + 1].line.dy, + )) < o.lineCornerTurn + ) j++; + if (j > i) out.push({ kind: "C", a: pieces[i].a, b: pieces[j].b }); + else out.push(pieces[i]); + i = j + 1; + } + return out; +} + +/** Point on `ln` closest to (x,y). */ +function projectOn(ln, x, y) { + const t = (x - ln.ox) * ln.dx + (y - ln.oy) * ln.dy; + return [ln.ox + ln.dx * t, ln.oy + ln.dy * t]; +} + +/** Where two lines cross, or null if they are too close to parallel to trust. */ +function intersectLines(p, q) { + const den = p.dx * q.dy - p.dy * q.dx; + if (Math.abs(den) < 0.12) return null; // < ~7 degrees apart + const t = ((q.ox - p.ox) * q.dy - (q.oy - p.oy) * q.dx) / den; + return [p.ox + p.dx * t, p.oy + p.dy * t]; +} + +/** + * Pull a straight run's ends back to where the contour is genuinely parallel to + * the fitted line. + * + * The run's extent comes from an RMS test, and RMS is an average: a run may + * absorb the first few pixels of a bend at each end and still pass, because the + * long straight middle dominates. The line then chords across the start of the + * curve, and the neighbouring cubic is left too little room to bend, so the + * handover shows as a kink -- which on a stem-to-shoulder junction is the most + * conspicuous place on a letterform to put one. + * + * Trimming on DIRECTION rather than on distance is what makes this work: at the + * moment the contour starts to turn it is still sitting almost exactly on the + * line, so a distance test cannot see it, while the direction has already moved. + * What is trimmed off is not lost -- it becomes part of the adjacent curve, + * which is where it belonged. + * + * The direction has to come from a LINE FITTED OVER A WINDOW, not from a + * point tangent. A tangent taken over a few pixels of a grainy edge carries + * several degrees of noise all by itself, so testing it against a few-degree + * tolerance trims at random and eats genuinely straight edges -- measured, that + * put the ripple back up from 0.19px to 0.48px. A window of a few dozen pixels + * averages that noise down to well under a degree. + */ +function trimToTangent(q, s, M, a, b, line, o, holdHead, holdTail) { + // Undirected angle between a fitted window and the run's own line. + const skew = (ln) => { + const c = Math.abs(ln.dx * line.dy - ln.dy * line.dx); + const d = Math.abs(ln.dx * line.dx + ln.dy * line.dy); + return Math.atan2(c, d); + }; + // The window's own direction is only known to within the grain it contains: + // a least-squares slope over n samples spanning W with noise sigma is + // uncertain by about sigma*sqrt(12/n)/W. Testing against a FIXED angle + // therefore means something different on a clean print than on a rough one, + // and on the rough ones it trims edges that are perfectly straight -- + // measured, that alone put Camping In Alaska's ripple from 0.16px to 0.44px + // while barely touching Hot Toddy's. So the floor rises with the run's own + // measured noise, and the same threshold means the same thing everywhere. + const angNoise = (ln, n) => + (3 * ln.rms * Math.sqrt(12 / Math.max(4, n))) / Math.max(1e-6, o.trimWin); + let a2 = a, b2 = b; + // An end sitting ON a detected corner does not get the tangent trim: the + // contour runs straight into a corner by definition, so there is no bend to + // back out of, and that end's direction is exactly what the corner is + // reconstructed from. + // + // (The blur-rounded apex is still kept out of the DIRECTION the line is + // fitted from -- see `lineForRun` -- but that is a narrower fit range, not a + // shorter run, so it costs no straight edge.) + for (let guard = 0; guard < 24 && !holdHead; guard++) { + let e = a2; + while (e < b2 && s[e] - s[a2] < o.trimWin) e++; + if (e >= b2 - 2 || s[e] - s[a2] < o.trimWin * 0.5) break; + const ln = lineOver(M, a2, e); + if (!ln) break; + if (skew(ln) <= Math.max(o.lineTangentTol, angNoise(ln, e - a2 + 1))) break; + a2 += Math.max(1, (e - a2) >> 1); + } + for (let guard = 0; guard < 24 && !holdTail; guard++) { + let e = b2; + while (e > a2 && s[b2] - s[e] < o.trimWin) e--; + if (e <= a2 + 2 || s[b2] - s[e] < o.trimWin * 0.5) break; + const ln = lineOver(M, e, b2); + if (!ln) break; + if (skew(ln) <= Math.max(o.lineTangentTol, angNoise(ln, b2 - e + 1))) break; + b2 -= Math.max(1, (b2 - e) >> 1); + } + return [a2, b2]; +} + +/** + * The line a straight run is drawn along, fitted from the part of the run that + * is clean. + * + * Where the run ends at a CORNER, the last couple of pixels are the apex the + * blur rounded into a little arc. Those points still belong to the run -- the + * edge really does go that far -- but letting them into the FIT tilts the + * direction, and the corner is reconstructed by intersecting two of these + * directions. At an acute wedge the two lines cross at a shallow angle, where a + * fraction of a degree of tilt moves the crossing a long way past the true + * apex: a spike sticking out of the letterform, which is the most visible thing + * this pass can get wrong. Measured on the apex of Mildlife's "M", leaving the + * rounding in put the reconstructed point 11px out. + * + * Narrowing the fit range costs nothing, because a line piece's endpoints come + * from its joints, not from the samples it was fitted to. + */ +function lineForRun(q, s, M, a, b, o, headCorner, tailCorner) { + let fa = a, fb = b; + if (headCorner) while (fa < fb - 3 && s[fa] - s[a] < o.cornerGuard) fa++; + if (tailCorner) while (fb > fa + 3 && s[b] - s[fb] < o.cornerGuard) fb--; + return orient(lineOver(M, fa, fb) ?? lineOver(M, a, b), q[a], q[b]); +} + +/** + * Maximal straight runs within pts[from..to], as [a, b] index pairs. + * + * Greedy from `from`, which is a corner (or the seam), so runs start where the + * geometry actually starts. Growth doubles then binary-refines, and each probe + * is O(1) against the moment table, so a 400px straight edge costs a couple of + * dozen probes rather than a scan per candidate endpoint. + */ +function straightRuns(pts, s, M, from, to, wrap, N, o, atCorner, toCorner) { + const out = []; + let a = from; + const fits = (i, j) => { + const ln = lineOver(M, i, j); + return ln && ln.rms <= o.lineRms ? ln : null; + }; + while (a < to - 2) { + let b = Math.min(to, a + 3); + if (!fits(a, b)) { a++; continue; } + let step = 4; + while (b + step <= to && fits(a, b + step)) { b += step; step *= 2; } + for (let jump = step >> 1; jump >= 1; jump >>= 1) + if (b + jump <= to && fits(a, b + jump)) b += jump; + let keep = null, head = a; + if (s[b] - s[a] >= o.minLineLen && !isChordOfArc(pts, s, a, b, from, to, wrap, N, o)) { + // Pull the ends back out of the bends the RMS average absorbed, then + // re-fit: the shortened run has a slightly different direction, and + // it is that direction the neighbouring curve will be handed. + const ln0 = orient(lineOver(M, a, b), pts[a], pts[b]); + if (ln0) { + const [ta, tb] = trimToTangent( + pts, s, M, a, b, ln0, o, + a === from && atCorner, b === to && toCorner, + ); + head = ta; + // A run that trims away was never a straight edge -- it was the + // gentle part of a bend that happened to average out flat. + if (s[tb] - s[ta] >= o.minLineLen) keep = [ta, tb]; + } + } + if (keep) { + out.push(keep); + a = keep[1]; + } else if (head > a) { + // The scan began inside a bend and the straight part starts later. + // Resume exactly there rather than stepping blindly past it, or the + // straight edge that follows is never found and gets drawn as curve. + a = head; + } else { + // Not a straight edge. Step past most of what was examined; retrying + // from a+1 would be quadratic for no gain. + a += Math.max(1, (b - a) >> 1); + } + } + return out; +} + +const B0 = (u) => (1 - u) ** 3; +const B1 = (u) => 3 * u * (1 - u) ** 2; +const B2 = (u) => 3 * u * u * (1 - u); +const B3 = (u) => u * u * u; + +function bezPoint(p0, p1, p2, p3, u) { + const a = B0(u), b = B1(u), c = B2(u), d = B3(u); + return [ + a * p0[0] + b * p1[0] + c * p2[0] + d * p3[0], + a * p0[1] + b * p1[1] + c * p2[1] + d * p3[1], + ]; +} + +/** + * Least-squares cubic through `pts[i0..i1]` with the endpoints pinned and the + * end tangents fixed (Schneider, Graphics Gems 1990). Solving only for the two + * tangent MAGNITUDES is what keeps neighbouring segments joining smoothly: + * direction is decided once, by the tangent estimate, and shared across the join. + */ +function fitOne(pts, i0, i1, t0, t1, u, p0, p3) { + let c00 = 0, c01 = 0, c11 = 0, x0 = 0, x1 = 0; + for (let k = i0; k <= i1; k++) { + const uu = u[k - i0]; + const a0x = t0[0] * B1(uu), a0y = t0[1] * B1(uu); + const a1x = t1[0] * B2(uu), a1y = t1[1] * B2(uu); + c00 += a0x * a0x + a0y * a0y; + c01 += a0x * a1x + a0y * a1y; + c11 += a1x * a1x + a1y * a1y; + const base = B0(uu) + B1(uu), tail = B2(uu) + B3(uu); + const rx = pts[k][0] - (p0[0] * base + p3[0] * tail); + const ry = pts[k][1] - (p0[1] * base + p3[1] * tail); + x0 += a0x * rx + a0y * ry; + x1 += a1x * rx + a1y * ry; + } + const det = c00 * c11 - c01 * c01; + const chord = Math.hypot(p3[0] - p0[0], p3[1] - p0[1]); + let a0, a1; + if (Math.abs(det) < 1e-12) { + a0 = a1 = chord / 3; + } else { + a0 = (x0 * c11 - x1 * c01) / det; + a1 = (c00 * x1 - c01 * x0) / det; + } + // Wu/Barsky fallback: a non-positive or wildly long handle means the normal + // equations found a solution that is not a sane arc, usually because the run + // doubles back. The chord/3 heuristic is the standard, well-behaved answer. + const lim = chord * 1.5; + if (!(a0 > 1e-6) || !(a1 > 1e-6) || a0 > lim || a1 > lim) a0 = a1 = chord / 3; + return [ + p0, + [p0[0] + t0[0] * a0, p0[1] + t0[1] * a0], + [p3[0] + t1[0] * a1, p3[1] + t1[1] * a1], + p3, + ]; +} + +/** + * How well the curve fits, and where to cut if it does not. + * + * Judging on the worst single point is a trap at these sample counts. A run of + * 250 points carrying s.d. sigma of grain has an EXPECTED worst point around + * 3.3 sigma just from the noise -- about 0.8px here -- so a fit that is already + * as good as the data allows keeps failing a max-residual test and splitting. + * Worse, the split then lands on whichever point the noise spiked at, planting a + * knot in the middle of a smooth arc. + * + * So acceptance is on RMS, with the max kept only as an outlier ceiling, and the + * split point is taken from the residual after a short moving average -- which + * finds where the curve is SYSTEMATICALLY off and ignores single spikes. + */ +function fitError(pts, i0, i1, bez, u) { + let max = 0, sum2 = 0; + for (let k = i0; k <= i1; k++) { + const q = bezPoint(bez[0], bez[1], bez[2], bez[3], u[k - i0]); + const dx = q[0] - pts[k][0], dy = q[1] - pts[k][1]; + const d2 = dx * dx + dy * dy; + sum2 += d2; + if (d2 > max && k > i0 && k < i1) max = d2; + } + return { max: Math.sqrt(max), rms: Math.sqrt(sum2 / (i1 - i0 + 1)) }; +} + +/** + * Where to cut a cubic that missed tolerance: the peak of the residual after a + * short moving average, so the knot lands where the curve is SYSTEMATICALLY off + * rather than wherever a single sample spiked. Only worth computing on the + * fits that are actually going to split, which is a minority of them. + */ +function splitAt(pts, s, i0, i1, bez, u, win) { + const n = i1 - i0 + 1; + const res = new Float64Array(n); + for (let k = i0; k <= i1; k++) { + const q = bezPoint(bez[0], bez[1], bez[2], bez[3], u[k - i0]); + res[k - i0] = Math.hypot(q[0] - pts[k][0], q[1] - pts[k][1]); + } + let best = -1, at = (i0 + i1) >> 1, lo = 0, hi = 0, run = 0; + for (let k = 0; k < n; k++) { + const a = s[i0 + k] - win, b = s[i0 + k] + win; + while (hi < n && s[i0 + hi] <= b) { run += res[hi]; hi++; } + while (lo < n && s[i0 + lo] < a) { run -= res[lo]; lo++; } + const cnt = hi - lo; + if (cnt <= 0) continue; + const v = run / cnt; + if (v > best && k > 1 && k < n - 2) { best = v; at = i0 + k; } + } + return at; +} + +/** + * One Newton step per point, pulling each parameter toward the curve's closest + * point. Chord-length parameterisation is only an estimate of arc position, and + * on a strongly curved run the error it leaves behind can force a split that the + * geometry does not actually need. + */ +function reparam(pts, i0, i1, bez, u) { + const [p0, p1, p2, p3] = bez; + const d1 = [ + [3 * (p1[0] - p0[0]), 3 * (p1[1] - p0[1])], + [3 * (p2[0] - p1[0]), 3 * (p2[1] - p1[1])], + [3 * (p3[0] - p2[0]), 3 * (p3[1] - p2[1])], + ]; + const d2 = [ + [2 * (d1[1][0] - d1[0][0]), 2 * (d1[1][1] - d1[0][1])], + [2 * (d1[2][0] - d1[1][0]), 2 * (d1[2][1] - d1[1][1])], + ]; + const out = new Float64Array(u.length); + for (let k = i0; k <= i1; k++) { + const uu = u[k - i0]; + const q = bezPoint(p0, p1, p2, p3, uu); + const v = 1 - uu; + const qx = v * v * d1[0][0] + 2 * v * uu * d1[1][0] + uu * uu * d1[2][0]; + const qy = v * v * d1[0][1] + 2 * v * uu * d1[1][1] + uu * uu * d1[2][1]; + const qxx = v * d2[0][0] + uu * d2[1][0]; + const qyy = v * d2[0][1] + uu * d2[1][1]; + const dx = q[0] - pts[k][0], dy = q[1] - pts[k][1]; + const num = dx * qx + dy * qy; + const den = qx * qx + qy * qy + dx * qxx + dy * qyy; + let nu = Math.abs(den) < 1e-12 ? uu : uu - num / den; + if (!(nu >= 0)) nu = 0; + if (nu > 1) nu = 1; + out[k - i0] = nu; + } + // Parameters must stay increasing or the next solve is garbage. + for (let k = 1; k < out.length; k++) + if (out[k] <= out[k - 1]) out[k] = Math.min(1, out[k - 1] + 1e-6); + return out; +} + +/** + * Fit `pts[i0..i1]` with as few cubics as the tolerance allows, splitting at the + * worst point and recursing. `t0`/`t1` are the unit tangents the ends must obey. + */ +function fitRun(pts, s, i0, i1, t0, t1, o, out, depth, P0, P3) { + const p0 = P0 || pts[i0]; + const p3 = P3 || pts[i1]; + const n = i1 - i0 + 1; + if (n < 3) { out.push([p3]); return; } + + // Chord-length parameterisation, the standard starting estimate. + const span = s[i1] - s[i0]; + let u = new Float64Array(n); + if (span > 1e-9) for (let k = 0; k < n; k++) u[k] = (s[i0 + k] - s[i0]) / span; + else for (let k = 0; k < n; k++) u[k] = k / (n - 1); + + const ok = (e) => e.rms <= o.fitRms && e.max <= o.fitMax; + let bez = fitOne(pts, i0, i1, t0, t1, u, p0, p3); + let err = fitError(pts, i0, i1, bez, u); + // Reparameterise when the fit is short of tolerance but within reach. The + // gate is on RMS, not on the worst sample: keyed to the max it would be a + // function of the grain and of how many points the run happens to hold, and + // would skip the long strongly-curved runs where chord-length + // parameterisation is worst and Newton is worth most. + const reach = o.fitRms * o.reparamGate; + for (let it = 0; it < o.reparamIters && !ok(err) && err.rms <= reach; it++) { + const nu = reparam(pts, i0, i1, bez, u); + const nb = fitOne(pts, i0, i1, t0, t1, nu, p0, p3); + const ne = fitError(pts, i0, i1, nb, nu); + if (ne.rms >= err.rms) break; + u = nu; bez = nb; err = ne; + } + if (ok(err) || depth >= o.maxDepth) { out.push(bez.slice(1)); return; } + + // Split at the worst point, but never so close to an end that the sub-run is + // too short to estimate a tangent from. + let at = splitAt(pts, s, i0, i1, bez, u, o.tanWin); + const minPts = 2; + if (at - i0 < minPts) at = i0 + minPts; + if (i1 - at < minPts) at = i1 - minPts; + if (at <= i0 || at >= i1) { out.push(bez.slice(1)); return; } + + const tm = tangentAt(pts, s, at, o.tanWin, 0, false); + fitRun(pts, s, i0, at, t0, [-tm[0], -tm[1]], o, out, depth + 1, p0, pts[at]); + fitRun(pts, s, at, i1, tm, t1, o, out, depth + 1, pts[at], p3); +} + +/** + * Fit one closed contour into a list of path segments. + * + * Returns `{ start, segs }` where each seg is `[p]` for a line or + * `[c1, c2, p]` for a cubic, or null for a loop with no area to draw -- + * one bad subpath would otherwise poison the whole `d` string. + */ +function fitLoop(pts, o) { + let n = pts.length; + if (n < 3) return null; + + // Marching squares emits a pair of near-coincident points wherever the + // isoline clips a cell corner. They carry no information, and a zero-length + // step makes every tangent and parameterisation below divide by nothing. + const clean = [pts[0]]; + for (let i = 1; i < n; i++) { + const p = pts[i], q = clean[clean.length - 1]; + if (Math.abs(p[0] - q[0]) > 1e-4 || Math.abs(p[1] - q[1]) > 1e-4) clean.push(p); + } + while (clean.length > 2) { + const a = clean[0], b = clean[clean.length - 1]; + if (Math.abs(a[0] - b[0]) < 1e-4 && Math.abs(a[1] - b[1]) < 1e-4) clean.pop(); + else break; + } + pts = clean; + n = pts.length; + // Everything that survived the dedupe collapsed onto one place: no shape. + if (n < 3) return null; + if (n < 8) return { start: pts[0], segs: pts.slice(1).map((p) => [p]) }; + + const s0 = arcLengths(pts); + const perim = s0[n - 1] + Math.hypot(pts[0][0] - pts[n - 1][0], pts[0][1] - pts[n - 1][1]); + + // Tiny loops (a dot on an i, a punched counter) carry too few samples for a + // window-based tangent to mean anything. Round-trip them through RDP, which + // is the honest answer at that size. + if (perim < o.minPerim) { + const sp = simplify(pts, o.eps); + return { start: sp[0], segs: sp.slice(1).concat([sp[0]]).map((p) => [p]) }; + } + + const corners = findCorners(pts, s0, o.cornerWin, o.minTurn); + + // A closed contour has no natural seam. Cut at a corner when there is one -- + // the fit breaks there anyway. With none (a bowl, an "O", a rounded + // rectangle) cut at the TIGHTEST bend. + // + // The flattest point is the tempting choice and it is the wrong one: flattest + // means the middle of the longest straight run, which is the one place a cut + // costs something. Measured on a rounded rectangle, seaming there splits a + // 248px side into two independently fitted lines that meet at a slight + // angle. The tightest bend is where the fit was going to place a knot + // anyway, and both ends take the same centred tangent so the join stays + // smooth. + let seam = 0; + if (corners.length) seam = corners[0]; + else { + let bend = -1; + const stride = Math.max(1, n >> 8); + for (let i = 0; i < n; i += stride) { + const tb = tangentAt(pts, s0, i, o.cornerWin, -1, true); + const tf = tangentAt(pts, s0, i, o.cornerWin, +1, true); + const t = Math.abs(Math.atan2(tb[0] * tf[1] - tb[1] * tf[0], tb[0] * tf[0] + tb[1] * tf[1])); + if (t > bend) { bend = t; seam = i; } + } + } + + // Work on the loop opened at the seam, with the seam point repeated at the + // end so index N is the same place as index 0. + const q = pts.slice(seam).concat(pts.slice(0, seam)); + q.push(q[0]); + const N = q.length - 1; + const s = arcLengths(q); + const M = momentTable(q); + + const cut = corners.map((c) => (c - seam + n) % n).filter((c) => c > 0 && c < N); + cut.sort((a, b) => a - b); + const bounds = [0, ...cut, N]; + const isCorner = new Set(cut); + if (corners.length) isCorner.add(0); + + // --- split each corner-bounded stretch into line runs and curved gaps ---- + const pieces = []; + for (let bi = 0; bi + 1 < bounds.length; bi++) { + const from = bounds[bi], to = bounds[bi + 1]; + if (to <= from) continue; + // With no corners the loop is one unbroken stretch, so the arc test may + // look across the seam; with corners, each stretch ends at a real break + // and looking past it would read the wrong geometry. + const runs = straightRuns( + q, s, M, from, to, !cut.length, N, o, + isCorner.has(from) || (bi === 0 && isCorner.has(0)), + isCorner.has(to) || (to === N && isCorner.has(0)), + ); + const local = []; + let at = from; + for (const [a, b] of runs) { + if (a > at) local.push({ kind: "C", a: at, b: a }); + local.push({ + kind: "L", a, b, + line: lineForRun(q, s, M, a, b, o, isCorner.has(a), isCorner.has(b)), + }); + at = b; + } + if (at < to) local.push({ kind: "C", a: at, b: to }); + pieces.push(...unpolygonise(local, isCorner, o)); + } + if (!pieces.length) return null; + + // --- resolve the joints ------------------------------------------------- + // Every vertex is where two primitives MEET, not a sample point. For two + // lines that is their intersection, which is what puts a blur-rounded corner + // back; for anything involving a curve the line is the trustworthy side, so + // the sample point is projected onto it. + const m = pieces.length; + const V = new Array(m); // vertex at the START of piece i + const T = new Array(m); // unit tangent leaving that vertex, for curve fits + for (let i = 0; i < m; i++) { + const cur = pieces[i]; + const prev = pieces[(i - 1 + m) % m]; + const k = cur.a; + const corner = isCorner.has(k) || (i === 0 && isCorner.has(0)); + + if (prev.kind === "L" && cur.kind === "L" && prev.line && cur.line) { + const x = intersectLines(prev.line, cur.line); + // How far the apex may legitimately sit beyond the contour is set by + // physics, not by taste: a corner of interior angle theta rounded to + // radius r has its true apex r*(1/sin(theta/2) - 1) past the arc. At + // a right angle that is a fraction of the rounding; at an acute wedge + // it grows fast, and so does the error in it, because the two lines + // cross at a shallow angle where a fractional tilt moves the crossing + // a long way. So the shift is bounded by what the rounding can + // explain, and CLAMPED rather than discarded -- a corner pulled back + // to the limit is still sharper than the blurred contour, whereas + // falling back to the contour point throws the reconstruction away. + let v = null; + if (x) { + const dx = x[0] - q[k][0], dy = x[1] - q[k][1]; + const d = Math.hypot(dx, dy); + // The path turns by phi, so the wedge's interior angle is + // pi - phi and sin(interior/2) is cos(phi/2), which comes + // straight from the dot product via the half-angle identity. + const dot = prev.line.dx * cur.line.dx + prev.line.dy * cur.line.dy; + const cosHalf = Math.sqrt(Math.max(1e-3, (1 + dot) / 2)); + const lim = Math.min(o.cornerReach, o.cornerGuard * (1 / cosHalf - 1) + 1); + v = d <= lim || d < 1e-9 + ? x + : [q[k][0] + (dx / d) * lim, q[k][1] + (dy / d) * lim]; + } + V[i] = v ?? projectOn(cur.line, q[k][0], q[k][1]); + T[i] = [cur.line.dx, cur.line.dy]; + } else if (cur.kind === "L" && cur.line) { + V[i] = projectOn(cur.line, q[k][0], q[k][1]); + T[i] = [cur.line.dx, cur.line.dy]; + } else if (prev.kind === "L" && prev.line) { + // Curve leaving a line: start it ON the line and, unless this is a + // corner, leave along the line so the join has no visible break. + V[i] = projectOn(prev.line, q[k][0], q[k][1]); + T[i] = corner + ? tangentAt(q, s, k, o.tanWin, +1, false) + : [prev.line.dx, prev.line.dy]; + } else { + V[i] = q[k]; + T[i] = corner + ? tangentAt(q, s, k, o.tanWin, +1, false) + : tangentAt(q, s, k, o.tanWin, 0, false); + } + } + + // A line's direction is fixed, so its own start tangent must agree with it; + // for the run's END the incoming tangent is the next vertex's outgoing one + // when that vertex belongs to a line, and a backward estimate otherwise. + const endTangent = (i) => { + const nxt = pieces[(i + 1) % m]; + if (nxt.kind === "L" && nxt.line && !isCorner.has(pieces[i].b)) + return [nxt.line.dx, nxt.line.dy]; + return tangentAt(q, s, pieces[i].b, o.tanWin, -1, false); + }; + + const segs = []; + for (let i = 0; i < m; i++) { + const p = pieces[i]; + const start = V[i]; + const end = V[(i + 1) % m]; + if (p.kind === "L") { + segs.push([end]); + } else { + const te = endTangent(i); + fitRun(q, s, p.a, p.b, T[i], [-te[0], -te[1]], o, segs, 0, start, end); + } + } + return { start: V[0], segs }; +} + +/** + * Fitting tolerances, in source pixels of the raster being traced. + * + * They have to scale with the raster or they mean different things on different + * paths: a cover wordmark is traced at 2400px on its long edge, a Fanart PNG at + * its native 800. The fractions below are calibrated on the 2400px cover masks + * (where `tol` lands at 0.6px, comfortably above the measured 0.15-0.3px noise + * floor and far below the ~2px a real curve departs over the same run) and carry + * over by proportion. + */ +function fitOptions(opts) { + const edge = Math.max(opts.w || 0, opts.h || 0) || 2400; + // Cubic fit tolerance. Raising it past this buys steadily less: with the + // straight runs already fitted as lines, the 99th-percentile error barely + // moves between 0.4 and 0.9 (0.92px -> 0.93px) because the lines, not the + // cubics, set it -- while the segment count and the residual wobble both + // fall. So it sits at the loose end of the range where fidelity is flat. + const tol = opts.fitTol ?? Math.max(0.35, 0.000375 * edge); + const cornerWin = opts.cornerWin ?? Math.max(3, 0.0025 * edge); + return { + tol, + // The straight-run test, as an RMS perpendicular residual: grain on a + // straight edge contributes only its own variance, while a real bend + // contributes a systematic term, so one threshold separates them. + // + // 0.46px at a 2400px raster, chosen off the sweep. Below it the grain + // starts breaking runs that should be single lines (ripple rises from + // 0.24 to 0.40px by 0.25); above it real curves start being flattened + // into chords -- at 0.8 the worst straight-edge error jumps to 6px and + // Hot Toddy loses three quarters of its curves. 0.6 already costs more + // fidelity than the small ripple gain is worth. + lineRms: opts.lineRms ?? Math.max(0.08, 0.00019 * edge), + // Shorter than this is not worth calling a straight edge -- it just + // fragments a curve into chords. + minLineLen: opts.minLineLen ?? Math.max(6, 0.01 * edge), + // How far a surrounding arc may bow across a candidate straight run + // before the run is judged a chord of that arc rather than an edge. + // 0.36px at a 2400px raster is about a tenth of a pixel once the hero + // has scaled the mask down -- below notice, and well under the sagitta + // at which a chain of chords starts to read as a polygon. + arcSagitta: opts.arcSagitta ?? Math.max(0.1, 0.00015 * edge), + // How far the contour's local direction may stray from a fitted line + // before that stretch stops counting as part of the straight edge. 4 + // degrees: below it the difference is inside the tangent estimate's own + // noise, above it the bend is already under way and belongs to the curve. + lineTangentTol: ((opts.lineTangentDeg ?? 3) * Math.PI) / 180, + // Arc length of the window whose fitted direction is compared to the + // run. Long enough that its own angular noise is well under the + // tolerance, short enough to locate where the bend actually starts. + trimWin: opts.trimWin ?? Math.max(6, 0.0067 * edge), + // The radius the blur rounds a corner to, in source px. Does two jobs: + // it is kept out of the line fit at a corner so the rounding cannot tilt + // the direction, and it sets how far past the contour the reconstructed + // apex may legitimately sit for a given wedge angle. + cornerGuard: opts.cornerGuard ?? Math.max(1.5, 0.00125 * edge), + // Hard ceiling on how far a reconstructed corner may sit from the + // contour point it replaces, on top of the angle-derived bound. + // + // It exists for the case the angle bound cannot cover: as the wedge + // closes, the legitimate apex offset grows without limit, and so does + // the error in it, because two nearly anti-parallel lines cross at a + // point that a fraction of a degree can move a long way. At 12px a + // narrow "V" -- the inside of Mildlife's M -- reconstructed 11px past + // its real apex and cut a white spike into the letterform. Measured over + // the acceptance set, corner sharpness is IDENTICAL anywhere from 2px to + // 12px (Hot Toddy 0.25px turn-arc, 0 of 34 softened, at every value), so + // this ceiling costs nothing at all and only removes the failure. + cornerReach: opts.cornerReach ?? Math.max(3, 0.002 * edge), + cornerWin, + tanWin: opts.tanWin ?? cornerWin * 0.7, + // 38 degrees. Lower starts finding "corners" on smooth bowls, and a + // false corner is a visible kink -- worse than the facets being removed. + minTurn: ((opts.minTurnDeg ?? 38) * Math.PI) / 180, + // The same question asked of two FITTED LINES, where the evidence is far + // stronger, so the bar is far lower: below this the two are one curve. + lineCornerTurn: ((opts.lineCornerTurnDeg ?? 16) * Math.PI) / 180, + minPerim: opts.minPerim ?? Math.max(8, 0.005 * edge), + eps: opts.eps ?? 0.22, + // Acceptance is on RMS; `fitMax` is only an outlier ceiling. Splitting a + // cubic on its worst single sample makes the segment count a function of + // the grain rather than of the shape. + // tol/2 reproduces what a max-residual test at `tol` used to accept on + // these rasters, but without inheriting its dependence on how many + // samples the run happens to contain. + fitRms: opts.fitRms ?? tol / 2, + fitMax: opts.fitMax ?? tol * 2, + // Past three passes Newton has converged; the remaining segment-count + // gain is a couple of percent for a third more time. + // Beyond this the fit is too far out for Newton to rescue and splitting + // is both cheaper and more likely to be what it needs. + reparamGate: opts.reparamGate ?? 3, + reparamIters: opts.reparamIters ?? 3, + maxDepth: opts.maxDepth ?? 24, + }; +} + /** Signed area -- used to discard specks by size rather than point count. */ function area(pts) { let s = 0; @@ -293,18 +1305,21 @@ export function traceToSvg(px, w, h, opts = {}) { const { eps = 0.22, // simplify tolerance in source px; sub-pixel, so keep small minArea = 6, // drop specks (source px^2) + fit = true, // fit curves; false emits the old RDP polyline } = opts; - const loops = contours(inkField(px, w, h), w, h) - .filter((l) => area(l) >= minArea) - .map((l) => simplify(l, eps)); - if (!loops.length) return null; + const raw = contours(inkField(px, w, h), w, h).filter((l) => area(l) >= minArea); + if (!raw.length) return null; + // Everything below works in the ink's own box, so measure it on the source + // points -- fitting only ever moves a vertex by a fraction of the tolerance, + // and taking the box from the fitted output would let the tolerance jitter + // the viewBox between runs. let minX = Infinity; let minY = Infinity; let maxX = -Infinity; let maxY = -Infinity; - for (const l of loops) + for (const l of raw) for (const [x, y] of l) { if (x < minX) minX = x; if (x > maxX) maxX = x; @@ -316,12 +1331,33 @@ export function traceToSvg(px, w, h, opts = {}) { if (!(bw > 0 && bh > 0)) return null; const n = (v) => Math.round(v * 100) / 100; - const d = loops - .map( - (l) => - "M" + l.map(([x, y]) => `${n(x - minX)} ${n(y - minY)}`).join("L") + "Z", - ) - .join(""); + let d; + if (fit) { + const o = fitOptions({ ...opts, w, h }); + d = raw + .map((l) => fitLoop(l, o)) + .filter(Boolean) + .map(({ start, segs }) => { + let s = `M${n(start[0] - minX)} ${n(start[1] - minY)}`; + for (const sg of segs) { + const pt = ([x, y]) => `${n(x - minX)} ${n(y - minY)}`; + s += sg.length === 3 + ? `C${pt(sg[0])} ${pt(sg[1])} ${pt(sg[2])}` + : `L${pt(sg[0])}`; + } + return s + "Z"; + }) + .join(""); + if (!d) return null; + } else { + d = raw + .map((l) => simplify(l, eps)) + .map( + (l) => + "M" + l.map(([x, y]) => `${n(x - minX)} ${n(y - minY)}`).join("L") + "Z", + ) + .join(""); + } // Filled black: a CSS mask consumes the SVG's alpha, so only coverage // matters. evenodd keeps letter counters (the hole in an "o") open.