31:T4b3,function fitMobius(pairs) { if (pairs.length < 3) throw new Error("need >=3 pairs"); const a = pairs[0], b = pairs[1], c = pairs[2]; // guard degeneracies if (a.x === b.x || a.x === c.x || b.x === c.x) throw new Error("x must be distinct"); if (b.k === a.k || b.k === c.k) throw new Error("k must satisfy B≠A and B≠C"); const A = a.k, B = b.k, C = c.k; const alpha = (B - C) / (B - A); // M(z) = CR(z; a.x,b.x,c.x) = ((b-c)z - (b-c)a) / ((b-a)z - (b-a)c) const m11 = b.x - c.x; const m12 = -a.x * (b.x - c.x); const m21 = b.x - a.x; const m22 = -c.x * (b.x - a.x); // L(w) = (C w - alpha A) / (w - alpha) maps {0,1,∞} -> {A,B,C} const l11 = C, l12 = -alpha * A, l21 = 1, l22 = -alpha; // Composition T = L ∘ M ⇒ [p q; r s] = L · M (up to scale) let p = l11 * m11 + l12 * m21; let q = l11 * m12 + l12 * m22; let r = l21 * m11 + l22 * m21; let s = l21 * m12 + l22 * m22; // normalize scale if (Math.abs(s) > 1e-12) { p /= s; q /= s; r /= s; s = 1; } else { const norm = Math.hypot(p, q, r, s) || 1; p /= norm; q /= norm; r /= norm; s /= norm; } return { p, q, r, s }; }32:T1120,dpAlign = function dpAlign(template, scan, opts = {}) { // Needleman-Wunsch style dynamic programming alignment for 1D positions. // template: array of numbers or objects with .x and optional .s (polarity) // scan: array of numbers or objects with .x and optional .s (polarity) // opts: // gapPenalty (number) default 20 // polarityPenalty (number) default 50 // matchCost (fn) optional custom cost fn (xt, xs) -> number; defaults to abs difference // matchThreshold (number) if provided, matches with cost > threshold are discouraged (treated normally by DP) const { gapPenalty = 20, polarityPenalty = 50, matchCost = (a, b) => Math.abs(a - b), matchThreshold = Infinity } = opts; const N = template.length, M = scan.length; // helpers to extract x and polarity const getX = (v) => typeof v === "number" ? v : v && ("x" in v ? v.x : v.x === 0 ? 0 : v); const getS = (v) => typeof v === "number" ? v : v && "s" in v ? v.s : v && "sign" in v ? v.sign : 0; // DP matrices: use flat arrays for speed const rows = N + 1, cols = M + 1; const D = new Float64Array(rows * cols); // cost matrix const PTR = new Int8Array(rows * cols); // 0 = diag, 1 = up (delete template / gap in scan), 2 = left (insert / gap in template) const idx = (i, j) => i * cols + j; // initialize for (let i = 1; i <= N; i++) { D[idx(i, 0)] = i * gapPenalty; PTR[idx(i, 0)] = 1; // from up } for (let j = 1; j <= M; j++) { D[idx(0, j)] = j * gapPenalty; PTR[idx(0, j)] = 2; // from left } // fill for (let i = 1; i <= N; i++) { const xt = getX(template[i - 1]); const st = getS(template[i - 1]); for (let j = 1; j <= M; j++) { const xs = getX(scan[j - 1]); const ss = getS(scan[j - 1]); let mcost = matchCost(xt, xs); // polarity penalty (if either side has a defined sign and they differ) if (st && ss && Math.sign(st) !== Math.sign(ss)) mcost += polarityPenalty; // match threshold can be used by caller to interpret results (we don't hard forbid >threshold here) const diag = D[idx(i - 1, j - 1)] + mcost; const up = D[idx(i - 1, j)] + gapPenalty; // skip template[i-1] const left = D[idx(i, j - 1)] + gapPenalty; // skip scan[j-1] // choose minimum (prefer diag on ties, then up) let best = diag, ptr = 0; if (up < best || (up === best && ptr !== 0 && ptr !== 1)) { best = up; ptr = 1; } if (left < best) { best = left; ptr = 2; } D[idx(i, j)] = best; PTR[idx(i, j)] = ptr; } } // backtrace let i = N, j = M; const pairs = []; // aligned pairs (ti, sj, cost) while (i > 0 || j > 0) { const p = PTR[idx(i, j)]; if (p === 0) { // diagonal => match template[i-1] with scan[j-1] const xt = getX(template[i - 1]), xs = getX(scan[j - 1]); let c = matchCost(xt, xs); const st = getS(template[i - 1]), ss = getS(scan[j - 1]); if (st && ss && Math.sign(st) !== Math.sign(ss)) c += polarityPenalty; pairs.push({ ti: i - 1, sj: j - 1, cost: c, xt, xs, st, ss }); i--; j--; } else if (p === 1) { // up => gap in scan (template element skipped) pairs.push({ ti: i - 1, sj: null, cost: gapPenalty, xt: getX(template[i - 1]), xs: null, st: getS(template[i - 1]), ss: null }); i--; } else { // left => gap in template (scan element skipped) pairs.push({ ti: null, sj: j - 1, cost: gapPenalty, xt: null, xs: getX(scan[j - 1]), st: null, ss: getS(scan[j - 1]) }); j--; } } pairs.reverse(); // Build template->scan mapping (one-to-one or null). If multiple matches to same scan occur (shouldn't for NW), // prefer the diagonal matches. const templateToScan = new Array(N).fill(null); const scanToTemplate = new Array(M).fill(null); for (const p of pairs) { if (p.ti != null && p.sj != null) { templateToScan[p.ti] = p.sj; scanToTemplate[p.sj] = p.ti; } } return { score: D[idx(N, M)], pairs, templateToScan, scanToTemplate, // expose DP dimensions so caller can inspect if desired N, M }; }33:T9c1,computeAnchorCandidatesFast = function computeAnchorCandidatesFast( template_edges, scanEdges, opts = {} ) { const kSorted = template_edges.slice().sort((a, b) => a - b); if (kSorted.length < 4) throw new Error("need >=4 template edges"); const kMin = kSorted[0], kMax = kSorted[kSorted.length - 1], kLeftInner = kSorted[1], kRightInner = kSorted[kSorted.length - 2]; const sLen = scanEdges.length; if (sLen < 4) return []; const innerExpected = kSorted.length - 2; const minBetween = Math.max(2, opts.minBetween ?? innerExpected - 1); const maxBetween = Math.min(sLen - 2, opts.maxBetween ?? innerExpected + 1); const topK = opts.topK ?? 10; const scanX = scanEdges.map((e) => (typeof e === "number" ? e : e.x)); const candidates = []; for (let i = 0; i < sLen - 1; i++) { const jMin = i + 1 + minBetween; const jMax = Math.min(sLen - 1, i + 1 + maxBetween); if (jMin >= sLen) continue; for (let j = jMin; j <= jMax; j++) { if (i + 1 >= j) continue; const pairs4 = [ { x: scanX[i], k: kMin }, { x: scanX[i + 1], k: kLeftInner }, { x: scanX[j - 1], k: kRightInner }, { x: scanX[j], k: kMax } ]; let mob; try { mob = fitMobiusLS(pairs4); if ( !( isFinite(mob.p) && isFinite(mob.q) && isFinite(mob.r) && isFinite(mob.s) ) ) continue; } catch { continue; } const kErrs = new Array(pairs4.length); let bad = false; for (let t = 0; t < pairs4.length; t++) { const x = pairs4[t].x, k = pairs4[t].k; const den = mob.r * x + mob.s; if (Math.abs(den) < 1e-12) { bad = true; kErrs[t] = Infinity; } else { const khat = (mob.p * x + mob.q) / den; kErrs[t] = khat - k; } } const finiteErrs = kErrs.filter((v) => isFinite(v)); const kRMSE = finiteErrs.length ? Math.sqrt( finiteErrs.reduce((a, b) => a + b * b, 0) / finiteErrs.length ) : Infinity; candidates.push({ startIndex: i, endIndex: j, pairs4, mobiusInitial: mob, kRMSE }); } } candidates.sort((a, b) => { const ka = a.kRMSE === undefined ? Infinity : a.kRMSE; const kb = b.kRMSE === undefined ? Infinity : b.kRMSE; return ka - kb; }); return candidates.slice(0, topK); }34:T948,function refineAnchorCandidates( candidates, template_edges, scanEdges, opts = {} ) { const topK = opts.topK ?? candidates.length; if (!Array.isArray(candidates) || candidates.length === 0) return []; const kSorted = template_edges.slice().sort((a, b) => a - b); const scanX = scanEdges.map((e) => (typeof e === "number" ? e : e.x)); const dpOpts = { gapPenalty: opts.gapPenalty ?? 20, polarityPenalty: opts.polarityPenalty ?? 0, matchCost: (a, b) => Math.abs(a - b), matchThreshold: Infinity }; const out = []; for (const cand of candidates) { const mobInit = cand.mobiusInitial ?? cand.mobius; if (!mobInit) continue; const proj = kSorted.map((k) => { const x = xFromK(mobInit, k); return x; }); const alignment = dpAlign(proj, scanX, dpOpts); const map = alignment.templateToScan || []; const matchedPairs = []; for (let t = 0; t < kSorted.length; t++) { const sIdx = map[t]; if (sIdx != null) matchedPairs.push({ x: scanX[sIdx], k: kSorted[t] }); } if (matchedPairs.length < 3) continue; let mobRefined; try { mobRefined = fitMobiusLS(matchedPairs); } catch { continue; } const { p, q, r, s } = mobRefined; const kErrs = []; const xErrs = []; for (let t = 0; t < kSorted.length; t++) { const sIdx = map[t]; if (sIdx == null) continue; const x = scanX[sIdx]; const kTrue = kSorted[t]; const den = r * x + s; if (!isFinite(den) || Math.abs(den) < 1e-12) continue; const khat = (p * x + q) / den; kErrs.push(khat - kTrue); const predX = xFromK(mobRefined, kTrue); if (isFinite(predX)) xErrs.push(predX - x); } if (kErrs.length === 0) continue; const kRMSE = Math.sqrt( kErrs.reduce((a, b) => a + b * b, 0) / kErrs.length ); const xRMSE = xErrs.length ? Math.sqrt(xErrs.reduce((a, b) => a + b * b, 0) / xErrs.length) : Infinity; out.push({ startIndex: cand.startIndex, endIndex: cand.endIndex, pairs4: cand.pairs4, mobiusInitial: cand.mobiusInitial, mobiusRefined: mobRefined, kRMSE_full: kRMSE, xRMSE_full: xRMSE, pairsUsed: matchedPairs.length }); } out.sort( (a, b) => a.kRMSE_full - b.kRMSE_full || a.xRMSE_full - b.xRMSE_full ); return out.slice(0, topK); }35:Ta0b,autoFitScanline = function autoFitScanline(scanLine, template, opts = {}) { if (!scanLine.length) throw new Error("scanLine must be an array"); if (!Array.isArray(template)) throw new Error( "template must be an array (either binary mask or template k-values)" ); const edgeThr = opts.edgeThreshold ?? 6; const topK = opts.topK ?? 10; const topKRefine = opts.topKRefine ?? 5; const minBetween = opts.minBetween; const maxBetween = opts.maxBetween; const gapPenalty = opts.gapPenalty; const polarityPenalty = opts.polarityPenalty; let template_edges; if (opts.templateIsBinary === true) { template_edges = binaryToEdges(template); } else { const allZerosOnes = template.every((v) => v === 0 || v === 1); if (allZerosOnes) template_edges = binaryToEdges(template); else template_edges = template.slice(); } const scanEdges = edges1D(scanLine, edgeThr); if (!scanEdges || scanEdges.length < 3) return { success: false, reason: "not_enough_scan_edges", scanEdges, template_edges }; const candidates = computeAnchorCandidatesFast(template_edges, scanEdges, { topK, minBetween, maxBetween }); const refined = refineAnchorCandidates( candidates, template_edges, scanEdges, { topK: topKRefine, gapPenalty, polarityPenalty } ); const chosen = refined && refined.length ? refined[0] : candidates && candidates.length ? candidates[0] : null; const mob = chosen ? chosen.mobiusRefined ?? chosen.mobiusInitial ?? chosen.mobius : null; const projectedTemplate = mob ? template_edges.map((k) => ({ k, x: xFromK(mob, k) })) : []; const tplX = projectedTemplate.map((t) => t.x); const scanX = scanEdges.map((e) => (typeof e === "number" ? e : e.x)); const alignment = dpAlign(tplX, scanX); const map = alignment.templateToScan || []; const residuals = projectedTemplate.map((t, i) => { const sj = map[i]; if (sj != null) { const actual = scanX[sj]; return { k: t.k, predictedX: t.x, actualX: actual, cost: Math.abs(t.x - actual), matched: true, sj }; } else { return { k: t.k, predictedX: t.x, actualX: null, cost: 20, matched: false, sj: null }; } }); return { success: true, template_edges, scanEdges, candidates, refined, chosen, mobius: mob, projectedTemplate, alignment, residuals }; }36:T1080,function canny(src, w, opts = {}) { const { sigma = 1.0, low = 20, high = 60 } = opts; if (!Number.isInteger(w) || w <= 0) throw new Error("w must be positive integer"); const h = (src.length / w) | 0; if (w * h !== src.length) throw new Error("buffer length not divisible by w"); // 1) Gaussian blur, separable const k = gaussianKernel1D(sigma); const tmp = new Float32Array(w * h); const blur = new Float32Array(w * h); // horizontal for (let y = 0; y < h; y++) { const base = y * w; for (let x = 0; x < w; x++) { let acc = 0; for (let i = -k.half; i <= k.half; i++) { const xx = clamp(x + i, 0, w - 1); acc += k.data[i + k.half] * src[base + xx]; } tmp[base + x] = acc; } } // vertical for (let x = 0; x < w; x++) { for (let y = 0; y < h; y++) { let acc = 0; for (let i = -k.half; i <= k.half; i++) { const yy = clamp(y + i, 0, h - 1); acc += k.data[i + k.half] * tmp[yy * w + x]; } blur[y * w + x] = acc; } } // 2) Sobel gradients const gx = new Float32Array(w * h); const gy = new Float32Array(w * h); for (let y = 0; y < h; y++) { const ym1 = Math.max(0, y - 1), yp1 = Math.min(h - 1, y + 1); for (let x = 0; x < w; x++) { const xm1 = Math.max(0, x - 1), xp1 = Math.min(w - 1, x + 1); const a = blur[ym1 * w + xm1], b = blur[ym1 * w + x], c = blur[ym1 * w + xp1]; const d = blur[y * w + xm1], /* e = blur[y*w + x] */ f = blur[y * w + xp1]; const g = blur[yp1 * w + xm1], hh = blur[yp1 * w + x], i = blur[yp1 * w + xp1]; gx[y * w + x] = c + 2 * f + i - (a + 2 * d + g); gy[y * w + x] = g + 2 * hh + i - (a + 2 * b + c); } } // gradient magnitude (L2), scaled to 0..255 for thresholding const mag = new Float32Array(w * h); let maxMag = 0; for (let idx = 0; idx < mag.length; idx++) { const m = Math.hypot(gx[idx], gy[idx]); mag[idx] = m; if (m > maxMag) maxMag = m; } const scale = maxMag > 0 ? 255 / maxMag : 0; // 3) Non-maximum suppression along quantized directions (0,45,90,135) const thin = new Uint8Array(w * h); // holds scaled magnitude at local maxima for (let y = 1; y < h - 1; y++) { for (let x = 1; x < w - 1; x++) { const idx = y * w + x; const gxx = gx[idx], gyy = gy[idx]; const m = mag[idx] * scale; // direction sector const angle = Math.atan2(gyy, gxx) * (180 / Math.PI); const a = angle < 0 ? angle + 180 : angle; let n1 = 0, n2 = 0; if ((a >= 0 && a < 22.5) || (a >= 157.5 && a < 180)) { n1 = mag[idx - 1] * scale; n2 = mag[idx + 1] * scale; // horizontal } else if (a >= 22.5 && a < 67.5) { n1 = mag[idx - w - 1] * scale; n2 = mag[idx + w + 1] * scale; // 45° } else if (a >= 67.5 && a < 112.5) { n1 = mag[idx - w] * scale; n2 = mag[idx + w] * scale; // vertical } else { // 112.5..157.5 n1 = mag[idx - w + 1] * scale; n2 = mag[idx + w - 1] * scale; // 135° } thin[idx] = m >= n1 && m >= n2 ? m | 0 : 0; } } // 4) Hysteresis thresholding via stack flood fill const STRONG = 255; const WEAK = 128; const out = new Uint8Array(w * h); const stack = new Int32Array(w * h); let sp = 0; for (let i = 0; i < thin.length; i++) { const v = thin[i]; if (v >= high) { out[i] = STRONG; stack[sp++] = i; } else if (v >= low) { out[i] = WEAK; } } // promote weak connected to strong (8-connectivity) while (sp > 0) { const idx = stack[--sp]; const y = (idx / w) | 0, x = idx % w; for (let dy = -1; dy <= 1; dy++) { for (let dx = -1; dx <= 1; dx++) { if (dx === 0 && dy === 0) continue; const xx = x + dx, yy = y + dy; if (xx < 0 || xx >= w || yy < 0 || yy >= h) continue; const j = yy * w + xx; if (out[j] === WEAK) { out[j] = STRONG; stack[sp++] = j; } } } } // suppress remaining weak for (let i = 0; i < out.length; i++) out[i] = out[i] === STRONG ? 255 : 0; return out; }37:[{"id":"761b4331b51bb160fe71de3b320d5fd4b3afab5f04b079862051dc2a60264ad23e154d0e2deb3a18f9c7e982b5e2a8bc0837599de53e888fe7744a23a842b0d7","href":"https://static.observableusercontent.com/files/761b4331b51bb160fe71de3b320d5fd4b3afab5f04b079862051dc2a60264ad23e154d0e2deb3a18f9c7e982b5e2a8bc0837599de53e888fe7744a23a842b0d7","name":"image.png","type":"image/png","lastModified":1760541582576.459,"size":47223}] 2f:{"title":"Fast 1D Circular Barcode Matching","cells":[{"id":0,"value":"# Fast 1D Circular Barcode Matching\n## with the Möbius transform\n\n\nThe aim of this notebook is to prototype a way of quickly detecting the center and tilt of a circular barcode passing through a scan line. This will be a building block of a localization algorithm.\n\nPart II of a long quest to create a cheap high speed optical localization system. Part I is [here](https://observablehq.com/@tomlarkworthy/circular-barcode-simulator).","pinned":false,"hidden":false,"mode":"md"},{"id":4,"value":"import {\n viewof scanY,\n canvas,\n pixelBuffer,\n renders,\n scanline,\n overlay,\n template\n} from \"@tomlarkworthy/circular-barcode-simulator\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":638,"value":"### Final Matched Barcode projection\nIt works! If the scanline passed through the centre of a circular barcode, this matches the pattern.","pinned":false,"hidden":false,"mode":"md"},{"id":632,"value":"Plot.plot({\n width,\n height: 100,\n y: { axis: false },\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(autoFitProjectedTemplate, {\n x: (d) => d.x\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":364,"value":"viewof scanY","pinned":false,"hidden":false,"mode":"ojs"},{"id":661,"value":"

\n ${canvas}\n

","pinned":false,"hidden":false,"mode":"html"},{"id":655,"value":"## Implementation","pinned":false,"hidden":false,"mode":"md"},{"id":517,"value":"#### Pixel intensity through scan line","pinned":false,"hidden":false,"mode":"md"},{"id":114,"value":"viewof scanLinePlot = Plot.plot({\n width,\n marks: [Plot.tickX(scanLineBuffer, { x: (d, i) => i, stroke: (d) => d })]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":102,"value":"scanLineBuffer = new Uint8Array(scanline.map((p) => p.v))","pinned":false,"hidden":false,"mode":"ojs"},{"id":122,"value":"scanEdges = edges1D(scanLineBuffer, 20)","pinned":false,"hidden":false,"mode":"ojs"},{"id":423,"value":"#### (manual) Find three known edges","pinned":false,"hidden":false,"mode":"md"},{"id":327,"value":"In the three charts below, click the leftmost side of the a barcode's circular pattern, the rightmost side and then the right of the middle as the three rings to fit mobius to. The larger point of this notebook is to figure out an efficient way to initialize the algorithm *without* manual initialization, but to get there it is useful to first concentrate on the easier problem of an already initialized fit.\n\n

example

\n ![image.png](${await FileAttachment(\"image.png\").url()})\n

","pinned":false,"hidden":false,"mode":"md"},{"id":165,"value":"viewof leftmost = Plot.plot({\n width,\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(scanEdges, { x: (d) => d.x, stroke: (d) => d.s, tip: true }),\n Plot.dot(scanEdges, Plot.pointer({ x: (d) => d.x }))\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":259,"value":"viewof rightmost = Plot.plot({\n width,\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(scanEdges, { x: (d) => d.x, stroke: (d) => d.s, tip: true }),\n Plot.dot(scanEdges, Plot.pointer({ x: (d) => d.x }))\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":261,"value":"viewof middle_right = Plot.plot({\n width,\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(scanEdges, { x: (d) => d.x, stroke: (d) => d.s, tip: true }),\n Plot.dot(scanEdges, Plot.pointer({ x: (d) => d.x }))\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":421,"value":"#### Initial Mobius Fit to 3 anchors","pinned":false,"hidden":false,"mode":"md"},{"id":352,"value":"Predicted template positions given pqrs fit to anchors.","pinned":false,"hidden":false,"mode":"md"},{"id":302,"value":"Plot.plot({\n width,\n height: 100,\n y: { axis: false },\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(template_projection, {\n x: (d) => d.x\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":465,"value":"##### Residuals on initial template fit","pinned":false,"hidden":false,"mode":"md"},{"id":452,"value":"Plot.plot({\n width,\n height: 100,\n color: {\n domain: [0, 1]\n },\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.ruleX(initialAlignment.pairs, {\n x: (d) => d.xt,\n y: (d) => (d.cost == 20 ? 0 : d.cost),\n stroke: (d) => (d.cost == 20 ? 0 : d.cost)\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":414,"value":"### Least squares Möbius fit to template match","pinned":false,"hidden":false,"mode":"md"},{"id":429,"value":"Plot.plot({\n width,\n height: 100,\n y: { axis: false },\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.tickX(mobiusLSFit.fittedProjection, {\n x: (d) => d.x\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":485,"value":"##### Residuals on LS template fit: kRMSE ${mobiusLSFit.kRMSE}","pinned":false,"hidden":false,"mode":"md"},{"id":473,"value":"Plot.plot({\n width,\n height: 100,\n color: {\n domain: [0, 1]\n },\n x: {\n domain: [0, width]\n },\n marks: [\n Plot.ruleX(residualsLS, {\n x: (d) => d.x,\n y: (d) => (d.cost == 20 ? 0 : d.cost),\n stroke: (d) => (d.cost == 20 ? 0 : d.cost)\n })\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":512,"value":"# Code","pinned":false,"hidden":false,"mode":"md"},{"id":345,"value":"### fit a Möbius transform to three points","pinned":false,"hidden":false,"mode":"md"},{"id":125,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":254,"value":"pqrs = fitMobius([\n { ...leftmost, k: -27 },\n { ...middle_right, k: 5 },\n { ...rightmost, k: 27 }\n])","pinned":false,"hidden":false,"mode":"ojs"},{"id":442,"value":"### Predict template edges given fit","pinned":false,"hidden":false,"mode":"md"},{"id":332,"value":"template_edges = binaryToEdges(template)","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"function xFromK(pqrs, k) {\n const denom = k * pqrs.r - pqrs.p;\n if (Math.abs(denom) < 1e-12) return NaN; // parallel / undefined\n return (pqrs.q - k * pqrs.s) / denom;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":290,"value":"template_projection = template_edges.map((k) => ({ k, x: xFromK(pqrs, k) }))","pinned":false,"hidden":false,"mode":"ojs"},{"id":411,"value":"### Obtain closest matching edges to initial Möbius estimate","pinned":false,"hidden":false,"mode":"md"},{"id":386,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":388,"value":"initialAlignment = dpAlign(\n template_projection.map((t) => t.x),\n scanEdges.map((s) => s.x)\n)","pinned":false,"hidden":false,"mode":"ojs"},{"id":418,"value":"### Least squares mobius fit to paired edges","pinned":false,"hidden":false,"mode":"md"},{"id":354,"value":"function fitMobiusLS(pairs) {\n const N = pairs.length;\n if (N < 3) throw new Error(\"need >=3 points\");\n\n // Build A as N×4 array\n const A = new Array(N);\n for (let i = 0; i < N; i++) {\n const { x, k } = pairs[i];\n A[i] = [x, 1, -k * x, -k];\n }\n\n // SVD: A = U Σ Vᵀ\n const { v } = SVD(A);\n\n // v is 4×4; columns are right singular vectors\n // The column corresponding to the smallest singular value is the nullspace vector\n const smallestColIndex = 3; // svd-js returns singular values in descending order\n let p = v[0][smallestColIndex];\n let q = v[1][smallestColIndex];\n let r = v[2][smallestColIndex];\n let s = v[3][smallestColIndex];\n\n // normalize so s = 1\n if (Math.abs(s) > 1e-12) {\n p /= s;\n q /= s;\n r /= s;\n s = 1;\n }\n\n return { p, q, r, s };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":406,"value":"mobiusLSFit = {\n const matches = initialAlignment.pairs.filter(\n (p) => p.ti != null && p.sj != null\n );\n const pairs = matches.map((p) => ({\n x: scanEdges[p.sj].x,\n k: template_projection[p.ti].k\n }));\n if (pairs.length < 3) throw new Error(\"need >=3 matched pairs for mobiusLS\");\n const mobiusLS = fitMobiusLS(pairs);\n // compute residuals in k-space (mapped k from x) and attach diagnostics\n const { p, q, r, s } = mobiusLS;\n const kErrs = pairs.map(({ x, k }) => {\n const khat = (p * x + q) / (r * x + s);\n return khat - k;\n });\n const kRMSE = Math.sqrt(kErrs.reduce((a, b) => a + b * b, 0) / kErrs.length);\n const fittedProjection = template_projection.map((t) => ({\n k: t.k,\n x: xFromK(mobiusLS, t.k)\n }));\n mobiusLS.kRMSE = kRMSE;\n mobiusLS.pairsUsed = pairs.length;\n mobiusLS.fittedProjection = fittedProjection;\n return mobiusLS;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":476,"value":"function calculatLSResiduals(fitted) {\n const N = fitted.length;\n const map = initialAlignment.templateToScan || [];\n const out = new Array(N);\n for (let i = 0; i < N; i++) {\n const tplX = template_projection[i].x;\n const predX = fitted[i] ? fitted[i].x : NaN;\n const sj = map[i];\n if (sj != null) {\n const actualX = scanEdges[sj].x;\n out[i] = {\n x: tplX,\n cost: Math.abs(predX - actualX),\n matched: true,\n predictedX: predX,\n actualX\n };\n } else {\n out[i] = {\n x: tplX,\n cost: 20,\n matched: false,\n predictedX: predX,\n actualX: null\n };\n }\n }\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":495,"value":"residualsLS = calculatLSResiduals(mobiusLSFit.fittedProjection)","pinned":false,"hidden":false,"mode":"ojs"},{"id":520,"value":"## Auto fit anchors","pinned":false,"hidden":false,"mode":"md"},{"id":571,"value":"$33","pinned":true,"hidden":false,"mode":"ojs"},{"id":575,"value":"anchorCandidatesFast = computeAnchorCandidatesFast(template_edges, scanEdges, {\n topK: 20\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":608,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":612,"value":"autoFit = refineAnchorCandidates(\n anchorCandidatesFast,\n template_edges,\n scanEdges,\n {\n topK: 1\n }\n)[0]","pinned":false,"hidden":false,"mode":"ojs"},{"id":620,"value":"autoFitProjectedTemplate = {\n const mob =\n autoFit &&\n (autoFit.mobiusRefined ?? autoFit.mobiusInitial ?? autoFit.mobius);\n if (!mob) return [];\n return template_edges.map((k) => ({ k, x: xFromK(mob, k) }));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":622,"value":"autoFitResiduals = {\n const projectedTemplate = autoFitProjectedTemplate;\n if (!projectedTemplate || projectedTemplate.length === 0) return [];\n const tplX = projectedTemplate.map((t) => t.x);\n const scanX = scanEdges.map((e) => (typeof e === \"number\" ? e : e.x));\n const alignment = dpAlign(tplX, scanX);\n const map = alignment.templateToScan || [];\n const out = new Array(projectedTemplate.length);\n for (let i = 0; i < projectedTemplate.length; i++) {\n const pred = projectedTemplate[i].x;\n const sj = map[i];\n if (sj != null) {\n const actual = scanX[sj];\n out[i] = {\n k: projectedTemplate[i].k,\n predictedX: pred,\n actualX: actual,\n cost: Math.abs(pred - actual),\n matched: true,\n sj\n };\n } else {\n out[i] = {\n k: projectedTemplate[i].k,\n predictedX: pred,\n actualX: null,\n cost: 20,\n matched: false,\n sj: null\n };\n }\n }\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":678,"value":"## End to end function","pinned":false,"hidden":false,"mode":"md"},{"id":683,"value":"$35","pinned":true,"hidden":false,"mode":"ojs"},{"id":685,"value":"autoFitScanline(scanLineBuffer, template)","pinned":true,"hidden":false,"mode":"ojs"},{"id":561,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":335,"value":"function binaryToEdges(arr) {\n const edges = [];\n const n = arr.length;\n\n // transitions\n for (let i = 1; i < n; i++) {\n if (arr[i] !== arr[i - 1]) edges.push(i);\n }\n\n return edges.map((i) => i - arr.length / 2);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"grey = {\n renders;\n return rgbaToGray(pixelBuffer);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":41,"value":"edges = {\n renders;\n return canny(grey, width);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":133,"value":"test_fitMobius = {\n // Ground truth params\n const p = 0.5,\n q = 2,\n r = 0.05,\n s = 1;\n\n // Make 4 points from the model\n const pts = [0, 10, 20, 30].map((x) => ({ x, k: (p * x + q) / (r * x + s) }));\n\n return fitMobius(pts);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":127,"value":"function solve3(A, b) {\n const [a, b1, c, d, e, f, g, h, i] = [\n A[0][0],\n A[0][1],\n A[0][2],\n A[1][0],\n A[1][1],\n A[1][2],\n A[2][0],\n A[2][1],\n A[2][2]\n ];\n const D = a * (e * i - f * h) - b1 * (d * i - f * g) + c * (d * h - e * g);\n if (Math.abs(D) < 1e-8) return [0, 0, 0];\n const dx =\n (b[0] * (e * i - f * h) -\n b1 * (b[1] * i - f * b[2]) +\n c * (b[1] * h - e * b[2])) /\n D;\n const dy =\n (a * (b[1] * i - f * b[2]) -\n b[0] * (d * i - f * g) +\n c * (d * b[2] - b[1] * g)) /\n D;\n const dz =\n (a * (e * b[2] - b[1] * h) -\n b1 * (d * b[2] - b[1] * g) +\n b[0] * (d * h - e * g)) /\n D;\n return [dx, dy, dz];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":78,"value":"function sampleLine(gray, w, cx, cy, theta, len) {\n const h = (gray.length / w) | 0;\n const dx = Math.cos(theta),\n dy = Math.sin(theta);\n const half = (len / 2) | 0;\n const out = new Float32Array(len);\n for (let t = -half, j = 0; j < len; t++, j++) {\n const xf = cx + t * dx,\n yf = cy + t * dy;\n const x0 = Math.floor(xf),\n y0 = Math.floor(yf);\n if (x0 < 0 || x0 >= w - 1 || y0 < 0 || y0 >= h - 1) {\n out[j] = 0;\n continue;\n }\n const a = xf - x0,\n b = yf - y0;\n const i = y0 * w + x0;\n const p00 = gray[i],\n p10 = gray[i + 1],\n p01 = gray[i + w],\n p11 = gray[i + w + 1];\n out[j] =\n (1 - a) * (1 - b) * p00 +\n a * (1 - b) * p10 +\n (1 - a) * b * p01 +\n a * b * p11;\n }\n return out;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"function rgbaToGray(rgba) {\n if (rgba.length % 4 !== 0) throw new Error(\"RGBA length not multiple of 4\");\n const gray = new Uint8Array(rgba.length / 4);\n for (let i = 0, j = 0; i < rgba.length; i += 4, j++) {\n const r = rgba[i],\n g = rgba[i + 1],\n b = rgba[i + 2];\n gray[j] = (0.299 * r + 0.587 * g + 0.114 * b) | 0;\n }\n return gray;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":111,"value":"function edges1D(sig, thr = 6) {\n const n = sig.length;\n const d = new Float32Array(n);\n for (let i = 1; i < n; i++) d[i] = sig[i] - sig[i - 1];\n const idx = [];\n for (let i = 2; i < n - 2; i++) {\n const v = d[i];\n if (Math.abs(v) < thr) continue;\n if (\n (v > 0 && d[i] >= d[i - 1] && d[i] >= d[i + 1]) ||\n (v < 0 && d[i] <= d[i - 1] && d[i] <= d[i + 1])\n ) {\n idx.push({ x: i, s: Math.sign(v) });\n }\n }\n return idx;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":45,"value":"function canvasFromGrayU8(gray, w) {\n if (!Number.isInteger(w) || w <= 0) throw new Error(\"w invalid\");\n const h = (gray.length / w) | 0;\n if (w * h !== gray.length) throw new Error(\"buffer size mismatch\");\n\n const canvas = document.createElement(\"canvas\");\n canvas.width = w;\n canvas.height = h;\n const ctx = canvas.getContext(\"2d\");\n const img = ctx.createImageData(w, h);\n const dst = img.data; // RGBA\n\n for (let i = 0, j = 0; i < gray.length; i++, j += 4) {\n const v = gray[i];\n dst[j] = v;\n dst[j + 1] = v;\n dst[j + 2] = v;\n dst[j + 3] = 255;\n }\n ctx.putImageData(img, 0, 0);\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"function gaussianKernel1D(sigma) {\n const s = Math.max(0.3, sigma);\n const half = Math.max(1, Math.round(s * 3));\n const size = 2 * half + 1;\n const data = new Float32Array(size);\n const a = 1 / (Math.sqrt(2 * Math.PI) * s);\n const twoSigma2 = 2 * s * s;\n let sum = 0;\n for (let i = -half; i <= half; i++) {\n const v = a * Math.exp(-(i * i) / twoSigma2);\n data[i + half] = v;\n sum += v;\n }\n for (let i = 0; i < size; i++) data[i] /= sum;\n return { data, half };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":382,"value":"md`Here are several algorithms / strategies you can use to assign every template_edge to a subset of scan edges. Which is best depends on constraints you want to enforce (ordering, polarity, one‑to‑one), the expected amount of outliers/missing edges, and how big the problem is.\n\n...\n\n4) Dynamic programming (sequence alignment) — preserves order, allows gaps\n- If template edges must map in order to scan edges (monotone mapping) and you allow insertions/deletions, use DP (Needleman‑Wunsch style):\n cost(i,j) = min( cost(i-1,j-1) + match_cost(x_ti,x_sj),\n cost(i-1,j) + gap_penalty, // skip template\n cost(i,j-1) + gap_penalty ) // skip scan\n- match_cost could be |x_ti - x_sj| (possibly squared, and include polarity mismatch penalty).\n- Complexity O(N*M) (N = #template edges, M = #scan edges).\n- Pros: optimal under the chosen cost, enforces monotonicity, handles missing edges.\n- Cons: quadratic; choose gap_penalty carefully.\n\n...\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":34,"value":"function clamp(v, lo, hi) {\n return v < lo ? lo : v > hi ? hi : v;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":91,"value":"overlay","pinned":false,"hidden":false,"mode":"ojs"},{"id":356,"value":"SVD = {\n return (await require(\"svd-js\")).SVD;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":378,"value":"import { robocoop } from \"@tomlarkworthy/robocoop-2\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":380,"value":"robocoop","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W37","stdlib":"1","titling":"auto"}