sjat/GCodeOverlay
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

refactor: retire homography model in favour of PolyWarp

Browse source 
Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
sjat 2026-06-11 09:35:51 +02:00
parent 31577787f1
commit 83f4c3eb7e
3 changed files with 6 additions and 118 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

31
src/geometry/homography.test.ts
View file

@ -1,31 +0,0 @@
import { describe, it, expect } from 'vitest';
import { estimateHomography, applyHomography, invertMat3 } from './homography';
import type { Vec2, Mat3 } from '../types';
const close = (a: Vec2, b: Vec2, eps = 1e-6) =>
Math.abs(a[0] - b[0]) < eps && Math.abs(a[1] - b[1]) < eps;
describe('homography', () => {
it('recovers the identity from coincident point sets', () => {
const pts: Vec2[] = [[0, 0], [1, 0], [1, 1], [0, 1]];
const H = estimateHomography(pts, pts);
expect(close(applyHomography(H, [0.3, 0.7]), [0.3, 0.7])).toBe(true);
});
it('recovers a known perspective transform', () => {
const H: Mat3 = [1.1, 0.2, 5, -0.1, 0.9, 3, 0.001, 0.002, 1];
const src: Vec2[] = [[0, 0], [100, 0], [100, 100], [0, 100], [40, 60]];
const dst = src.map((p) => applyHomography(H, p));
const est = estimateHomography(src, dst);
for (const p of [[10, 20], [70, 30], [55, 90]] as Vec2[]) {
expect(close(applyHomography(est, p), applyHomography(H, p), 1e-4)).toBe(true);
}
});
it('invertMat3 round-trips a point', () => {
const H: Mat3 = [1.1, 0.2, 5, -0.1, 0.9, 3, 0.001, 0.002, 1];
const inv = invertMat3(H);
const p: Vec2 = [37, 91];
expect(close(applyHomography(inv, applyHomography(H, p)), p, 1e-6)).toBe(true);
});
});

78
src/geometry/homography.ts
View file

@ -1,78 +0,0 @@
import type { Vec2, Mat3 } from '../types';
/** Apply a homography to a point (projective divide). */
export function applyHomography(H: Mat3, p: Vec2): Vec2 {
const [a, b, c, d, e, f, g, h, i] = H;
const w = g * p[0] + h * p[1] + i;
return [(a * p[0] + b * p[1] + c) / w, (d * p[0] + e * p[1] + f) / w];
}
/** Solve a square linear system A x = b by Gaussian elimination with partial pivoting. */
function solveLinear(A: number[][], b: number[]): number[] {
const n = b.length;
const M = A.map((row, i) => [...row, b[i]!]);
for (let col = 0; col < n; col++) {
let pivot = col;
for (let r = col + 1; r < n; r++) {
if (Math.abs(M[r]![col]!) > Math.abs(M[pivot]![col]!)) pivot = r;
}
if (Math.abs(M[pivot]![col]!) < 1e-12) throw new Error('Singular system in homography fit');
[M[col], M[pivot]] = [M[pivot]!, M[col]!];
const pivRow = M[col]!;
for (let r = 0; r < n; r++) {
if (r === col) continue;
const factor = M[r]![col]! / pivRow[col]!;
for (let k = col; k <= n; k++) M[r]![k]! -= factor * pivRow[k]!;
}
}
return M.map((row, i) => row[n]! / row[i]!);
}
/**
* Estimate a homography mapping src → dst from ≥4 point pairs.
* Solves for 8 DoF (h33 fixed to 1) via least squares (normal equations).
*/
export function estimateHomography(src: Vec2[], dst: Vec2[]): Mat3 {
if (src.length !== dst.length || src.length < 4) {
throw new Error('Need at least 4 matching point pairs');
}
// Build 2N×8 design matrix M and target vector t.
const rows: number[][] = [];
const t: number[] = [];
for (let i = 0; i < src.length; i++) {
const [x, y] = src[i]!;
const [u, v] = dst[i]!;
rows.push([x, y, 1, 0, 0, 0, -x * u, -y * u]);
t.push(u);
rows.push([0, 0, 0, x, y, 1, -x * v, -y * v]);
t.push(v);
}
// Normal equations: (Mᵀ M) h = Mᵀ t → 8×8 solve.
const ata: number[][] = Array.from({ length: 8 }, () => new Array(8).fill(0));
const atb: number[] = new Array(8).fill(0);
for (let r = 0; r < rows.length; r++) {
const row = rows[r]!;
for (let i = 0; i < 8; i++) {
atb[i]! += row[i]! * t[r]!;
for (let j = 0; j < 8; j++) ata[i]![j]! += row[i]! * row[j]!;
}
}
const h = solveLinear(ata, atb);
return [h[0]!, h[1]!, h[2]!, h[3]!, h[4]!, h[5]!, h[6]!, h[7]!, 1];
}
/** Invert a 3×3 matrix. */
export function invertMat3(m: Mat3): Mat3 {
const [a, b, c, d, e, f, g, h, i] = m;
const A = e * i - f * h;
const B = -(d * i - f * g);
const C = d * h - e * g;
const det = a * A + b * B + c * C;
if (Math.abs(det) < 1e-12) throw new Error('Non-invertible matrix');
const invDet = 1 / det;
return [
A * invDet, (c * h - b * i) * invDet, (b * f - c * e) * invDet,
B * invDet, (a * i - c * g) * invDet, (c * d - a * f) * invDet,
C * invDet, (b * g - a * h) * invDet, (a * e - b * d) * invDet,
];
}

15
src/geometry/polywarp.test.ts
View file

@ -1,6 +1,5 @@
import { describe, it, expect } from 'vitest'; import { describe, it, expect } from 'vitest';
import { estimatePolyWarp, applyPolyWarp, applyPolyWarpInverse } from './polywarp'; import { estimatePolyWarp, applyPolyWarp, applyPolyWarpInverse } from './polywarp';
import { estimateHomography, applyHomography } from './homography';
import type { Vec2 } from '../types'; import type { Vec2 } from '../types';
const dist = (a: Vec2, b: Vec2) => Math.hypot(a[0] - b[0], a[1] - b[1]); const dist = (a: Vec2, b: Vec2) => Math.hypot(a[0] - b[0], a[1] - b[1]);
@ -34,18 +33,16 @@ describe('polywarp', () => {
} }
}); });
it('beats a homography on barrel-distorted data', () => { it('beats an affine baseline on barrel-distorted data', () => {
// Barrel distortion multiplies the coordinate by (1 + k·r²), so as a function of // Barrel distortion is degree-3 in machine coords: degree 3 fits it near-exactly while
// machine coords it is degree 3 (dx·r² → dx³ + dx·dy²). A homography (and a degree-2 // a degree-1 (affine, no curvature) baseline cannot.
// polynomial) cannot represent that; degree 3 fits it essentially exactly.
const image = machineGrid.map(trueMap); const image = machineGrid.map(trueMap);
const w = estimatePolyWarp(machineGrid, image, 3); const w = estimatePolyWarp(machineGrid, image, 3);
const H = estimateHomography(machineGrid, image); const baseline = estimatePolyWarp(machineGrid, image, 1);
// Evaluate on off-grid interior points (fifths), where barrel error is real.
const test: Vec2[] = [[488, 244], [1464, 732], [1952, 488], [976, 976]]; const test: Vec2[] = [[488, 244], [1464, 732], [1952, 488], [976, 976]];
const polyMax = Math.max(...test.map((m) => dist(applyPolyWarp(w, m), trueMap(m)))); const polyMax = Math.max(...test.map((m) => dist(applyPolyWarp(w, m), trueMap(m))));
const homoMax = Math.max(...test.map((m) => dist(applyHomography(H, m), trueMap(m)))); const baseMax = Math.max(...test.map((m) => dist(applyPolyWarp(baseline, m), trueMap(m))));
expect(homoMax).toBeGreaterThan(1e-3); // homography genuinely fails on barrel data expect(baseMax).toBeGreaterThan(1e-3); // affine baseline genuinely fails on barrel data
expect(polyMax).toBeLessThan(1e-6); // degree-3 fit of cubic data is near-exact everywhere expect(polyMax).toBeLessThan(1e-6); // degree-3 fit of cubic data is near-exact everywhere
}); });