From 3d345ff95bb280adac3d24555c7dba07c481af44 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Fri, 30 Oct 2026 23:07:50 -0700 Subject: [PATCH] Geometry: a perspective point that overflows maps to nothing instead of infinity Homography::apply checked that w is in front of the horizon but not the quotients, so a point with w just above 1e-9 and a large x or y, or a NaN w, mapped to an infinite or NaN point. apply now returns None unless the mapped point is finite. inverse, conjugated and affine_at return None when their result has an entry that is not finite. --- crates/geom/src/projective.rs | 59 ++++++++++++++++++++++++++++++----- 1 file changed, 52 insertions(+), 7 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 349ddb80f..6264df67d 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -40,7 +40,8 @@ impl Homography { [m[0][0], m[0][1], m[0][2], m[1][0], m[1][1], m[1][2], m[2][0], m[2][1], m[2][2]] } - /// Map a point; `None` when it lands on or beyond the horizon (w ≤ 0). + /// Map a point; `None` when it lands on or beyond the horizon (w ≤ 0), or when it does not map + /// to a finite page point. pub fn apply(&self, p: Point) -> Option { let m = &self.0; let x = m[0][0] * p.x + m[0][1] * p.y + m[0][2]; @@ -49,7 +50,8 @@ impl Homography { if w <= 1e-9 || !x.is_finite() || !y.is_finite() { return None; } - Some(Point::new(x / w, y / w)) + let q = Point::new(x / w, y / w); + q.is_finite().then_some(q) } pub fn det(&self) -> f64 { @@ -58,7 +60,8 @@ impl Homography { + m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) } - /// Inverse. + /// Inverse; `None` when the map is not invertible or its inverse has an entry that is not + /// finite. pub fn inverse(&self) -> Option { let m = &self.0; let d = self.det(); @@ -77,6 +80,9 @@ impl Homography { out[r][k] = v / d; } } + if out.iter().flatten().any(|v| !v.is_finite()) { + return None; + } // H⁻¹·(x, y, 1) = (u, v, 1)/w_H, so points in front of the camera keep w > 0. Some(Self(out)) } @@ -94,16 +100,18 @@ impl Homography { } /// `self` seen through `a`: the map `a ∘ self ∘ a⁻¹` (what `self` becomes when everything it - /// maps from and to is moved by `a`). `None` when `a` isn't invertible. + /// maps from and to is moved by `a`). `None` when `a` isn't invertible or the result has an + /// entry that is not finite. pub fn conjugated(&self, a: Affine) -> Option { if a.determinant().abs() < 1e-15 { return None; } - Some(Self::from_affine(a).then_after(self).then_after(&Self::from_affine(a.inverse()))) + let c = Self::from_affine(a).then_after(self).then_after(&Self::from_affine(a.inverse())); + c.0.iter().flatten().all(|v| v.is_finite()).then_some(c) } /// The affine map that matches `self` to first order at `p` (its derivative there), `None` - /// beyond the horizon. + /// beyond the horizon or when the derivative is not finite. pub fn affine_at(&self, p: Point) -> Option { let m = &self.0; let w = m[2][0] * p.x + m[2][1] * p.y + m[2][2]; @@ -111,7 +119,8 @@ impl Homography { // d(x/w) = (dx − (x/w)·dw) / w, and the same for y. let j = |r: usize, c: usize| (m[r][c] - [q.x, q.y][r] * m[2][c]) / w; let lin = Affine::new([j(0, 0), j(1, 0), j(0, 1), j(1, 1), 0.0, 0.0]); - Some(Affine::translate(q.to_vec2()) * lin * Affine::translate(-p.to_vec2())) + let out = Affine::translate(q.to_vec2()) * lin * Affine::translate(-p.to_vec2()); + out.is_finite().then_some(out) } /// The bounding box of the image of `r` (lines stay lines: its corners bound it), `None` when @@ -175,6 +184,42 @@ mod tests { } } + #[test] + fn a_point_that_overflows_maps_to_nothing() { + // w = 1e-8 is in front of the horizon, but 1e305 / 1e-8 overflows. + let near = Homography([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1e-8]]); + assert_eq!(near.apply(Point::new(1e305, 0.0)), None); + assert_eq!(near.apply(Point::new(0.0, -1e305)), None); + assert_eq!(near.map_rect_bbox(Rect::new(0.0, 0.0, 1e305, 1.0)), None); + assert!(near.apply(Point::new(2.0, 3.0)).unwrap().distance(Point::new(2e8, 3e8)) < 1e-6); + // A NaN w is not caught by the horizon test. + let nan = Homography([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, f64::NAN]]); + assert_eq!(nan.apply(Point::new(1.0, 1.0)), None); + } + + #[test] + fn an_inverse_that_overflows_is_none() { + let big = Homography([[1e300, 0.0, 0.0], [0.0, 1e300, 0.0], [0.0, 0.0, 1e-310]]); + assert_eq!(big.inverse(), None); + assert!(h().inverse().is_some()); + } + + #[test] + fn a_derivative_that_overflows_is_none() { + // The point maps to 1e300, but the derivative there is 1e302 / 1e-8. + let steep = Homography([[1e302, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1e-8]]); + let p = Point::new(1e-10, 0.0); + assert!(steep.apply(p).is_some()); + assert_eq!(steep.affine_at(p), None); + } + + #[test] + fn a_conjugate_that_overflows_is_none() { + let far = Homography([[1.0, 0.0, 1e200], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); + assert_eq!(far.conjugated(Affine::scale(1e150)), None); + assert!(far.conjugated(Affine::scale(2.0)).is_some()); + } + #[test] fn rect_bbox_is_the_corners_box() { let b = h().map_rect_bbox(Rect::new(0.0, 0.0, 10.0, 20.0)).unwrap();