diff --git a/src/components/ResultViewer.vue b/src/components/ResultViewer.vue
index a92839e..64f44b7 100644
--- a/src/components/ResultViewer.vue
+++ b/src/components/ResultViewer.vue
@@ -2,7 +2,7 @@
import { ref, computed, onMounted, watch } from "vue"
import { useAppStore } from "@/stores/app"
import { deskewImage, waitForOpenCV } from "@/lib/deskew"
-import type { RectDatum } from "@/types"
+import type { Datum } from "@/types"
import { DEFAULT_SCALE_PX_PER_MM } from "@/types"
import { Button } from "@/components/ui/button"
import { Input } from "@/components/ui/input"
@@ -59,39 +59,67 @@ watch(scaleInput, (v) => {
const MAX_AUTO_SCALE_DIM = 8192
+/** Estimate the image-pixels-per-mm implied by a single datum. Picks the
+ * best datum by type priority (rect > line > ellipse) and then confidence.
+ * Returns null if no datum gives a usable scale. */
+function pickScaleRef(): { srcPxPerMm: number } | null {
+ const best = [...store.datums].sort((a, b) => {
+ const rank = (d: Datum) =>
+ d.type === "rectangle" ? 0 : d.type === "ellipse" ? 1 : 2
+ const r = rank(a) - rank(b)
+ if (r !== 0) return r
+ return b.confidence - a.confidence
+ })[0]
+ if (!best) return null
+ if (best.type === "rectangle") {
+ if (best.widthMm <= 0 || best.heightMm <= 0) return null
+ const c = best.corners
+ const srcW = Math.max(
+ Math.hypot(c[1].x - c[0].x, c[1].y - c[0].y),
+ Math.hypot(c[2].x - c[3].x, c[2].y - c[3].y),
+ )
+ const srcH = Math.max(
+ Math.hypot(c[3].x - c[0].x, c[3].y - c[0].y),
+ Math.hypot(c[2].x - c[1].x, c[2].y - c[1].y),
+ )
+ const sx = srcW / best.widthMm
+ const sy = srcH / best.heightMm
+ return { srcPxPerMm: Math.max(sx, sy) }
+ }
+ if (best.type === "line") {
+ if (best.lengthMm <= 0) return null
+ const L = Math.hypot(
+ best.endpoints[1].x - best.endpoints[0].x,
+ best.endpoints[1].y - best.endpoints[0].y,
+ )
+ return { srcPxPerMm: L / best.lengthMm }
+ }
+ if (best.diameterMm <= 0) return null
+ // Approximate the ellipse's "diameter" as max of the two semi-axis lengths × 2
+ const vA = Math.hypot(
+ best.axisEndA.x - best.center.x,
+ best.axisEndA.y - best.center.y,
+ )
+ const vB = Math.hypot(
+ best.axisEndB.x - best.center.x,
+ best.axisEndB.y - best.center.y,
+ )
+ return { srcPxPerMm: (2 * Math.max(vA, vB)) / best.diameterMm }
+}
+
function computeAutoScale(): number {
const img = store.loadedImage
- const primary = store.datums.find(
- (d): d is RectDatum => d.type === "rectangle",
- )
- if (!img || !primary) return DEFAULT_SCALE_PX_PER_MM
+ const ref = pickScaleRef()
+ if (!img || !ref || ref.srcPxPerMm <= 0) return DEFAULT_SCALE_PX_PER_MM
- // Approximate source-pixel size of the datum
- const c = primary.corners
- const datumSrcW = Math.max(
- Math.hypot(c[1].x - c[0].x, c[1].y - c[0].y),
- Math.hypot(c[2].x - c[3].x, c[2].y - c[3].y),
- )
- const datumSrcH = Math.max(
- Math.hypot(c[3].x - c[0].x, c[3].y - c[0].y),
- Math.hypot(c[2].x - c[1].x, c[2].y - c[1].y),
- )
-
- // Scale that would make the datum the same pixel size as in source
- const sx =
- datumSrcW > 0 ? datumSrcW / primary.widthMm : 0
- const sy =
- datumSrcH > 0 ? datumSrcH / primary.heightMm : 0
- let autoScale = Math.max(sx, sy)
+ let autoScale = ref.srcPxPerMm
// Clamp so the full output doesn't exceed MAX_AUTO_SCALE_DIM
- const estW = img.naturalWidth * autoScale / Math.max(datumSrcW / primary.widthMm, 0.001)
- const estH = img.naturalHeight * autoScale / Math.max(datumSrcH / primary.heightMm, 0.001)
- if (estW > MAX_AUTO_SCALE_DIM || estH > MAX_AUTO_SCALE_DIM) {
- autoScale *= MAX_AUTO_SCALE_DIM / Math.max(estW, estH)
+ const estMax = Math.max(img.naturalWidth, img.naturalHeight)
+ if (estMax > MAX_AUTO_SCALE_DIM) {
+ autoScale *= MAX_AUTO_SCALE_DIM / estMax
}
- // Round to a clean number
return Math.max(1, Math.round(autoScale * 10) / 10)
}
@@ -134,39 +162,18 @@ const progressPercent = computed(() =>
// Estimated output size — accounts for full warped image, not just datum
const MAX_RGBA_MB = 512
const estimatedOutput = computed(() => {
- const primary = store.datums.find(
- (d): d is RectDatum => d.type === "rectangle",
- )
+ const ref = pickScaleRef()
const img = store.loadedImage
- if (!primary || !img || store.scalePxPerMm <= 0) return null
+ if (!ref || !img || store.scalePxPerMm <= 0 || ref.srcPxPerMm <= 0)
+ return null
- // Datum dimensions in output pixels
- const datumOutW = primary.widthMm * store.scalePxPerMm
- const datumOutH = primary.heightMm * store.scalePxPerMm
+ // source-pixels-per-mm implied by the datum vs. requested output px/mm
+ const avgScale = store.scalePxPerMm / ref.srcPxPerMm
- // Datum dimensions in source pixels (approximate from corner spread)
- const c = primary.corners
- const datumSrcW = Math.max(
- Math.hypot(c[1].x - c[0].x, c[1].y - c[0].y),
- Math.hypot(c[2].x - c[3].x, c[2].y - c[3].y),
- )
- const datumSrcH = Math.max(
- Math.hypot(c[3].x - c[0].x, c[3].y - c[0].y),
- Math.hypot(c[2].x - c[1].x, c[2].y - c[1].y),
- )
-
- // Scale factor from source to output (per-axis average)
- const sx =
- datumSrcW > 0 ? datumOutW / datumSrcW : store.scalePxPerMm
- const sy =
- datumSrcH > 0 ? datumOutH / datumSrcH : store.scalePxPerMm
- const avgScale = (sx + sy) / 2
-
- // Estimated full warped output = source image scaled
const w = Math.round(img.naturalWidth * avgScale)
const h = Math.round(img.naturalHeight * avgScale)
const mb = (w * h * 4) / (1024 * 1024)
- return { w, h, mb, datumW: Math.round(datumOutW), datumH: Math.round(datumOutH) }
+ return { w, h, mb }
})
const tooLarge = computed(
@@ -343,9 +350,6 @@ async function download() {
URL.revokeObjectURL(url)
}
-function hasRects(): boolean {
- return store.datums.some((d) => d.type === "rectangle")
-}
@@ -451,7 +455,9 @@ function hasRects(): boolean {
{{
datum.type === "rectangle"
? `${datum.widthMm}\u00D7${datum.heightMm}mm`
- : `${datum.lengthMm}mm`
+ : datum.type === "line"
+ ? `${datum.lengthMm}mm`
+ : `⌀${datum.diameterMm}mm`
}}conf {{ datum.confidence }}/5
- At least one rectangle datum is required for perspective
- correction.
+ Add at least one datum (rectangle, line, or circle) to run
+ the correction.
- The primary rectangle datum
- defines a
- homography
- — a 3×3 projective
- transform mapping the
- quadrilateral in the source
- image to a true rectangle at
- the specified real-world
- dimensions.
+ The highest-confidence rectangle
+ gives a closed-form warm start via
+ cv::getPerspectiveTransform, fixing the output
+ orientation.
- Secondary datums (additional
- rectangles or line segments)
- provide
+ Each datum is turned into
weighted correction
- factorsshape-based point
+ correspondences
- for the X and Y axes. Each
- secondary datum's contribution
- is weighted by its confidence
- score, refining the scale in
- each axis independently.
+ whose target positions are
+ recomputed from the current
+ homography on every outer pass:
+ Procrustes-fit ideal rectangles,
+ midpoint-preserving line rescales,
+ and radially-snapped ellipse
+ samples that force circles to stay
+ circular.
- The final correction is applied
- as a single
+ cv::findHomography
+ refines the homography by
+ Levenberg–Marquardt on those
+ correspondences; confidence drives
+ per-datum replication. The loop
+ stops once the homography stops
+ moving.
+
+
+ A single
cv::warpPerspective
- call via OpenCV WASM, producing
- the output image at the
+ produces the output at the
requested px/mm scale.
diff --git a/src/lib/datums.ts b/src/lib/datums.ts
index 1bf6531..7115c05 100644
--- a/src/lib/datums.ts
+++ b/src/lib/datums.ts
@@ -1,12 +1,27 @@
import { nanoid } from "nanoid"
-import type { LineDatum, Point, RectDatum, RectPreset } from "@/types"
+import type {
+ CirclePreset,
+ EllipseDatum,
+ LineDatum,
+ Point,
+ RectDatum,
+ RectPreset,
+} from "@/types"
export const RECT_PRESETS: RectPreset[] = [
{ label: "A3", widthMm: 297, heightMm: 420 },
{ label: "A4", widthMm: 210, heightMm: 297 },
{ label: "A5", widthMm: 148, heightMm: 210 },
{ label: "A6", widthMm: 105, heightMm: 148 },
- { label: "15\u00D710 cm", widthMm: 150, heightMm: 100 },
+ { label: "15×10 cm", widthMm: 150, heightMm: 100 },
+]
+
+export const CIRCLE_PRESETS: CirclePreset[] = [
+ { label: "€1", diameterMm: 23.25 },
+ { label: "€2", diameterMm: 25.75 },
+ { label: "US 25¢", diameterMm: 24.26 },
+ { label: "UK 1p", diameterMm: 20.3 },
+ { label: "CD", diameterMm: 120 },
]
const DATUM_COLORS = [
@@ -73,3 +88,21 @@ export function createLineDatum(center: Point, index: number): LineDatum {
label: `Line ${String(index)}`,
}
}
+
+export function createEllipseDatum(
+ center: Point,
+ index: number,
+ preset?: CirclePreset,
+): EllipseDatum {
+ const spread = 80
+ return {
+ id: nanoid(),
+ type: "ellipse",
+ center: { x: center.x, y: center.y },
+ axisEndA: { x: center.x + spread, y: center.y },
+ axisEndB: { x: center.x, y: center.y + spread },
+ diameterMm: preset?.diameterMm ?? 0,
+ confidence: 3,
+ label: preset?.label ?? `Circle ${String(index)}`,
+ }
+}
diff --git a/src/lib/deskew.ts b/src/lib/deskew.ts
index 1a09750..81a0609 100644
--- a/src/lib/deskew.ts
+++ b/src/lib/deskew.ts
@@ -1,126 +1,51 @@
/**
- * deskew.ts — Browser-based perspective correction using OpenCV.js (WASM)
+ * deskew.ts — perspective correction pipeline.
*
- * Accepts N datums (rectangles and/or lines), each with known real-world
- * dimensions and a confidence score (1–5). Minimum: one rectangle.
- *
- * Algorithm:
- * 1. Pick the highest-confidence rectangle as primary reference.
- * 2. getPerspectiveTransform from its 4 corners → initial correction.
- * 3. Project all other datums through that transform and measure them.
- * 4. Compute per-axis weighted scale corrections from all secondary datums.
- * 5. Fold corrections into the destination rectangle, recompute
- * getPerspectiveTransform → single clean perspective matrix.
- * 6. warpPerspective the image.
+ * Delegates the homography solve to src/lib/solver.ts (alternating
+ * minimization around cv.findHomography, which internally runs
+ * Levenberg-Marquardt). This file handles I/O: image loading,
+ * output-bounds computation, warp, PNG export, progress reporting.
*/
import cv from "@techstark/opencv-js"
import type {
- AxisCorrection,
- Datum,
- DatumReport,
DeskewDiagnostics,
DeskewInput,
DeskewResult,
Point,
- RectDatum,
} from "@/types"
+import { solveHomographyForDatums, type Mat3 } from "@/lib/solver"
// Max output dimension in pixels to avoid WASM OOM
// 12288 = ~576MB RGBA at square, but actual images are rarely square
const MAX_OUTPUT_DIM = 12288
-// ─── OpenCV helpers ──────────────────────────────────────────────────────────
+// ─── Small helpers ──────────────────────────────────────────────────────────
-function pointsToMat(points: Point[]): InstanceType {
- const flat = points.flatMap((p) => [p.x, p.y])
- return cv.matFromArray(points.length, 1, cv.CV_32FC2, flat)
-}
-
-function transformPoints(
- points: Point[],
- M: InstanceType,
-): Point[] {
- const src = pointsToMat(points)
- const dst = new cv.Mat()
- cv.perspectiveTransform(src, dst, M)
- const result: Point[] = []
- const data = dst.data32F
- for (let i = 0; i < points.length; i++) {
- const x = data[i * 2]
- const y = data[i * 2 + 1]
- if (x === undefined || y === undefined) continue
- result.push({ x, y })
- }
- src.delete()
- dst.delete()
- return result
-}
-
-function dist(a: Point, b: Point): number {
- return Math.hypot(b.x - a.x, b.y - a.y)
-}
-
-function readMat3x3(M: InstanceType): number[] {
- const d: number[] = []
- for (let r = 0; r < 3; r++) {
- for (let c = 0; c < 3; c++) {
- d.push(M.doubleAt(r, c))
+function projectPoints(h: Mat3, pts: Point[]): Point[] {
+ return pts.map((p) => {
+ const w = h[6] * p.x + h[7] * p.y + h[8]
+ return {
+ x: (h[0] * p.x + h[1] * p.y + h[2]) / w,
+ y: (h[3] * p.x + h[4] * p.y + h[5]) / w,
}
- }
- return d
+ })
}
-/** Row-major 3x3 matrix multiply */
-function mul3x3(A: number[], B: number[]): number[] {
- const R = Array(9).fill(0)
+function mul3x3(A: Mat3, B: Mat3): Mat3 {
+ const R: number[] = Array(9).fill(0)
for (let r = 0; r < 3; r++) {
for (let c = 0; c < 3; c++) {
let sum = 0
for (let k = 0; k < 3; k++) {
- sum +=
- (A[r * 3 + k] ?? 0) * (B[k * 3 + c] ?? 0)
+ sum += (A[r * 3 + k] ?? 0) * (B[k * 3 + c] ?? 0)
}
R[r * 3 + c] = sum
}
}
- return R
+ return R as Mat3
}
-// ─── Validation ──────────────────────────────────────────────────────────────
-
-function pickPrimary(datums: Datum[]): RectDatum {
- const rects = datums.filter(
- (d): d is RectDatum => d.type === "rectangle",
- )
- if (rects.length === 0) {
- throw new Error(
- "At least one rectangle datum is required for perspective correction.",
- )
- }
- rects.sort((a, b) => {
- if (b.confidence !== a.confidence)
- return b.confidence - a.confidence
- const area = (r: RectDatum) =>
- dist(r.corners[0], r.corners[1]) *
- dist(r.corners[0], r.corners[3])
- return area(b) - area(a)
- })
- return rects[0] as RectDatum
-}
-
-/**
- * Convert our app corner order (TL, TR, BR, BL) to the algorithm's
- * expected order (TL, TR, BL, BR) for getPerspectiveTransform.
- */
-function cornersToAlgoOrder(
- corners: [Point, Point, Point, Point],
-): [Point, Point, Point, Point] {
- return [corners[0], corners[1], corners[3], corners[2]]
-}
-
-// ─── Canvas → Blob helper ───────────────────────────────────────────────────
-
function canvasToBlob(
canvas: HTMLCanvasElement,
type = "image/png",
@@ -138,33 +63,32 @@ function canvasToBlob(
})
}
-// ─── Core ────────────────────────────────────────────────────────────────────
-
const log = (tag: string, ...args: unknown[]) => {
console.log(`[deskew:${tag}]`, ...args)
}
+// ─── Core ───────────────────────────────────────────────────────────────────
+
export async function deskewImage(
input: DeskewInput,
): Promise {
const { image, datums, scalePxPerMm: scale, onProgress } = input
- log("start", `${String(datums.length)} datums, scale=${String(scale)} px/mm`)
+ log(
+ "start",
+ `${String(datums.length)} datums, scale=${String(scale)} px/mm`,
+ )
- const TOTAL_STEPS = 7
+ const TOTAL_STEPS = 5
const progress = async (step: number, label: string) => {
- log(`progress`, `[${String(step + 1)}/${String(TOTAL_STEPS)}] ${label}`)
+ log("progress", `[${String(step + 1)}/${String(TOTAL_STEPS)}] ${label}`)
onProgress?.(step, TOTAL_STEPS, label)
- // Yield to let the browser repaint
await new Promise((r) => {
requestAnimationFrame(r)
})
}
if (datums.length === 0) throw new Error("No datums provided.")
- const primary = pickPrimary(datums)
- log("primary", primary.label, `${String(primary.widthMm)}×${String(primary.heightMm)}mm`, `conf=${String(primary.confidence)}`)
-
- // Load source image into OpenCV
+ // Load source image into a canvas
let srcCanvas: HTMLCanvasElement
if (image instanceof HTMLCanvasElement) {
srcCanvas = image
@@ -173,7 +97,10 @@ export async function deskewImage(
srcCanvas = document.createElement("canvas")
srcCanvas.width = image.naturalWidth
srcCanvas.height = image.naturalHeight
- log("input", `img ${String(image.naturalWidth)}×${String(image.naturalHeight)}, drawing to canvas`)
+ log(
+ "input",
+ `img ${String(image.naturalWidth)}×${String(image.naturalHeight)}`,
+ )
const ctx = srcCanvas.getContext("2d")
if (!ctx) throw new Error("Failed to get 2d context")
ctx.drawImage(image, 0, 0)
@@ -181,7 +108,6 @@ export async function deskewImage(
await progress(0, "Loading image into OpenCV")
- // All OpenCV mats to clean up
const mats: InstanceType[] = []
const track = >(m: T): T => {
mats.push(m)
@@ -189,184 +115,41 @@ export async function deskewImage(
}
try {
- log("cv.imread", "reading source canvas into cv.Mat")
const src = track(cv.imread(srcCanvas))
const imgW = src.cols
const imgH = src.rows
- log("cv.imread", `done: ${String(imgW)}×${String(imgH)}, type=${String(src.type())}, channels=${String(src.channels())}`)
-
- // ============================================================
- // STEP 1 — Initial perspective correction from primary rect
- // ============================================================
- await progress(1, "Computing initial homography")
- const pw = primary.widthMm * scale
- const ph = primary.heightMm * scale
- log("step1", `dest rect: ${pw.toFixed(1)}×${ph.toFixed(1)} px`)
-
- const algoCorners = cornersToAlgoOrder(primary.corners)
- log("step1", `corners (algo order): ${JSON.stringify(algoCorners)}`)
- const srcPts = track(pointsToMat(algoCorners))
-
- const dstInit = track(
- pointsToMat([
- { x: 0, y: 0 },
- { x: pw, y: 0 },
- { x: 0, y: ph },
- { x: pw, y: ph },
- ]),
+ log(
+ "cv.imread",
+ `${String(imgW)}×${String(imgH)}, channels=${String(src.channels())}`,
)
- log("step1", "calling getPerspectiveTransform (initial)")
- const mInit = track(
- cv.getPerspectiveTransform(srcPts, dstInit),
+
+ // ========================================================
+ // STEP 1 — Solve homography (outer loop around findHomography)
+ // ========================================================
+ await progress(1, "Solving homography")
+ const solved = solveHomographyForDatums(datums, scale)
+ const H = solved.H
+ log(
+ "solve",
+ `primary=${solved.primaryLabel} (${solved.primaryType}), iters=${String(solved.iterations)}, rms=${solved.rmsPercent.toFixed(3)}%`,
)
- log("step1", `mInit type=${String(mInit.type())}, rows=${String(mInit.rows)}, cols=${String(mInit.cols)}`)
- // ============================================================
- // STEP 2 — Measure secondary datums, accumulate corrections
- // ============================================================
- await progress(2, "Measuring secondary datums")
- let xWSum = 0,
- xWTotal = 0
- let yWSum = 0,
- yWTotal = 0
- const reports: DatumReport[] = []
-
- for (const datum of datums) {
- const w = datum.confidence
-
- if (datum === primary) {
- reports.push({
- label: datum.label,
- type: "rectangle",
- measuredMm: datum.widthMm,
- expectedMm: datum.widthMm,
- errorPercent: 0,
- axisContribution: "both",
- })
- continue
- }
-
- if (datum.type === "line") {
- const pts = transformPoints(
- datum.endpoints as Point[],
- mInit,
- )
- const s = pts[0]
- const e = pts[1]
- if (!s || !e) continue
- const dx = Math.abs(e.x - s.x)
- const dy = Math.abs(e.y - s.y)
- const measured = dist(s, e)
- const expected = datum.lengthMm * scale
- const ratio = expected / measured
-
- const total = dx + dy
- if (total > 1e-6) {
- const xFrac = dx / total
- const yFrac = dy / total
- xWSum += ratio * w * xFrac
- xWTotal += w * xFrac
- yWSum += ratio * w * yFrac
- yWTotal += w * yFrac
- }
-
- reports.push({
- label: datum.label,
- type: "line",
- measuredMm: measured / scale,
- expectedMm: datum.lengthMm,
- errorPercent: Math.abs(1 - ratio) * 100,
- axisContribution: dx > dy ? "x" : "y",
- })
- } else {
- const ac = cornersToAlgoOrder(datum.corners)
- const pts = transformPoints(
- [ac[0], ac[1], ac[2]],
- mInit,
- )
- const tl = pts[0]
- const tr = pts[1]
- const bl = pts[2]
- if (!tl || !tr || !bl) continue
- const mW = dist(tl, tr)
- const mH = dist(tl, bl)
- const xR = (datum.widthMm * scale) / mW
- const yR = (datum.heightMm * scale) / mH
-
- xWSum += xR * w
- xWTotal += w
- yWSum += yR * w
- yWTotal += w
-
- reports.push({
- label: datum.label,
- type: "rectangle",
- measuredMm: mW / scale,
- expectedMm: datum.widthMm,
- errorPercent:
- (Math.abs(1 - xR) + Math.abs(1 - yR)) * 50,
- axisContribution: "both",
- })
- }
- }
-
- // ============================================================
- // STEP 3 — Weighted corrections (1.0 = no secondary data)
- // ============================================================
- await progress(3, "Computing axis corrections")
- const xCorr: AxisCorrection = {
- ratio: xWTotal > 0 ? xWSum / xWTotal : 1.0,
- totalWeight: xWTotal,
- }
- const yCorr: AxisCorrection = {
- ratio: yWTotal > 0 ? yWSum / yWTotal : 1.0,
- totalWeight: yWTotal,
- }
- log("step3", `xCorr=${xCorr.ratio.toFixed(4)} (w=${xCorr.totalWeight.toFixed(1)}), yCorr=${yCorr.ratio.toFixed(4)} (w=${yCorr.totalWeight.toFixed(1)})`)
-
- // ============================================================
- // STEP 4 — Fold corrections, recompute transform
- // ============================================================
- await progress(4, "Recomputing final transform")
- const pwFinal = pw * xCorr.ratio
- const phFinal = ph * yCorr.ratio
- log("step4", `final dest rect: ${pwFinal.toFixed(1)}×${phFinal.toFixed(1)} px`)
-
- const dstFinal = track(
- pointsToMat([
- { x: 0, y: 0 },
- { x: pwFinal, y: 0 },
- { x: 0, y: phFinal },
- { x: pwFinal, y: phFinal },
- ]),
- )
- log("step4", "calling getPerspectiveTransform (final)")
- const mFinal = track(
- cv.getPerspectiveTransform(srcPts, dstFinal),
- )
- log("step4", `mFinal type=${String(mFinal.type())}, rows=${String(mFinal.rows)}, cols=${String(mFinal.cols)}`)
-
- // ============================================================
- // STEP 5 — Output bounds + translation shift
- // ============================================================
- await progress(5, "Computing output bounds")
+ // ========================================================
+ // STEP 2 — Compute output bounds and translation shift
+ // ========================================================
+ await progress(2, "Computing output bounds")
const imgCorners: Point[] = [
{ x: 0, y: 0 },
{ x: imgW, y: 0 },
{ x: 0, y: imgH },
{ x: imgW, y: imgH },
]
- const warped = transformPoints(imgCorners, mFinal)
- if (warped.length < 4) {
- throw new Error(
- "Perspective transform produced invalid bounds",
- )
- }
+ const warped = projectPoints(H, imgCorners)
- let xMin = Infinity,
- yMin = Infinity,
- xMax = -Infinity,
- yMax = -Infinity
+ let xMin = Infinity
+ let yMin = Infinity
+ let xMax = -Infinity
+ let yMax = -Infinity
for (const c of warped) {
xMin = Math.min(xMin, c.x)
yMin = Math.min(yMin, c.y)
@@ -376,10 +159,10 @@ export async function deskewImage(
let outW = Math.ceil(xMax - xMin)
let outH = Math.ceil(yMax - yMin)
- log("step5", `bounds: x=[${xMin.toFixed(1)}, ${xMax.toFixed(1)}], y=[${yMin.toFixed(1)}, ${yMax.toFixed(1)}]`)
- log("step5", `raw output: ${String(outW)}×${String(outH)} px`)
-
- // Guard against absurd output sizes that crash WASM
+ log(
+ "bounds",
+ `x=[${xMin.toFixed(1)},${xMax.toFixed(1)}], y=[${yMin.toFixed(1)},${yMax.toFixed(1)}]`,
+ )
if (outW <= 0 || outH <= 0) {
throw new Error(
`Invalid output dimensions: ${String(outW)}×${String(outH)}`,
@@ -388,30 +171,31 @@ export async function deskewImage(
let downscale = 1
if (outW > MAX_OUTPUT_DIM || outH > MAX_OUTPUT_DIM) {
downscale = MAX_OUTPUT_DIM / Math.max(outW, outH)
- log("step5", `CLAMPING from ${String(outW)}×${String(outH)} by factor ${downscale.toFixed(4)}`)
outW = Math.ceil(outW * downscale)
outH = Math.ceil(outH * downscale)
+ log("bounds", `clamped by ${downscale.toFixed(4)} → ${String(outW)}×${String(outH)}`)
}
- log("step5", `final output: ${String(outW)}×${String(outH)} px (${String(Math.round(outW * outH * 4 / 1024 / 1024))} MB RGBA)`)
- const mData: number[] = readMat3x3(mFinal)
- // Translate so the top-left warped corner is at (0,0),
- // then scale down if we clamped the output size.
- const tShift: number[] = [
- downscale, 0, -xMin * downscale,
- 0, downscale, -yMin * downscale,
- 0, 0, 1,
+ // Compose a shift (translation + optional downscale) with H so the
+ // top-left corner of the warped image lands at (0, 0).
+ const tShift: Mat3 = [
+ downscale,
+ 0,
+ -xMin * downscale,
+ 0,
+ downscale,
+ -yMin * downscale,
+ 0,
+ 0,
+ 1,
]
- const mOutData: number[] = mul3x3(tShift, mData)
- const mOut = track(
- cv.matFromArray(3, 3, cv.CV_64FC1, mOutData),
- )
+ const mOutData = mul3x3(tShift, H)
+ const mOut = track(cv.matFromArray(3, 3, cv.CV_64FC1, mOutData))
- // ============================================================
- // STEP 6 — Warp
- // ============================================================
- await progress(6, "Warping image (this may take a moment)")
- log("step6", "calling warpPerspective...")
+ // ========================================================
+ // STEP 3 — Warp
+ // ========================================================
+ await progress(3, "Warping image")
const dstMat = track(new cv.Mat())
cv.warpPerspective(
src,
@@ -423,36 +207,33 @@ export async function deskewImage(
new cv.Scalar(0, 0, 0, 0),
)
- log("step6", `warpPerspective done, dstMat: ${String(dstMat.cols)}×${String(dstMat.rows)}, type=${String(dstMat.type())}`)
-
- log("export", "cv.imshow to canvas")
+ // ========================================================
+ // STEP 4 — Export
+ // ========================================================
+ await progress(4, "Encoding output")
const outCanvas = document.createElement("canvas")
outCanvas.width = outW
outCanvas.height = outH
cv.imshow(outCanvas, dstMat)
-
- log("export", "canvas.toBlob (PNG)")
const blob = await canvasToBlob(outCanvas, "image/png", 0.95)
- log("export", `blob size: ${String(Math.round(blob.size / 1024))} KB`)
+ log("export", `blob ${String(Math.round(blob.size / 1024))} KB`)
const diagnostics: DeskewDiagnostics = {
- primaryDatum: primary.label,
- xCorrection: xCorr,
- yCorrection: yCorr,
- perDatum: reports,
+ primaryDatum: solved.primaryLabel,
+ iterations: solved.iterations,
+ finalRMSPercent: solved.rmsPercent,
+ perDatum: solved.reports,
outputWidthPx: outW,
outputHeightPx: outH,
}
- log("done", "success")
return { correctedImageBlob: blob, diagnostics }
} finally {
- // Always clean up all OpenCV mats, even on error
for (const m of mats) {
try {
m.delete()
} catch {
- // already deleted or invalid — ignore
+ // already deleted — ignore
}
}
}
@@ -464,31 +245,22 @@ let cvReady = false
/** Wait for OpenCV WASM to initialize. Call once at app startup. */
export function waitForOpenCV(): Promise {
- log("opencv", "waitForOpenCV called, cvReady=" + String(cvReady))
return new Promise((resolve) => {
if (cvReady) {
- log("opencv", "already ready")
resolve()
return
}
-
- // Test if WASM is actually functional by trying to create a mat
try {
- log("opencv", "probing cv.Mat()...")
const test = new cv.Mat()
test.delete()
cvReady = true
- log("opencv", "probe succeeded, WASM ready")
resolve()
return
} catch {
- log("opencv", "probe failed, waiting for onRuntimeInitialized")
- // Not ready yet, wait for callback
+ // Runtime not ready yet
}
-
cv.onRuntimeInitialized = () => {
cvReady = true
- log("opencv", "onRuntimeInitialized fired, WASM ready")
resolve()
}
})
diff --git a/src/lib/solver.ts b/src/lib/solver.ts
new file mode 100644
index 0000000..0cce662
--- /dev/null
+++ b/src/lib/solver.ts
@@ -0,0 +1,843 @@
+/**
+ * solver.ts — alternating-minimization solver around cv.findHomography.
+ *
+ * Each datum is converted to point correspondences (src → dst) whose dst
+ * positions are recomputed from the current H and the datum's shape
+ * constraint on every outer iteration. findHomography then refines H (DLT
+ * initial estimate + internal Levenberg-Marquardt). We loop until H stops
+ * changing.
+ *
+ * Shape constraints per datum type:
+ * rectangle (primary) — hard-anchored axis-aligned destination
+ * rectangle (other) — Procrustes-fit ideal (w × h) rect to projected corners
+ * line — preserve projected midpoint + direction, rescale to L
+ * ellipse — sample N points; radially snap projections to a circle
+ * of known diameter centred on the projected
+ * user-marked centre point
+ *
+ * Confidence is the per-datum weight; we realise it as correspondence
+ * replication (findHomography has no native weighting).
+ */
+
+import cv from "@techstark/opencv-js"
+import type {
+ Datum,
+ DatumReport,
+ DatumType,
+ EllipseDatum,
+ LineDatum,
+ Point,
+ RectDatum,
+} from "@/types"
+
+// ─── Tunables ───────────────────────────────────────────────────────────────
+
+/** Extra weight on primary-datum anchor correspondences, on top of its
+ * confidence. Keeps the output gauge (orientation/position) stable across
+ * iterations while still letting consistent secondary data nudge H. */
+const PRIMARY_GAUGE_BOOST = 3
+
+/** Points sampled per ellipse datum (more points = tighter fit but more
+ * replicated correspondences). 12 gives good angular coverage. */
+const ELLIPSE_SAMPLES = 12
+
+const MAX_OUTER_ITERS = 30
+
+/** Convergence threshold: max entrywise change in H between successive
+ * iterations, relative to the largest entry of H. A relative metric
+ * treats the small perspective entries (h[6], h[7] ≈ 1e-4) and the O(1)
+ * affine entries on the same footing. */
+const CONVERGENCE_TOL = 1e-6
+
+// ─── 3×3 matrix helpers ─────────────────────────────────────────────────────
+
+export type Mat3 = [number, number, number, number, number, number, number, number, number]
+
+function readMat3x3(M: InstanceType): Mat3 {
+ const d: number[] = []
+ for (let r = 0; r < 3; r++) {
+ for (let c = 0; c < 3; c++) {
+ d.push(M.doubleAt(r, c))
+ }
+ }
+ return d as Mat3
+}
+
+function projectPoint(h: Mat3, p: Point): Point {
+ const w = h[6] * p.x + h[7] * p.y + h[8]
+ return {
+ x: (h[0] * p.x + h[1] * p.y + h[2]) / w,
+ y: (h[3] * p.x + h[4] * p.y + h[5]) / w,
+ }
+}
+
+function normalized(h: Mat3): Mat3 {
+ const s = h[8] !== 0 ? h[8] : 1
+ return h.map((v) => v / s) as Mat3
+}
+
+/** Relative max-entry change between two homographies, normalised by the
+ * larger-magnitude entry in either matrix. Returns a dimensionless fraction
+ * so a single convergence threshold is meaningful across affine and
+ * perspective parameters. */
+function relativeMaxDiff(a: Mat3, b: Mat3): number {
+ let diff = 0
+ let scale = 0
+ for (let i = 0; i < 9; i++) {
+ const av = Math.abs(a[i] ?? 0)
+ const bv = Math.abs(b[i] ?? 0)
+ if (av > scale) scale = av
+ if (bv > scale) scale = bv
+ const d = Math.abs((a[i] ?? 0) - (b[i] ?? 0))
+ if (d > diff) diff = d
+ }
+ return scale > 0 ? diff / scale : diff
+}
+
+// ─── Geometric helpers ──────────────────────────────────────────────────────
+
+function dist(a: Point, b: Point): number {
+ return Math.hypot(b.x - a.x, b.y - a.y)
+}
+
+function centroid(pts: Point[]): Point {
+ let sx = 0
+ let sy = 0
+ for (const p of pts) {
+ sx += p.x
+ sy += p.y
+ }
+ return { x: sx / pts.length, y: sy / pts.length }
+}
+
+/**
+ * Best rigid transform (rotation + translation, no scale) aligning src to dst.
+ * Closed-form 2D Procrustes. Returns (R, t) with R @ src_i + t ≈ dst_i.
+ */
+function procrustes2D(
+ src: Point[],
+ dst: Point[],
+): { cos: number; sin: number; tx: number; ty: number } {
+ const cs = centroid(src)
+ const cd = centroid(dst)
+ let hxx = 0
+ let hxy = 0
+ let hyx = 0
+ let hyy = 0
+ for (let i = 0; i < src.length; i++) {
+ const s = src[i]
+ const d = dst[i]
+ if (!s || !d) continue
+ const sx = s.x - cs.x
+ const sy = s.y - cs.y
+ const dx = d.x - cd.x
+ const dy = d.y - cd.y
+ hxx += sx * dx
+ hxy += sx * dy
+ hyx += sy * dx
+ hyy += sy * dy
+ }
+ // argmax over θ of tr(R · H^T) = cos·(hxx+hyy) + sin·(hyx−hxy)
+ const theta = Math.atan2(hyx - hxy, hxx + hyy)
+ const cos = Math.cos(theta)
+ const sin = Math.sin(theta)
+ const tx = cd.x - (cos * cs.x - sin * cs.y)
+ const ty = cd.y - (sin * cs.x + cos * cs.y)
+ return { cos, sin, tx, ty }
+}
+
+function applyRT(
+ p: Point,
+ R: { cos: number; sin: number; tx: number; ty: number },
+): Point {
+ return {
+ x: R.cos * p.x - R.sin * p.y + R.tx,
+ y: R.sin * p.x + R.cos * p.y + R.ty,
+ }
+}
+
+/**
+ * Compute the 3×3 symmetric conic matrix E for an ellipse parameterised by
+ * center and two (not necessarily orthogonal) conjugate semi-axes vA, vB.
+ * Parametrically: p(t) = center + vA cos t + vB sin t. The ellipse's
+ * quadratic form is (p − c)^T Q (p − c) = 1 with Q = (M M^T)^{-1} where
+ * M = [vA | vB] (2×2). The homogeneous conic matrix is then assembled so
+ * that [x y 1] E [x y 1]^T = 0 on the ellipse.
+ */
+function ellipseMatrix(
+ center: Point,
+ axisEndA: Point,
+ axisEndB: Point,
+): number[][] {
+ const ax = axisEndA.x - center.x
+ const ay = axisEndA.y - center.y
+ const bx = axisEndB.x - center.x
+ const by = axisEndB.y - center.y
+ // M M^T = [[ax²+bx², ax·ay+bx·by], [·, ay²+by²]]
+ const m00 = ax * ax + bx * bx
+ const m01 = ax * ay + bx * by
+ const m11 = ay * ay + by * by
+ const det = m00 * m11 - m01 * m01
+ if (Math.abs(det) < 1e-12) {
+ throw new Error("Ellipse is degenerate (axes are collinear).")
+ }
+ // Q = (MM^T)^{-1}
+ const q00 = m11 / det
+ const q01 = -m01 / det
+ const q11 = m00 / det
+ const cx = center.x
+ const cy = center.y
+ // E (homogeneous): p^T E p = 0 with p = (x, y, 1)
+ const e02 = -(q00 * cx + q01 * cy)
+ const e12 = -(q01 * cx + q11 * cy)
+ const e22 = q00 * cx * cx + 2 * q01 * cx * cy + q11 * cy * cy - 1
+ return [
+ [q00, q01, e02],
+ [q01, q11, e12],
+ [e02, e12, e22],
+ ]
+}
+
+/** Sample N image-space points along the user-drawn ellipse. */
+function sampleEllipse(
+ center: Point,
+ axisEndA: Point,
+ axisEndB: Point,
+ n: number,
+): Point[] {
+ const vAx = axisEndA.x - center.x
+ const vAy = axisEndA.y - center.y
+ const vBx = axisEndB.x - center.x
+ const vBy = axisEndB.y - center.y
+ const out: Point[] = []
+ for (let i = 0; i < n; i++) {
+ const t = (2 * Math.PI * i) / n
+ const c = Math.cos(t)
+ const s = Math.sin(t)
+ out.push({
+ x: center.x + vAx * c + vBx * s,
+ y: center.y + vAy * c + vBy * s,
+ })
+ }
+ return out
+}
+
+// ─── Primary selection (gauge) ──────────────────────────────────────────────
+
+/** The "primary" datum fixes the output gauge (rotation, translation, scale)
+ * and provides the 4-point warm-start for the homography. We prefer
+ * rectangles (strongest gauge: 4 corners, aspect ratio), then ellipses
+ * (4 conjugate-axis samples fully pin H), then lines (weakest: only 2 real
+ * correspondences + synthetic perpendiculars). Within a type class, higher
+ * confidence and larger image size break ties. */
+type Primary =
+ | { kind: "rect"; datum: RectDatum }
+ | { kind: "line"; datum: LineDatum }
+ | { kind: "ellipse"; datum: EllipseDatum }
+
+function pickPrimary(datums: Datum[]): Primary {
+ if (datums.length === 0) throw new Error("No datums provided.")
+ const typeRank = (d: Datum): number =>
+ d.type === "rectangle" ? 0 : d.type === "ellipse" ? 1 : 2
+ const sizeKey = (d: Datum): number => {
+ if (d.type === "rectangle")
+ return (
+ dist(d.corners[0], d.corners[1]) *
+ dist(d.corners[0], d.corners[3])
+ )
+ if (d.type === "line") return dist(d.endpoints[0], d.endpoints[1])
+ return dist(d.center, d.axisEndA) * dist(d.center, d.axisEndB)
+ }
+ const sorted = [...datums].sort((a, b) => {
+ const tr = typeRank(a) - typeRank(b)
+ if (tr !== 0) return tr
+ if (b.confidence !== a.confidence)
+ return b.confidence - a.confidence
+ return sizeKey(b) - sizeKey(a)
+ })
+ const best = sorted[0] as Datum
+ if (best.type === "rectangle") return { kind: "rect", datum: best }
+ if (best.type === "line") return { kind: "line", datum: best }
+ return { kind: "ellipse", datum: best }
+}
+
+function primaryLabel(primary: Primary): string {
+ return primary.datum.label
+}
+
+// ─── Correspondence builders per datum type ─────────────────────────────────
+
+interface Correspondence {
+ src: Point
+ dst: Point
+ weight: number // integer replication count
+}
+
+function primaryRectCorrespondences(
+ rect: RectDatum,
+ scale: number,
+): Correspondence[] {
+ const w = rect.widthMm * scale
+ const h = rect.heightMm * scale
+ // Corner order TL, TR, BR, BL matches the RectDatum contract
+ const targets: [Point, Point, Point, Point] = [
+ { x: 0, y: 0 },
+ { x: w, y: 0 },
+ { x: w, y: h },
+ { x: 0, y: h },
+ ]
+ const weight = Math.max(
+ 1,
+ Math.round(rect.confidence * PRIMARY_GAUGE_BOOST),
+ )
+ return rect.corners.map((src, i) => ({
+ src,
+ dst: targets[i] as Point,
+ weight,
+ }))
+}
+
+/** Line primary: 2 real endpoints + 2 synthetic perpendicular points at
+ * the same image-space distance. Fixes the along-line scale exactly and
+ * assumes isotropic image scale perpendicular to the line (good warm-start;
+ * LM refines if other datums contradict it). */
+function primaryLineCorrespondences(
+ line: LineDatum,
+ scale: number,
+): Correspondence[] {
+ const L = line.lengthMm * scale
+ const p0 = line.endpoints[0]
+ const p1 = line.endpoints[1]
+ const dx = p1.x - p0.x
+ const dy = p1.y - p0.y
+ // 90° left-rotated copy of the line, same image length
+ const p2: Point = { x: p0.x - dy, y: p0.y + dx }
+ const p3: Point = { x: p1.x - dy, y: p1.y + dx }
+ const srcPts: [Point, Point, Point, Point] = [p0, p1, p2, p3]
+ const targets: [Point, Point, Point, Point] = [
+ { x: 0, y: 0 },
+ { x: L, y: 0 },
+ { x: 0, y: L },
+ { x: L, y: L },
+ ]
+ const weight = Math.max(
+ 1,
+ Math.round(line.confidence * PRIMARY_GAUGE_BOOST),
+ )
+ return srcPts.map((src, i) => ({
+ src,
+ dst: targets[i] as Point,
+ weight,
+ }))
+}
+
+/** Ellipse primary: 4 points on the user-drawn ellipse at t = 0, π/2, π,
+ * 3π/2 (= center ± vA, center ± vB). Anchored to a world-circle of the
+ * known diameter. Fixes translation + rotation + scale via the axisEndA
+ * direction convention. */
+function primaryEllipseCorrespondences(
+ ellipse: EllipseDatum,
+ scale: number,
+): Correspondence[] {
+ const r = (ellipse.diameterMm * scale) / 2
+ const c = ellipse.center
+ const vAx = ellipse.axisEndA.x - c.x
+ const vAy = ellipse.axisEndA.y - c.y
+ const vBx = ellipse.axisEndB.x - c.x
+ const vBy = ellipse.axisEndB.y - c.y
+ // Axes must not be (near-)collinear — otherwise the 4 anchor points are
+ // collinear and getPerspectiveTransform gives a garbage homography.
+ const cross = vAx * vBy - vAy * vBx
+ if (Math.abs(cross) < 1e-6) {
+ throw new Error(
+ `Ellipse "${ellipse.label}" has collinear axes; drag its handles apart before solving.`,
+ )
+ }
+ const srcPts: [Point, Point, Point, Point] = [
+ { x: c.x + vAx, y: c.y + vAy }, // t = 0 → (+r, 0)
+ { x: c.x + vBx, y: c.y + vBy }, // t = π/2 → (0, +r)
+ { x: c.x - vAx, y: c.y - vAy }, // t = π → (−r, 0)
+ { x: c.x - vBx, y: c.y - vBy }, // t = 3π/2 → (0, −r)
+ ]
+ const targets: [Point, Point, Point, Point] = [
+ { x: r, y: 0 },
+ { x: 0, y: r },
+ { x: -r, y: 0 },
+ { x: 0, y: -r },
+ ]
+ const weight = Math.max(
+ 1,
+ Math.round(ellipse.confidence * PRIMARY_GAUGE_BOOST),
+ )
+ return srcPts.map((src, i) => ({
+ src,
+ dst: targets[i] as Point,
+ weight,
+ }))
+}
+
+function primaryAnchors(primary: Primary, scale: number): Correspondence[] {
+ if (primary.kind === "rect")
+ return primaryRectCorrespondences(primary.datum, scale)
+ if (primary.kind === "line")
+ return primaryLineCorrespondences(primary.datum, scale)
+ return primaryEllipseCorrespondences(primary.datum, scale)
+}
+
+function secondaryRectCorrespondences(
+ rect: RectDatum,
+ H: Mat3,
+ scale: number,
+): Correspondence[] {
+ const w = rect.widthMm * scale
+ const h = rect.heightMm * scale
+ // Ideal (w × h) rect in local frame, same corner order
+ const ideal: Point[] = [
+ { x: 0, y: 0 },
+ { x: w, y: 0 },
+ { x: w, y: h },
+ { x: 0, y: h },
+ ]
+ // Project the user corners through current H to output space
+ const projected = rect.corners.map((c) => projectPoint(H, c))
+ // Best rigid placement of the ideal rect to those projections
+ const rt = procrustes2D(ideal, projected)
+ const targets = ideal.map((p) => applyRT(p, rt))
+ const weight = Math.max(1, rect.confidence)
+ return rect.corners.map((src, i) => ({
+ src,
+ dst: targets[i] as Point,
+ weight,
+ }))
+}
+
+function lineCorrespondences(
+ line: LineDatum,
+ H: Mat3,
+ scale: number,
+): Correspondence[] {
+ const expected = line.lengthMm * scale
+ const p0 = projectPoint(H, line.endpoints[0])
+ const p1 = projectPoint(H, line.endpoints[1])
+ const measured = dist(p0, p1)
+ if (measured < 1e-6) {
+ // Degenerate line, no useful correspondence this iteration
+ return []
+ }
+ const mx = (p0.x + p1.x) / 2
+ const my = (p0.y + p1.y) / 2
+ const ux = (p1.x - p0.x) / measured
+ const uy = (p1.y - p0.y) / measured
+ const halfL = expected / 2
+ const target0: Point = { x: mx - ux * halfL, y: my - uy * halfL }
+ const target1: Point = { x: mx + ux * halfL, y: my + uy * halfL }
+ const weight = Math.max(1, line.confidence)
+ return [
+ { src: line.endpoints[0], dst: target0, weight },
+ { src: line.endpoints[1], dst: target1, weight },
+ ]
+}
+
+function ellipseCorrespondences(
+ ellipse: EllipseDatum,
+ H: Mat3,
+ scale: number,
+): Correspondence[] {
+ const r = (ellipse.diameterMm * scale) / 2
+ const samples = sampleEllipse(
+ ellipse.center,
+ ellipse.axisEndA,
+ ellipse.axisEndB,
+ ELLIPSE_SAMPLES,
+ )
+ const projected = samples.map((p) => projectPoint(H, p))
+ // Snap radially around the projected center of the world circle.
+ // Under general perspective the centroid of boundary-point projections
+ // is NOT the image of the center (it drifts with the foreshortening),
+ // so we use the image of the user-marked center directly.
+ const c = projectPoint(H, ellipse.center)
+ const weight = Math.max(1, ellipse.confidence)
+ const out: Correspondence[] = []
+ for (let i = 0; i < samples.length; i++) {
+ const src = samples[i]
+ const q = projected[i]
+ if (!src || !q) continue
+ const dx = q.x - c.x
+ const dy = q.y - c.y
+ const d = Math.hypot(dx, dy)
+ if (d < 1e-6) continue
+ out.push({
+ src,
+ dst: { x: c.x + (dx / d) * r, y: c.y + (dy / d) * r },
+ weight,
+ })
+ }
+ return out
+}
+
+function buildCorrespondences(
+ datums: Datum[],
+ primary: Primary,
+ H: Mat3,
+ scale: number,
+): Correspondence[] {
+ const all: Correspondence[] = []
+ for (const d of datums) {
+ if (d === primary.datum) {
+ all.push(...primaryAnchors(primary, scale))
+ continue
+ }
+ if (d.type === "rectangle") {
+ all.push(...secondaryRectCorrespondences(d, H, scale))
+ } else if (d.type === "line") {
+ all.push(...lineCorrespondences(d, H, scale))
+ } else {
+ all.push(...ellipseCorrespondences(d, H, scale))
+ }
+ }
+ return all
+}
+
+// ─── findHomography wrapper ─────────────────────────────────────────────────
+
+function solveHomography(
+ correspondences: Correspondence[],
+): Mat3 | null {
+ // Replicate by weight to emulate per-correspondence weighting
+ const src: number[] = []
+ const dst: number[] = []
+ let n = 0
+ for (const c of correspondences) {
+ for (let i = 0; i < c.weight; i++) {
+ src.push(c.src.x, c.src.y)
+ dst.push(c.dst.x, c.dst.y)
+ n++
+ }
+ }
+ if (n < 4) return null
+ const srcMat = cv.matFromArray(n, 1, cv.CV_32FC2, src)
+ const dstMat = cv.matFromArray(n, 1, cv.CV_32FC2, dst)
+ let H: InstanceType | null = null
+ try {
+ H = cv.findHomography(srcMat, dstMat, 0)
+ if (H.rows !== 3 || H.cols !== 3) return null
+ return normalized(readMat3x3(H))
+ } finally {
+ srcMat.delete()
+ dstMat.delete()
+ if (H) H.delete()
+ }
+}
+
+// ─── Post-fit residual reporting ────────────────────────────────────────────
+
+interface RawReport {
+ label: string
+ type: DatumType
+ expectedMm: number
+ measuredMm: number
+ residuals: number[] // dimensionless fractions
+ details: string
+}
+
+function residualForLine(
+ line: LineDatum,
+ H: Mat3,
+ scale: number,
+ isPrimary: boolean,
+): RawReport {
+ const p0 = projectPoint(H, line.endpoints[0])
+ const p1 = projectPoint(H, line.endpoints[1])
+ const measuredMm = dist(p0, p1) / scale
+ const residual =
+ line.lengthMm > 0
+ ? (measuredMm - line.lengthMm) / line.lengthMm
+ : 0
+ const prefix = isPrimary ? "primary · " : ""
+ return {
+ label: line.label,
+ type: "line",
+ expectedMm: line.lengthMm,
+ measuredMm,
+ residuals: [residual],
+ details: `${prefix}length ${(residual * 100).toFixed(2)}%`,
+ }
+}
+
+function residualForRect(
+ rect: RectDatum,
+ H: Mat3,
+ scale: number,
+ isPrimary: boolean,
+): RawReport {
+ const p = rect.corners.map((c) => projectPoint(H, c))
+ const w = rect.widthMm
+ const h = rect.heightMm
+ const sides = [
+ { got: dist(p[0] as Point, p[1] as Point) / scale, exp: w }, // TL-TR
+ { got: dist(p[1] as Point, p[2] as Point) / scale, exp: h }, // TR-BR
+ { got: dist(p[2] as Point, p[3] as Point) / scale, exp: w }, // BR-BL
+ { got: dist(p[3] as Point, p[0] as Point) / scale, exp: h }, // BL-TL
+ ]
+ const edgeRes = sides.map((s) => (s.exp > 0 ? (s.got - s.exp) / s.exp : 0))
+ // Perpendicularity at TL and TR: cosine between adjacent edges
+ const e0x = (p[1] as Point).x - (p[0] as Point).x
+ const e0y = (p[1] as Point).y - (p[0] as Point).y
+ const e1x = (p[3] as Point).x - (p[0] as Point).x
+ const e1y = (p[3] as Point).y - (p[0] as Point).y
+ const e2x = (p[2] as Point).x - (p[1] as Point).x
+ const e2y = (p[2] as Point).y - (p[1] as Point).y
+ const cosTL =
+ (e0x * e1x + e0y * e1y) /
+ (Math.hypot(e0x, e0y) * Math.hypot(e1x, e1y) + 1e-12)
+ const cosTR =
+ (-e0x * e2x + -e0y * e2y) /
+ (Math.hypot(e0x, e0y) * Math.hypot(e2x, e2y) + 1e-12)
+ const residuals = [...edgeRes, cosTL, cosTR]
+ const avgEdge =
+ (Math.abs(edgeRes[0] ?? 0) +
+ Math.abs(edgeRes[1] ?? 0) +
+ Math.abs(edgeRes[2] ?? 0) +
+ Math.abs(edgeRes[3] ?? 0)) /
+ 4
+ // |cos θ| ≈ |90° − θ| in radians for small deviations; convert to degrees
+ const perpDeg =
+ ((Math.asin(Math.min(1, Math.abs(cosTL))) +
+ Math.asin(Math.min(1, Math.abs(cosTR)))) /
+ 2) *
+ (180 / Math.PI)
+ const measuredWidth =
+ ((sides[0]?.got ?? 0) + (sides[2]?.got ?? 0)) / 2
+ const prefix = isPrimary ? "primary · " : ""
+ return {
+ label: rect.label,
+ type: "rectangle",
+ expectedMm: w,
+ measuredMm: measuredWidth,
+ residuals,
+ details: `${prefix}edge ${(avgEdge * 100).toFixed(2)}%, perp Δ${perpDeg.toFixed(2)}°`,
+ }
+}
+
+function residualForEllipse(
+ ellipse: EllipseDatum,
+ H: Mat3,
+ scale: number,
+ isPrimary: boolean,
+): RawReport {
+ // Pull the image ellipse back to world space via C = H^T E H
+ const E = ellipseMatrix(ellipse.center, ellipse.axisEndA, ellipse.axisEndB)
+ const HT = [
+ [H[0], H[3], H[6]],
+ [H[1], H[4], H[7]],
+ [H[2], H[5], H[8]],
+ ]
+ const Hm = [
+ [H[0], H[1], H[2]],
+ [H[3], H[4], H[5]],
+ [H[6], H[7], H[8]],
+ ]
+ // C = HT · E · Hm
+ const EH: number[][] = [
+ [0, 0, 0],
+ [0, 0, 0],
+ [0, 0, 0],
+ ]
+ for (let i = 0; i < 3; i++) {
+ for (let j = 0; j < 3; j++) {
+ let s = 0
+ for (let k = 0; k < 3; k++) {
+ s += (E[i]?.[k] ?? 0) * (Hm[k]?.[j] ?? 0)
+ }
+ ;(EH[i] as number[])[j] = s
+ }
+ }
+ const C: number[][] = [
+ [0, 0, 0],
+ [0, 0, 0],
+ [0, 0, 0],
+ ]
+ for (let i = 0; i < 3; i++) {
+ for (let j = 0; j < 3; j++) {
+ let s = 0
+ for (let k = 0; k < 3; k++) {
+ s += (HT[i]?.[k] ?? 0) * (EH[k]?.[j] ?? 0)
+ }
+ ;(C[i] as number[])[j] = s
+ }
+ }
+ const a = C[0]?.[0] ?? 0
+ const b = C[0]?.[1] ?? 0
+ const c = C[1]?.[1] ?? 0
+ const d = C[0]?.[2] ?? 0
+ const e = C[1]?.[2] ?? 0
+ const f = C[2]?.[2] ?? 0
+ const sum = a + c
+ const iso = sum !== 0 ? (a - c) / sum : 0
+ const skew = sum !== 0 ? (2 * b) / sum : 0
+ const rExp = (ellipse.diameterMm * scale) / 2
+
+ // Geometric-mean radius from the general (possibly non-circular) conic.
+ // After centring, the conic is u^T A u = K with A = [[a, b], [b, c]] and
+ // K = a·x0² + 2b·x0·y0 + c·y0² − f, (x0, y0) = −A^{-1}(d, e).
+ // Semi-axes are √(K/λ_±), so √(semi_major · semi_minor) = √(K / √det A).
+ // This coincides with the circle radius when a = c and b = 0, and stays
+ // meaningful (= area-equivalent radius) while the solver drives the
+ // conic toward being a circle.
+ const det = a * c - b * b
+ let rMeasured = 0
+ if (det > 0) {
+ const x0 = (b * e - c * d) / det
+ const y0 = (b * d - a * e) / det
+ const K = a * x0 * x0 + 2 * b * x0 * y0 + c * y0 * y0 - f
+ if (K > 0) rMeasured = Math.sqrt(K / Math.sqrt(det))
+ }
+ // The dia residual drives the solver; we divide by rExp for scale-free.
+ // When the conic is still clearly non-circular the per-axis semi-axes
+ // differ, so we penalise by the geometric-mean radius — which is what
+ // we want to land at rExp once iso/skew have been driven to zero.
+ const diaRes = rExp > 0 ? (rMeasured - rExp) / rExp : 0
+ const prefix = isPrimary ? "primary · " : ""
+ return {
+ label: ellipse.label,
+ type: "ellipse",
+ expectedMm: ellipse.diameterMm,
+ measuredMm: (rMeasured * 2) / scale,
+ residuals: [iso, skew, diaRes],
+ details: `${prefix}iso ${(iso * 100).toFixed(2)}%, skew ${(skew * 100).toFixed(2)}%, dia ${(diaRes * 100).toFixed(2)}%`,
+ }
+}
+
+function buildReports(
+ datums: Datum[],
+ primary: Primary,
+ H: Mat3,
+ scale: number,
+): { reports: DatumReport[]; rmsPercent: number } {
+ const raw: RawReport[] = datums.map((d) => {
+ const isPrimary = d === primary.datum
+ if (d.type === "line")
+ return residualForLine(d, H, scale, isPrimary)
+ if (d.type === "rectangle")
+ return residualForRect(d, H, scale, isPrimary)
+ return residualForEllipse(d, H, scale, isPrimary)
+ })
+ const reports: DatumReport[] = raw.map((r) => {
+ const rms =
+ Math.sqrt(
+ r.residuals.reduce((s, x) => s + x * x, 0) /
+ Math.max(1, r.residuals.length),
+ ) * 100
+ return {
+ label: r.label,
+ type: r.type,
+ expectedMm: r.expectedMm,
+ measuredMm: r.measuredMm,
+ errorPercent: rms,
+ details: r.details,
+ }
+ })
+ let sumSq = 0
+ let n = 0
+ for (const r of raw) {
+ for (const x of r.residuals) {
+ sumSq += x * x
+ n++
+ }
+ }
+ const rmsPercent = n > 0 ? Math.sqrt(sumSq / n) * 100 : 0
+ return { reports, rmsPercent }
+}
+
+// ─── Public entry point ─────────────────────────────────────────────────────
+
+interface SolverResult {
+ /** 3×3 homography (row-major) mapping source-image px → output px. */
+ H: Mat3
+ /** Label of the datum used as the gauge reference (useful for UI). */
+ primaryLabel: string
+ /** Type of the primary, so callers can report it meaningfully. */
+ primaryType: DatumType
+ iterations: number
+ reports: DatumReport[]
+ rmsPercent: number
+}
+
+function warmStartH(primary: Primary, scale: number): Mat3 {
+ const anchors = primaryAnchors(primary, scale)
+ // Guaranteed 4 anchors for any primary kind
+ const src = cv.matFromArray(
+ 4,
+ 1,
+ cv.CV_32FC2,
+ anchors.flatMap((c) => [c.src.x, c.src.y]),
+ )
+ const dst = cv.matFromArray(
+ 4,
+ 1,
+ cv.CV_32FC2,
+ anchors.flatMap((c) => [c.dst.x, c.dst.y]),
+ )
+ try {
+ const M = cv.getPerspectiveTransform(src, dst)
+ const h = normalized(readMat3x3(M))
+ M.delete()
+ return h
+ } finally {
+ src.delete()
+ dst.delete()
+ }
+}
+
+export function solveHomographyForDatums(
+ datums: Datum[],
+ scale: number,
+): SolverResult {
+ const primary = pickPrimary(datums)
+ let H = warmStartH(primary, scale)
+
+ let iterations = 0
+ let Hprev: Mat3 | null = null
+ let oscillationWarned = false
+ for (let iter = 0; iter < MAX_OUTER_ITERS; iter++) {
+ const corrs = buildCorrespondences(datums, primary, H, scale)
+ const Hnew = solveHomography(corrs)
+ iterations = iter + 1
+ if (!Hnew) break
+ const delta = relativeMaxDiff(H, Hnew)
+
+ // Period-2 oscillation detection: if we're closer to two steps ago
+ // than to one step ago, the outer loop is cycling between two H
+ // values rather than converging. Alternating minimisation gives no
+ // monotone-decrease guarantee (no coherent global objective), so
+ // this can happen with adversarial datum combinations.
+ if (
+ !oscillationWarned &&
+ Hprev &&
+ iter >= 3 &&
+ delta > CONVERGENCE_TOL
+ ) {
+ const deltaSkip = relativeMaxDiff(Hprev, Hnew)
+ if (deltaSkip < delta * 0.25) {
+ console.warn(
+ `[solver] outer loop appears to be oscillating (iter=${String(iter + 1)}, δ=${delta.toExponential(2)}, δ_skip=${deltaSkip.toExponential(2)}). Result may not be optimal.`,
+ )
+ oscillationWarned = true
+ }
+ }
+
+ Hprev = H
+ H = Hnew
+ if (delta < CONVERGENCE_TOL) break
+ }
+
+ const { reports, rmsPercent } = buildReports(datums, primary, H, scale)
+ return {
+ H,
+ primaryLabel: primaryLabel(primary),
+ primaryType: primary.datum.type,
+ iterations,
+ reports,
+ rmsPercent,
+ }
+}
diff --git a/src/stores/app.ts b/src/stores/app.ts
index 269f2d3..e5f77f3 100644
--- a/src/stores/app.ts
+++ b/src/stores/app.ts
@@ -29,7 +29,8 @@ export const useAppStore = defineStore("app", () => {
if (!canProceedToStep3.value || datums.value.length === 0) return false
return datums.value.every((d) => {
if (d.type === "rectangle") return d.widthMm > 0 && d.heightMm > 0
- return d.lengthMm > 0
+ if (d.type === "line") return d.lengthMm > 0
+ return d.diameterMm > 0
})
})
diff --git a/src/types/index.ts b/src/types/index.ts
index 2727c54..5aed621 100644
--- a/src/types/index.ts
+++ b/src/types/index.ts
@@ -22,10 +22,27 @@ export interface LineDatum {
label: string
}
-export type Datum = RectDatum | LineDatum
+export interface EllipseDatum {
+ id: string
+ type: "ellipse"
+ /** Image-space ellipse as 3 free points: center + two conjugate
+ * semi-axis endpoints. axisEndA/axisEndB don't need to be perpendicular;
+ * together with center they give a full 5-DoF ellipse matrix. */
+ center: Point
+ axisEndA: Point
+ axisEndB: Point
+ /** Known real-world diameter of the circle being drawn. */
+ diameterMm: number
+ confidence: 1 | 2 | 3 | 4 | 5
+ label: string
+}
+
+export type Datum = RectDatum | LineDatum | EllipseDatum
export type ConfidenceScore = 1 | 2 | 3 | 4 | 5
+export type DatumType = Datum["type"]
+
export interface ExifData {
make?: string
model?: string
@@ -53,24 +70,26 @@ export interface DeskewInput {
onProgress?: (step: number, total: number, label: string) => void
}
-export interface AxisCorrection {
- ratio: number
- totalWeight: number
-}
-
export interface DatumReport {
label: string
- type: "rectangle" | "line"
- measuredMm: number
+ type: DatumType
+ /** Representative expected dimension in mm (widthMm / lengthMm / diameterMm). */
expectedMm: number
+ /** Representative measured dimension in mm under the solved H. */
+ measuredMm: number
+ /** Overall residual magnitude expressed as a percentage. */
errorPercent: number
- axisContribution: "x" | "y" | "both"
+ /** Free-form breakdown for debugging (e.g. "iso 0.2%, skew 0.1%, dia 0.8%"). */
+ details: string
}
export interface DeskewDiagnostics {
+ /** Label of the rectangle used to fix the output gauge. */
primaryDatum: string
- xCorrection: AxisCorrection
- yCorrection: AxisCorrection
+ /** Number of outer alternating-minimization iterations the solver ran. */
+ iterations: number
+ /** Final weighted RMS residual across all datums, as a percentage. */
+ finalRMSPercent: number
perDatum: DatumReport[]
outputWidthPx: number
outputHeightPx: number
@@ -91,3 +110,8 @@ export interface RectPreset {
widthMm: number
heightMm: number
}
+
+export interface CirclePreset {
+ label: string
+ diameterMm: number
+}