From e07ee9d20427e635a2a921655b99f6177c2c3379 Mon Sep 17 00:00:00 2001
From: Samuel Prevost <samuel@infinic.com>
Date: Fri, 24 Apr 2026 18:22:42 +0200
Subject: [PATCH] fix(solver): ellipse diagnostic used wrong H direction
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

residualForEllipse computed C = H^T·E·H, but H maps image → output, so
the correct output-space conic is C = H^{-T}·E·H^{-1}. The forward form
represents a geometrically meaningless conic and produced nonsensical
numbers in the per-datum table (a 210mm circle was reported as 2046mm).
The deskewed image itself was correct — only the diagnostic was wrong,
because these residuals are display-only and don't feed back into the
solver. Now invert H once and compute the conic in output space.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
---
 src/lib/solver.ts | 64 +++++++++++++++++++++++++++++++++++++++--------
 1 file changed, 54 insertions(+), 10 deletions(-)

diff --git a/src/lib/solver.ts b/src/lib/solver.ts
index a2353a4..999d1c5 100644
--- a/src/lib/solver.ts
+++ b/src/lib/solver.ts
@@ -94,6 +94,36 @@ function relativeMaxDiff(a: Mat3, b: Mat3): number {
     return scale > 0 ? diff / scale : diff
 }
 
+/** Inverse of a 3×3 homography (row-major). Returns null if singular. */
+function invertMat3(h: Mat3): Mat3 | null {
+    const a = h[0]
+    const b = h[1]
+    const c = h[2]
+    const d = h[3]
+    const e = h[4]
+    const f = h[5]
+    const g = h[6]
+    const hh = h[7]
+    const i = h[8]
+    const det =
+        a * (e * i - f * hh) -
+        b * (d * i - f * g) +
+        c * (d * hh - e * g)
+    if (Math.abs(det) < 1e-20) return null
+    const invDet = 1 / det
+    return [
+        (e * i - f * hh) * invDet,
+        (c * hh - b * i) * invDet,
+        (b * f - c * e) * invDet,
+        (f * g - d * i) * invDet,
+        (a * i - c * g) * invDet,
+        (c * d - a * f) * invDet,
+        (d * hh - e * g) * invDet,
+        (b * g - a * hh) * invDet,
+        (a * e - b * d) * invDet,
+    ]
+}
+
 // ─── Geometric helpers ──────────────────────────────────────────────────────
 
 function dist(a: Point, b: Point): number {
@@ -647,19 +677,33 @@ function residualForEllipse(
     scale: number,
     isPrimary: boolean,
 ): RawReport {
-    // Pull the image ellipse back to world space via C = H^T E H
+    // E is the image-space ellipse conic; H maps image → output. The
+    // output-space conic we want to check for circularity is therefore
+    //     C = H^{-T} · E · H^{-1}
+    // (points q on C iff H^{-1}·q lands on E). Using H^T·E·H computes the
+    // wrong conic and produces meaningless diagnostic numbers.
     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 Hi = invertMat3(H)
+    const zeroReport: RawReport = {
+        label: ellipse.label,
+        type: "ellipse",
+        expectedMm: ellipse.diameterMm,
+        measuredMm: 0,
+        residuals: [0, 0, 0],
+        details: "(H singular, cannot compute)",
+    }
+    if (!Hi) return zeroReport
+    const HiT = [
+        [Hi[0], Hi[3], Hi[6]],
+        [Hi[1], Hi[4], Hi[7]],
+        [Hi[2], Hi[5], Hi[8]],
     ]
     const Hm = [
-        [H[0], H[1], H[2]],
-        [H[3], H[4], H[5]],
-        [H[6], H[7], H[8]],
+        [Hi[0], Hi[1], Hi[2]],
+        [Hi[3], Hi[4], Hi[5]],
+        [Hi[6], Hi[7], Hi[8]],
     ]
-    // C = HT · E · Hm
+    // C = HiT · E · Hm    (= H^{-T} E H^{-1})
     const EH: number[][] = [
         [0, 0, 0],
         [0, 0, 0],
@@ -683,7 +727,7 @@ function residualForEllipse(
         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)
+                s += (HiT[i]?.[k] ?? 0) * (EH[k]?.[j] ?? 0)
             }
             ;(C[i] as number[])[j] = s
         }