`homography2d`: DLT normal equations accumulated in Float32, capping achievable precision

[kalwalt]

Summary

homography2d builds the 9×9 DLT normal matrix and solves for its smallest eigenvector entirely in single precision:

this.mLtL = new matrix_t(9, 9, JSFEAT_CONSTANTS.F32_t | JSFEAT_CONSTANTS.C1_t);
this.Evec = new matrix_t(9, 9, JSFEAT_CONSTANTS.F32_t | JSFEAT_CONSTANTS.C1_t);

and motion_model.T0 / T1 are F32_t as well.

OpenCV's HomographyEstimatorCallback::runKernel (modules/calib3d/src/fundam.cpp) uses doubles throughout:

double LtL[9][9], W[9][1], V[9][9];
Mat _LtL( 9, 9, CV_64F, &LtL[0][0] );
...
double invHnorm[9] = { 1./sm.x, 0, cm.x, 0, 1./sm.y, cm.y, 0, 0, 1 };

Two compounding problems:

1. Normal equations square the condition number. Forming LᵀL instead of solving on L directly costs roughly half the available significant digits. f32 has ~7 decimal digits to begin with, so this is expensive.
2. Every accumulation rounds. LtL[0] += X * X reads a Float32Array slot, adds in f64, and stores back rounded to f32 — 36 upper-triangle accumulators × N points, each rounded on every step.
   The structure of the hand-unrolled accumulation itself is correct — I checked all 45 upper-triangle slots against OpenCV's for(j) for(k=j) loop, and every omitted index is a genuine structural zero (Lx[0]*Lx[3] + Ly[0]*Ly[3] = X*0 + 0*X). This is purely a dtype issue.

Impact

Measured on webarkit/webarkit's cv-backend-jsfeatnext adapter: on a noise-free, mathematically exact 12-point translation, repeated ransac() calls returned homographies differing by up to ~0.36px in reprojection. The refit proposed in #185 will shrink that spread a lot but cannot remove it — f32 accumulation is the floor underneath it. Any sub-pixel acceptance criterion in #185 depends on this issue being fixed too.

Proposed change

eigenVV is already dtype-generic: it derives dt from A.type and sizes its scratch buffers with <<3 (8 bytes/element, f64-safe), and matrix_t already supports F64_t. matmath.multiply_3x3 operates on .data without dtype assumptions. So the change looks small:

• mLtL, Evec → F64_t | C1_t
• motion_model.T0, T1 → F64_t | C1_t (and AtA/AtB in affine2d's path, for consistency)
• Use an internal F64 scratch for the 3×3 model and copy into the caller's model once at the end, so the two multiply_3x3 denormalisation steps don't run at the caller's dtype (callers construct model as F32_t — see the existing tests)
• Guard the final normalisation the way OpenCV does. We currently do x = 1.0 / md[8] unconditionally; OpenCV uses scaleFor, which falls back to 1.0 when |x| is below float epsilon:

  static inline double scaleFor(double x){
      return (std::fabs(x) > std::numeric_limits<float>::epsilon()) ? 1./x : 1.;
  }

A near-zero md[8] (pure-perspective degenerate) currently produces Infinity silently.

Acceptance criteria

• DLT accumulation and eigen-decomposition run in f64
• Final normalisation guarded against md[8] ≈ 0
• tests/parity/linalg.test.ts and matmath.test.ts still pass (eigenVV is shared)
• tests/parity/motion_estimator.test.ts: the toBeCloseTo(..., 5) comparison against the f32 jsfeat oracle will now diverge in our favour — expect this and document it, or compare against a higher-precision reference instead
• Benchmark bench/motion_model.bench.ts before/after — f64 arithmetic in JS is normally not slower than f32, but Float64Array doubles the memory traffic on the 9×9, so confirm rather than assume

Out of scope

• Replacing normal equations with an SVD on L directly. linalg.svd_decompose exists, and it would be the numerically right answer, but it is a bigger change with WASM-portability implications for the downstream no_std work. Worth its own issue if f64 turns out not to be enough.
• affine2d's lu_solve path — same dtype question, but it is a 6×6 well-conditioned system and has not been measured as a problem.

Related

• motion_estimator.ransac: doc comment promises a refit over all inliers; jsfeatNext has no layer that performs one #185 (RANSAC refit) — the two together determine the achievable precision; neither is sufficient alone
• webarkit/webarkit#8 (loosened adapter test tolerance)

[kalwalt]

Fixed via jsfeatNext#192 (merged into dev): homography2d's DLT solve (mLtL, Evec, T0/T1, and a new F64 model64 scratch) now runs entirely in F64, only rounding to the caller's model dtype once at the very end. Guarded the final 1.0/md[8] normalization; a near-zero md[8] now reports degeneracy (return 0) rather than fabricating a corrupted model — this codebase's homography2d.error() hardcodes h33 as the literal constant 1.0, so a model that can't be scaled to h33=1 can't be represented correctly here, unlike OpenCV's own scaleFor() fallback which just leaves h33 unscaled. tests/parity/motion_estimator.test.ts's homography2d cases now compare inlier masks only (not model coefficients); tests/divergences.test.ts registers and measures the intentional divergence from jsfeat's F32 answer. Fixed a related bug found by Qodo's review during development (the near-zero fallback originally forced h33=1 without an equivalent homogeneous rescale).