From 4494d40c333474ad3b7ff33ccd697405e236561e Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Mon, 2 Nov 2026 07:55:19 -0800 Subject: [PATCH] Geometry: the inverse of an invertible projective map is kept however small its relative determinant adjugate_inverse no longer applies the relative determinant test to its result. That test is not symmetric between a map and its inverse: a map whose determinant is 1e-8 of its largest product can have a finite, accurate inverse whose own ratio is near 1e-16, so inverse() was None while is_invertible() was true. The result is now kept when every entry is finite and its determinant, summed from the six products in split form, is not 0. --- crates/geom/src/projective.rs | 59 ++++++++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 12 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index baf057065..2a9ff9093 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -90,13 +90,7 @@ impl Homography { if !self.is_finite() { return false; } - let f = self.0.map(|row| row.map(frexp)); - // The columns each of the six products takes from rows 0, 1 and 2, and its sign. - let cols = [(0, 1, 2, 1.0), (0, 2, 1, -1.0), (1, 0, 2, -1.0), (1, 2, 0, 1.0), (2, 0, 1, 1.0), (2, 1, 0, -1.0)]; - let terms = cols.map(|(a, b, c, sign)| { - let (a, b, c) = (f[0][a], f[1][b], f[2][c]); - (sign * a.0 * b.0 * c.0, a.1 + b.1 + c.1) - }); + let terms = self.det_terms(); let Some(e) = terms.iter().filter(|(m, _)| *m != 0.0).map(|(_, e)| *e).max() else { return false; }; @@ -106,6 +100,18 @@ impl Homography { d.abs() > 1e-15 * largest } + /// The six signed products the determinant sums, each as a mantissa product and an exponent + /// sum. + fn det_terms(&self) -> [(f64, i32); 6] { + let f = self.0.map(|row| row.map(frexp)); + // The columns each of the six products takes from rows 0, 1 and 2, and its sign. + let cols = [(0, 1, 2, 1.0), (0, 2, 1, -1.0), (1, 0, 2, -1.0), (1, 2, 0, 1.0), (2, 0, 1, 1.0), (2, 1, 0, -1.0)]; + cols.map(|(a, b, c, sign)| { + let (a, b, c) = (f[0][a], f[1][b], f[2][c]); + (sign * a.0 * b.0 * c.0, a.1 + b.1 + c.1) + }) + } + /// The same map with every entry divided by one power of two, chosen so that the largest of the /// six products its determinant sums is at least 1/8 and below 4. `None` when every one of /// those products is 0 or an entry is not finite. @@ -128,7 +134,7 @@ impl Homography { /// Inverse; `None` when [`Self::has_significant_det`] counts the map as not invertible. The /// inverse is formed from the map itself when its determinant is finite, and is kept when every - /// entry is finite and it passes [`Self::has_significant_det`]. Otherwise it is the inverse of + /// entry is finite and its own determinant is not 0. Otherwise it is the inverse of /// [`Self::scaled`] (the exact inverse multiplied by one power of two, which is the same map), /// and `None` when that one fails the same checks. pub fn inverse(&self) -> Option { @@ -144,9 +150,10 @@ impl Homography { } /// The adjugate divided by the determinant; `None` when an entry is not finite, the - /// determinant is 0, or the result fails [`Self::has_significant_det`]. Every cofactor and the - /// determinant are held as a mantissa and a power of two, so an entry of the result is lost - /// only when that entry itself is beyond the range of `f64`. + /// determinant is 0, or the result has an entry that is not finite or a determinant of 0 (its + /// six products, each held as a mantissa and a power of two, sum to 0, as when a row underflows + /// to 0). Every cofactor and the determinant are held as a mantissa and a power of two, so an + /// entry of the result is lost only when that entry itself is beyond the range of `f64`. fn adjugate_inverse(&self) -> Option { if !self.is_finite() { return None; @@ -176,7 +183,7 @@ impl Homography { } let out = Self(out); // H⁻¹·(x, y, 1) = (u, v, 1)/w_H, so points in front of the camera keep w > 0. - out.has_significant_det().then_some(out) + (out.is_finite() && split_sum(&out.det_terms()).0 != 0.0).then_some(out) } /// `self ∘ other`: `other` first. @@ -519,6 +526,34 @@ mod tests { } } + #[test] + fn an_invertible_map_has_an_inverse_at_any_scale() { + let times = |h: Homography, f: f64| Homography(h.0.map(|row| row.map(|v| v * f))); + // Determinant 1e-17 against a largest product of 1e-9 (and 1e-8 at 1e3 times the map): its + // inverse is finite, though that inverse's own determinant is near 1e-16 of its largest + // product. + let skew = Homography([[-1000.0, -1e-8, -1e-5], [0.0, 1e-4, -1e7], [1e-8, 0.0, 1e-8]]); + let exact = [1e5, 10.0, 1e16 + 1e8, -1e16, -1e12 + 1e4, -1e27, -1e5, -10.0, -1e16]; + for f in [1.0, 1e3] { + let s = times(skew, f); + assert_eq!(Homography::from_array(s.to_array()), Some(s), "{f}"); + let inv = s.inverse().unwrap_or_else(|| panic!("{f}")); + for (v, w) in inv.to_array().iter().zip(exact) { + assert!((v * f - w).abs() <= 1e-9 * w.abs(), "{f}: {inv:?}"); + } + // The map times its inverse is the identity, to within the rounding of the products + // each entry sums (they cancel by up to 1e23 in the second row). + let id = s.then_after(&inv); + for (r, row) in id.0.iter().enumerate() { + for (k, v) in row.iter().enumerate() { + let size: f64 = (0..3).map(|i| (s.0[r][i] * inv.0[i][k]).abs()).sum(); + let want = if r == k { 1.0 } else { 0.0 }; + assert!((v - want).abs() <= 1e-12 * size, "{f}: {id:?}"); + } + } + } + } + #[test] fn a_product_that_overflows_is_none() { // Each factor is finite and invertible, but the product's first entry is 1e400.