From 3d345ff95bb280adac3d24555c7dba07c481af44 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Fri, 30 Oct 2026 23:07:50 -0700 Subject: [PATCH 1/5] 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(); From db38e2e1144dd07e419f4668c35e70d41726e721 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Fri, 30 Oct 2026 23:27:44 -0700 Subject: [PATCH 2/5] Perspective: a projection that a transform overflows is scaled down or removed When moving type or a symbol instance in perspective makes the conjugated projection overflow, Homography::conjugated now forms the product from factors scaled to a largest entry of 1 and returns it at that scale, which maps every point to the same place. When even the scaled result is not invertible, the transform removes the projection and the art is drawn flat, instead of keeping the old projection while the flat art moves. --- crates/doc/src/node/tests.rs | 29 ++++++++++++++++++++++++ crates/doc/src/perspective.rs | 8 ++++--- crates/geom/src/projective.rs | 42 ++++++++++++++++++++++++++++++----- 3 files changed, 71 insertions(+), 8 deletions(-) diff --git a/crates/doc/src/node/tests.rs b/crates/doc/src/node/tests.rs index dc241841f..b4cc7b072 100644 --- a/crates/doc/src/node/tests.rs +++ b/crates/doc/src/node/tests.rs @@ -339,6 +339,35 @@ fn transform_group_recurses() { assert_eq!(child.appearance.stroke_width(), 2.0); } +/// A symbol instance in perspective keeps its picture moved by a transform, and loses its +/// projection when no invertible map of finite numbers draws the moved picture. +#[test] +fn a_projection_follows_the_transform_or_is_removed() { + use vectorcraft_geom::{Homography, Point}; + let instance = |h: Homography| { + let mut n = Node::new(NodeId(1), NodeKind::SymbolInstance { symbol: "S".into(), xf: Affine::IDENTITY }); + n.extras_mut().perspective = + Some(Box::new(crate::PerspectiveAttachment { projection: Some(h.to_array()), ..crate::PerspectiveAttachment::new("left", 0.0) })); + n + }; + let projection = |n: &Node| n.extras().perspective.as_deref().and_then(|p| p.projection); + // A move by 10 scaled by 1e307: moved by a scale by 10, the product overflows unless scaled down. + let big = Homography([[1e307, 0.0, 1e308], [0.0, 1e307, 0.0], [0.0, 0.0, 1e307]]); + let mut n = instance(big); + n.transform(Affine::scale(10.0), false); + let p = Point::new(3.0, 7.0); + assert!(n.projection().unwrap().apply(p).unwrap().distance(Point::new(103.0, 7.0)) < 1e-9); + // A move by 1e200 moved by a scale by 1e150 sends every point to about (1e350, 0). + let far = Homography([[1.0, 0.0, 1e200], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); + let mut n = instance(far); + n.transform(Affine::scale(1e150), false); + assert_eq!(projection(&n), None); + // A transform that is not invertible leaves the projection as it is. + let mut n = instance(far); + n.transform(Affine::scale(0.0), false); + assert_eq!(projection(&n), Some(far.to_array())); +} + #[test] fn unplaced_gradients_pin_only_when_a_refit_would_differ() { use vectorcraft_color::{Gradient, GradientGeom, GradientKind, GradientPaint, Paint}; diff --git a/crates/doc/src/perspective.rs b/crates/doc/src/perspective.rs index 4de3823db..6ec18bcfa 100644 --- a/crates/doc/src/perspective.rs +++ b/crates/doc/src/perspective.rs @@ -58,12 +58,14 @@ impl Node { } } - /// Keep a projection drawing the same picture moved by `a` (the flat art moves by `a` too). + /// Keep a projection drawing the same picture moved by `a` (the flat art moves by `a` too). An + /// invertible `a` that the projection cannot follow removes it, and the art is drawn flat. pub(crate) fn transform_projection(&mut self, a: vectorcraft_geom::Affine) { if let Some(p) = self.extras_if_any().and_then(|e| e.perspective.as_deref_mut()) - && let Some(h) = p.homography().and_then(|h| h.conjugated(a)) + && let Some(h) = p.homography() + && a.determinant().abs() >= 1e-15 { - p.projection = Some(h.to_array()); + p.projection = h.conjugated(a).map(|h| h.to_array()); } } } diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 6264df67d..e7754c131 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -100,14 +100,31 @@ 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 or the result has an - /// entry that is not finite. + /// maps from and to is moved by `a`). When the product overflows, it is formed from factors + /// scaled to a largest entry of 1 and returned at that scale (the same map). `None` when `a` + /// isn't invertible, or when the product overflows and the scaled result is not invertible. pub fn conjugated(&self, a: Affine) -> Option { if a.determinant().abs() < 1e-15 { return None; } - 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) + let (to, from) = (Self::from_affine(a), Self::from_affine(a.inverse())); + let c = to.then_after(self).then_after(&from); + if c.is_finite() { + return Some(c); + } + let c = to.normalized()?.then_after(&self.normalized()?).then_after(&from.normalized()?).normalized()?; + (c.det().abs() > 1e-15).then_some(c) + } + + fn is_finite(&self) -> bool { + self.0.iter().flatten().all(|v| v.is_finite()) + } + + /// `self` divided by its largest entry magnitude, a positive scale that maps every point to + /// the same place; `None` when an entry is not finite or every entry is 0. + fn normalized(&self) -> Option { + let k = self.0.iter().flatten().fold(0.0_f64, |k, v| k.max(v.abs())); + (self.is_finite() && k > 0.0).then(|| Self(self.0.map(|row| row.map(|v| v / k)))) } /// The affine map that matches `self` to first order at `p` (its derivative there), `None` @@ -214,10 +231,25 @@ mod tests { } #[test] - fn a_conjugate_that_overflows_is_none() { + fn a_conjugate_that_overflows_is_scaled_down() { + // A move by 10, scaled by 1e307: seen through a scale by 10 it moves by 100, but the + // product's corner entry is 1e309. + let big = Homography([[1e307, 0.0, 1e308], [0.0, 1e307, 0.0], [0.0, 0.0, 1e307]]); + let t = Affine::scale(10.0); + let c = big.conjugated(t).unwrap(); + assert!(c.0.iter().flatten().all(|v| v.is_finite() && v.abs() <= 1.0)); + assert_eq!(Homography::from_array(c.to_array()), Some(c)); + for p in [Point::new(3.0, 7.0), Point::new(-4.0, 1.5)] { + assert!(c.apply(t * p).unwrap().distance(t * big.apply(p).unwrap()) < 1e-9); + assert!(c.apply(p).unwrap().distance(p + kurbo::Vec2::new(100.0, 0.0)) < 1e-9); + } + // A move by 1e200 seen through a scale by 1e150 sends every point to about (1e350, 0), + // which no invertible map of finite numbers draws. 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()); + // A product that does not overflow is kept as it is. + assert_eq!(far.conjugated(Affine::scale(2.0)).unwrap().0[0][2], 2e200); } #[test] From 375b413ba6aa3d27fdef6f95346072f28940be20 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Fri, 30 Oct 2026 23:46:02 -0700 Subject: [PATCH 3/5] Perspective: a projection that a transform overflows keeps a w row of 1 When conjugating a projection overflows, the scaled-down result is now divided by the largest entry in its w row. A valid conjugate is no longer dropped by the determinant cutoff at the arbitrary scale, and every point of a kept projection passes the horizon test in apply instead of none. --- crates/geom/src/projective.rs | 32 ++++++++++++++++++++++++++++---- 1 file changed, 28 insertions(+), 4 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index e7754c131..94831fb5c 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -101,8 +101,10 @@ 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`). When the product overflows, it is formed from factors - /// scaled to a largest entry of 1 and returned at that scale (the same map). `None` when `a` - /// isn't invertible, or when the product overflows and the scaled result is not invertible. + /// scaled to a largest entry of 1 and returned divided by the largest magnitude in its `w` row + /// (the same map, with a `w` row whose largest entry is 1). `None` when `a` isn't invertible, + /// or when the product overflows and the result at that scale has an entry that is not finite + /// or is not invertible. pub fn conjugated(&self, a: Affine) -> Option { if a.determinant().abs() < 1e-15 { return None; @@ -113,7 +115,12 @@ impl Homography { return Some(c); } let c = to.normalized()?.then_after(&self.normalized()?).then_after(&from.normalized()?).normalized()?; - (c.det().abs() > 1e-15).then_some(c) + let k = c.0[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); + if k <= 0.0 { + return None; + } + let c = Self(c.0.map(|row| row.map(|v| v / k))); + (c.is_finite() && c.det().abs() > 1e-15).then_some(c) } fn is_finite(&self) -> bool { @@ -237,12 +244,29 @@ mod tests { let big = Homography([[1e307, 0.0, 1e308], [0.0, 1e307, 0.0], [0.0, 0.0, 1e307]]); let t = Affine::scale(10.0); let c = big.conjugated(t).unwrap(); - assert!(c.0.iter().flatten().all(|v| v.is_finite() && v.abs() <= 1.0)); + assert!(c.is_finite()); + assert_eq!(c.0[2], [0.0, 0.0, 1.0]); assert_eq!(Homography::from_array(c.to_array()), Some(c)); for p in [Point::new(3.0, 7.0), Point::new(-4.0, 1.5)] { assert!(c.apply(t * p).unwrap().distance(t * big.apply(p).unwrap()) < 1e-9); assert!(c.apply(p).unwrap().distance(p + kurbo::Vec2::new(100.0, 0.0)) < 1e-9); } + // A move by 1 scaled by 1e307, seen through a scale by 1e7: a move by 1e7, whose + // determinant is 1 at a w row of 1 but 1e-21 at a largest entry of 1. + let one = Homography([[1e307, 0.0, 1e307], [0.0, 1e307, 0.0], [0.0, 0.0, 1e307]]); + assert!(Homography::from_array(one.to_array()).is_some()); + let c = one.conjugated(Affine::scale(1e7)).unwrap(); + assert_eq!(Homography::from_array(c.to_array()), Some(c)); + let p = Point::new(3.0, -2.0); + assert!(c.apply(p).unwrap().distance(p + kurbo::Vec2::new(1e7, 0.0)) < 1e-6); + // A scale by 1e10 seen through a scale by 1e10: still a scale by 1e10, which every point + // reaches only when w is not left at 1e-10. + let wide = Homography([[1e300, 0.0, 0.0], [0.0, 1e300, 0.0], [0.0, 0.0, 1e290]]); + let c = wide.conjugated(Affine::scale(1e10)).unwrap(); + assert_eq!(Homography::from_array(c.to_array()), Some(c)); + for p in [Point::new(1.0, 2.0), Point::new(-5.0, 0.5)] { + assert!(c.apply(p).unwrap().distance(Point::new(p.x * 1e10, p.y * 1e10)) < 1e-3); + } // A move by 1e200 seen through a scale by 1e150 sends every point to about (1e350, 0), // which no invertible map of finite numbers draws. let far = Homography([[1.0, 0.0, 1e200], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); From c711ad42f0f2372a5bb5ad5d8783c235493a412f Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Fri, 30 Oct 2026 23:52:44 -0700 Subject: [PATCH 4/5] Perspective: a transform's overflow fallback scales each partial product When the conjugate of a projection overflows, it is now formed one factor at a time, and each partial product is divided by its largest entry magnitude. Scaling the affine factors up front let the w row underflow to 0, so a scale by 2 seen through a move by 1e308 returned None and the projection was removed, though its conjugate is finite. --- crates/geom/src/projective.rs | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 94831fb5c..4b17aac60 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -100,11 +100,11 @@ 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`). When the product overflows, it is formed from factors - /// scaled to a largest entry of 1 and returned divided by the largest magnitude in its `w` row - /// (the same map, with a `w` row whose largest entry is 1). `None` when `a` isn't invertible, - /// or when the product overflows and the result at that scale has an entry that is not finite - /// or is not invertible. + /// maps from and to is moved by `a`). When the product overflows, it is formed one factor at a + /// time, each partial product scaled to a largest entry of 1, and returned divided by the + /// largest magnitude in its `w` row (the same map, with a `w` row whose largest entry is 1). + /// `None` when `a` isn't invertible, or when the product overflows and the result at that + /// scale has an entry that is not finite or is not invertible. pub fn conjugated(&self, a: Affine) -> Option { if a.determinant().abs() < 1e-15 { return None; @@ -114,7 +114,7 @@ impl Homography { if c.is_finite() { return Some(c); } - let c = to.normalized()?.then_after(&self.normalized()?).then_after(&from.normalized()?).normalized()?; + let c = to.then_after_normalized(&self.normalized()?)?.then_after_normalized(&from)?; let k = c.0[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); if k <= 0.0 { return None; @@ -134,6 +134,17 @@ impl Homography { (self.is_finite() && k > 0.0).then(|| Self(self.0.map(|row| row.map(|v| v / k)))) } + /// `self ∘ other` divided by its largest entry magnitude; when that product overflows, it is + /// formed from both factors divided by their largest entry magnitudes. `None` when the result + /// has an entry that is not finite or every entry is 0. + fn then_after_normalized(&self, other: &Self) -> Option { + let c = self.then_after(other); + if c.is_finite() { + return c.normalized(); + } + self.normalized()?.then_after(&other.normalized()?).normalized() + } + /// The affine map that matches `self` to first order at `p` (its derivative there), `None` /// beyond the horizon or when the derivative is not finite. pub fn affine_at(&self, p: Point) -> Option { @@ -267,6 +278,16 @@ mod tests { for p in [Point::new(1.0, 2.0), Point::new(-5.0, 0.5)] { assert!(c.apply(p).unwrap().distance(Point::new(p.x * 1e10, p.y * 1e10)) < 1e-3); } + // A scale by 2 seen through a move by 1e308: p goes to 2p - 1e308, though the product's + // corner entry is 2 * -1e308 + 1e308 and the moves scaled to a largest entry of 1 leave + // their w rows at 1e-308. + let two = Homography([[2.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 1.0]]); + let c = two.conjugated(Affine::translate((1e308, 0.0))).unwrap(); + let want = [2.0, 0.0, -1e308, 0.0, 2.0, 0.0, 0.0, 0.0, 1.0]; + for (v, w) in c.to_array().iter().zip(want) { + assert!((v - w).abs() <= 1e-9 * w.abs().max(1.0), "{c:?}"); + } + assert_eq!(Homography::from_array(c.to_array()), Some(c)); // A move by 1e200 seen through a scale by 1e150 sends every point to about (1e350, 0), // which no invertible map of finite numbers draws. let far = Homography([[1.0, 0.0, 1e200], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]); From 8363c815ecb7222f1a04b664d414ad0a3f3a6265 Mon Sep 17 00:00:00 2001 From: dotCooCoo Date: Sat, 31 Oct 2026 00:05:46 -0700 Subject: [PATCH 5/5] Perspective: a transform's overflow fallback keeps each entry's powers of two apart MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit When the conjugate a ∘ self ∘ a⁻¹ overflows, conjugated now forms every entry of the product as a sum of terms whose factors are split into a mantissa and a power of two, and scales the whole product by one power of two from its w row at the end. Scaling each partial product to a largest entry of 1 could round an entry to 0 that the next factor would restore, so a projection with entries far apart in size, such as a scale by (1e160, 1e-160) seen through the same scale, was removed and the art drew flat. --- crates/geom/src/projective.rs | 89 +++++++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 26 deletions(-) diff --git a/crates/geom/src/projective.rs b/crates/geom/src/projective.rs index 4b17aac60..ddf836a23 100644 --- a/crates/geom/src/projective.rs +++ b/crates/geom/src/projective.rs @@ -100,11 +100,11 @@ 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`). When the product overflows, it is formed one factor at a - /// time, each partial product scaled to a largest entry of 1, and returned divided by the - /// largest magnitude in its `w` row (the same map, with a `w` row whose largest entry is 1). - /// `None` when `a` isn't invertible, or when the product overflows and the result at that - /// scale has an entry that is not finite or is not invertible. + /// maps from and to is moved by `a`). When the product overflows, every entry is formed with + /// its powers of two kept apart, and the product is returned divided by the largest magnitude + /// in its `w` row (the same map, with a `w` row whose largest entry is 1). `None` when `a` + /// isn't invertible, or when the product overflows and the result at that scale has an entry + /// that is not finite or is not invertible. pub fn conjugated(&self, a: Affine) -> Option { if a.determinant().abs() < 1e-15 { return None; @@ -114,12 +114,28 @@ impl Homography { if c.is_finite() { return Some(c); } - let c = to.then_after_normalized(&self.normalized()?)?.then_after_normalized(&from)?; - let k = c.0[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); + // Each entry of the product is a sum of nine terms to[r][i]·self[i][j]·from[j][k]. Every + // factor is split into a mantissa and a power of two, so each term is a mantissa product + // and an exponent sum, and the whole product is scaled by one power of two at the end. + let split = |h: &Self| h.0.map(|row| row.map(frexp)); + let (t, s, f) = (split(&to), split(self), split(&from)); + let term = |r: usize, i: usize, j: usize, k: usize| { + let (a, b, c) = (t[r][i], s[i][j], f[j][k]); + (a.0 * b.0 * c.0, a.1 + b.1 + c.1) + }; + let ijk = (0..27).map(|n| (n / 9, n / 3 % 3, n % 3)); + let e = ijk.map(|(i, j, k)| term(2, i, j, k)).filter(|(m, _)| *m != 0.0).map(|(_, e)| e).max()?; + let mut c = [[0.0; 3]; 3]; + for (r, row) in c.iter_mut().enumerate() { + for (k, v) in row.iter_mut().enumerate() { + *v = (0..9).map(|n| term(r, n / 3, n % 3, k)).map(|(m, x)| ldexp(m, x - e)).sum(); + } + } + let k = c[2].iter().fold(0.0_f64, |k, v| k.max(v.abs())); if k <= 0.0 { return None; } - let c = Self(c.0.map(|row| row.map(|v| v / k))); + let c = Self(c.map(|row| row.map(|v| v / k))); (c.is_finite() && c.det().abs() > 1e-15).then_some(c) } @@ -127,24 +143,6 @@ impl Homography { self.0.iter().flatten().all(|v| v.is_finite()) } - /// `self` divided by its largest entry magnitude, a positive scale that maps every point to - /// the same place; `None` when an entry is not finite or every entry is 0. - fn normalized(&self) -> Option { - let k = self.0.iter().flatten().fold(0.0_f64, |k, v| k.max(v.abs())); - (self.is_finite() && k > 0.0).then(|| Self(self.0.map(|row| row.map(|v| v / k)))) - } - - /// `self ∘ other` divided by its largest entry magnitude; when that product overflows, it is - /// formed from both factors divided by their largest entry magnitudes. `None` when the result - /// has an entry that is not finite or every entry is 0. - fn then_after_normalized(&self, other: &Self) -> Option { - let c = self.then_after(other); - if c.is_finite() { - return c.normalized(); - } - self.normalized()?.then_after(&other.normalized()?).normalized() - } - /// The affine map that matches `self` to first order at `p` (its derivative there), `None` /// beyond the horizon or when the derivative is not finite. pub fn affine_at(&self, p: Point) -> Option { @@ -171,6 +169,36 @@ impl Homography { } } +/// `v` as `m·2^e` with `0.5 ≤ |m| < 1`; 0, infinities and NaN come back as `(v, 0)`. +fn frexp(v: f64) -> (f64, i32) { + if v == 0.0 || !v.is_finite() { + return (v, 0); + } + // A subnormal is scaled into the normal range first. + let (v, bias) = if v.abs() < f64::MIN_POSITIVE { (v * pow2(64), -64) } else { (v, 0) }; + let bits = v.to_bits(); + let e = ((bits >> 52) & 0x7ff) as i32 - 1022; + (f64::from_bits((bits & !(0x7ff << 52)) | (1022 << 52)), e + bias) +} + +/// `m·2^e`, multiplied in steps of powers of two that are each within the normal range. +fn ldexp(mut m: f64, mut e: i32) -> f64 { + while e > 1023 && m.is_finite() { + m *= pow2(1023); + e -= 1023; + } + while e < -1022 && m != 0.0 { + m *= pow2(-1022); + e += 1022; + } + m * pow2(e.clamp(-1022, 1023)) +} + +/// `2^e` for `-1022 ≤ e ≤ 1023`. +fn pow2(e: i32) -> f64 { + f64::from_bits(((e + 1023) as u64) << 52) +} + #[cfg(test)] mod tests { use super::*; @@ -288,6 +316,15 @@ mod tests { assert!((v - w).abs() <= 1e-9 * w.abs().max(1.0), "{c:?}"); } assert_eq!(Homography::from_array(c.to_array()), Some(c)); + // A scale by (1e160, 1e-160) seen through the same scale is itself, though the product's + // first entry is 1e320 and a scale to a largest entry of 1 takes its second to 1e-480. + let thin = Homography([[1e160, 0.0, 0.0], [0.0, 1e-160, 0.0], [0.0, 0.0, 1.0]]); + assert!(Homography::from_array(thin.to_array()).is_some()); + let c = thin.conjugated(Affine::scale_non_uniform(1e160, 1e-160)).unwrap(); + for (v, w) in c.to_array().iter().zip(thin.to_array()) { + assert!((v - w).abs() <= 1e-9 * w.abs(), "{c:?}"); + } + assert_eq!(Homography::from_array(c.to_array()), Some(c)); // A move by 1e200 seen through a scale by 1e150 sends every point to about (1e350, 0), // which no invertible map of finite numbers draws. let far = Homography([[1.0, 0.0, 1e200], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]);