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 349ddb80f..ddf836a23 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,51 @@ 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`). 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; } - Some(Self::from_affine(a).then_after(self).then_after(&Self::from_affine(a.inverse()))) + 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); + } + // 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.map(|row| row.map(|v| v / k))); + (c.is_finite() && c.det().abs() > 1e-15).then_some(c) + } + + fn is_finite(&self) -> bool { + self.0.iter().flatten().all(|v| v.is_finite()) } /// 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 +152,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 @@ -127,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::*; @@ -175,6 +247,93 @@ 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_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.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 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 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]]); + 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] fn rect_bbox_is_the_corners_box() { let b = h().map_rect_bbox(Rect::new(0.0, 0.0, 10.0, 20.0)).unwrap();