fix(geom): weight-aware error floor

crates/geom/src/nurbs.rs

@@ -55,13 +55,15 @@ const MAX_GAUSS_DEPTH: usize = 10;
 
 /// `gauss` over [a, b], halving until the two halves agree with the whole: a fixed panel count
 /// is far off on a tight bend, and the arc-length inverse must agree with the measure.
-fn adaptive_gauss(a: f64, b: f64, whole: f64, f: &impl Fn(f64) -> f64, depth: usize) -> f64 {
+/// `floor` is an absolute error nobody can resolve (rounding noise of the control coordinates): on
+/// a constant spline the speed is pure noise, which a relative test alone never accepts.
+fn adaptive_gauss(a: f64, b: f64, whole: f64, f: &impl Fn(f64) -> f64, floor: f64, depth: usize) -> f64 {
     let m = (a + b) / 2.0;
     let (l, r) = (gauss(a, m, f), gauss(m, b, f));
-    if depth == 0 || !(l + r).is_finite() || (l + r - whole).abs() <= 1e-10 * (l + r).abs() {
+    if depth == 0 || !(l + r).is_finite() || (l + r - whole).abs() <= 1e-10 * (l + r).abs() + floor {
         return l + r;
     }
-    adaptive_gauss(a, m, l, f, depth - 1) + adaptive_gauss(m, b, r, f, depth - 1)
+    adaptive_gauss(a, m, l, f, floor, depth - 1) + adaptive_gauss(m, b, r, f, floor, depth - 1)
 }
 
 /// Clamped knots over Bézier pieces meeting at `breaks`, for degree `p`.
@@ -527,7 +529,22 @@ impl Spline {
     /// sample.
     fn panel_length(&self, a: f64, b: f64, derivative: fn(&Spline, f64) -> Vec2) -> f64 {
         let speed = |t: f64| derivative(self, t).len();
-        adaptive_gauss(a, b, gauss(a, b, speed), &speed, MAX_GAUSS_DEPTH)
+        adaptive_gauss(a, b, gauss(a, b, speed), &speed, self.length_floor(a, b), MAX_GAUSS_DEPTH)
+    }
+
+    /// The absolute length error below which the panel [a, b] counts as converged: rounding of the
+    /// coordinates of the control points that shape it (1e-15 of the largest), so a distant control
+    /// point does not relax the tolerance. Costs one span's control points, not the whole spline.
+    fn length_floor(&self, a: f64, b: f64) -> f64 {
+        let Some(n) = self.control.len().checked_sub(1) else { return 0.0 };
+        let span = find_span(n, self.degree, (a + b) / 2.0, &self.knots);
+        let first = span.saturating_sub(self.degree);
+        let last = span.min(n);
+        let local = self.control.get(first..=last).unwrap_or_default();
+        // A control point's pull on the curve is its coordinate times its share of the weight, so a
+        // far point of tiny weight does not loosen the floor.
+        let heaviest = (first..=last).map(|i| self.weight(i)).fold(f64::MIN_POSITIVE, f64::max);
+        1e-15 * local.iter().enumerate().map(|(k, c)| c.x.abs().max(c.y.abs()) * (self.weight(first + k) / heaviest).min(1.0)).fold(0.0, f64::max)
     }
 
     /// Arc length between parameters `a` and `b` (adaptive Gauss–Legendre over each knot span).
@@ -618,10 +635,11 @@ impl Ellipse {
         let h = sweep / n as f64;
         // Each panel halved where it needs it: a flat ellipse turns sharply at its ends.
         let speed = |t: f64| self.derivative(t).len();
+        let floor = 0.0;
         let l: f64 = (0..n)
             .map(|i| {
                 let (a, b) = (t0 + h * i as f64, t0 + h * (i + 1) as f64);
-                adaptive_gauss(a, b, gauss(a, b, speed), &speed, MAX_GAUSS_DEPTH)
+                adaptive_gauss(a, b, gauss(a, b, speed), &speed, floor, MAX_GAUSS_DEPTH)
             })
             .sum();
         if l.is_finite() { l } else { 0.0 }
@@ -659,6 +677,40 @@ impl Ellipse {
 mod tests {
     use super::*;
 
+    /// A constant spline's speed is rounding noise: the integration accepts it at once rather than
+    /// halving to the depth cap (1,024 panels per knot span).
+    #[test]
+    fn noise_speed_converges_with_the_absolute_floor() {
+        let calls = std::cell::Cell::new(0usize);
+        let noise = |t: f64| {
+            calls.set(calls.get() + 1);
+            1e-16 * (t * 977.0).sin().abs()
+        };
+        let l = adaptive_gauss(0.0, 1.0, gauss(0.0, 1.0, noise), &noise, 1e-15, MAX_GAUSS_DEPTH);
+        assert!(l < 1e-15 && calls.get() < 100, "{l} {}", calls.get());
+        let c = Spline::from_control(vec![Vec2::new(1e6, 1e6); 5], 3);
+        assert!(c.length() < 1e-6);
+    }
+
+    #[test]
+    fn a_low_weight_control_point_does_not_relax_the_floor() {
+        let mut s = Spline::from_control(vec![Vec2::new(0.0, 0.0), Vec2::new(1e13, 0.0), Vec2::new(0.0, 0.0)], 2);
+        s.knots = vec![0.0, 0.0, 0.0, 1.0, 1.0, 1.0];
+        s.weights = vec![1.0, 1e-13, 0.25];
+        let floor = s.length_floor(0.0, 1.0);
+        assert!(floor <= 1.5e-15, "{floor}");
+    }
+
+    /// A distant control point does not relax the tolerance of a span it does not shape.
+    #[test]
+    fn the_error_floor_is_local_to_the_span() {
+        let mut s =
+            Spline::from_control(vec![Vec2::new(0.0, 0.0), Vec2::new(1.0, 0.0), Vec2::new(0.2, 0.0), Vec2::new(1e13, 0.0), Vec2::new(1e13, 0.0)], 2);
+        s.knots = vec![0.0, 0.0, 0.0, 1.0, 1.0, 2.0, 2.0, 2.0];
+        let got = s.length_between(0.0, 1.0);
+        assert!((got - 41.0 / 45.0).abs() < 1e-9, "{got}");
+    }
+
     fn cubic() -> Spline {
         Spline::from_control(vec![Vec2::new(0.0, 0.0), Vec2::new(3.0, 3.0), Vec2::new(6.0, -3.0), Vec2::new(9.0, 0.0), Vec2::new(12.0, 2.0)], 3)
     }

crates/geom/src/region.rs

@@ -58,11 +58,18 @@ struct Graph {
     src: Vec<(usize, usize)>,
     /// The input polyline of every edge (lower vertex index first).
     edge_src: HashMap<(usize, usize), usize>,
+    /// Neighbour angles still to be measured by face tracing: a hub with thousands of spokes is
+    /// scanned once per face through it, so the scans are bounded in all ([`MAX_TRACE_WORK`]).
+    trace_work: std::cell::Cell<u64>,
 }
 
-/// Records the crossing of segments `i` and `j`, if any, as split parameters on both.
-fn split_pair(segs: &[(Vec2, Vec2)], splits: &mut [Vec<f64>], i: usize, j: usize) {
-    let (Some(si), Some(sj)) = (segs.get(i), segs.get(j)) else { return };
+/// Neighbour-angle measurements one graph allows its face traces in all (a few seconds' work); a
+/// star of more than about 14,000 spokes needs more and finds no region.
+const MAX_TRACE_WORK: u64 = 200_000_000;
+
+/// Records the crossing of segments `i` and `j`, if any, as split parameters on both; whether there was one.
+fn split_pair(segs: &[(Vec2, Vec2)], splits: &mut [Vec<f64>], i: usize, j: usize) -> bool {
+    let (Some(si), Some(sj)) = (segs.get(i), segs.get(j)) else { return false };
     let li = Line::new(si.0, si.1);
     let lj = Line::new(sj.0, sj.1);
     if let Some(x) = line_line(&li, &lj) {
@@ -72,7 +79,9 @@ fn split_pair(segs: &[(Vec2, Vec2)], splits: &mut [Vec<f64>], i: usize, j: usize
         if let Some(v) = splits.get_mut(j) {
             v.push(lj.param_of(x).clamp(0.0, 1.0));
         }
+        return true;
     }
+    false
 }
 
 /// Upper bound on candidate-pair *visits*. Every visit counts, including repeats of a pair that shares several grid
@@ -82,6 +91,9 @@ fn split_pair(segs: &[(Vec2, Vec2)], splits: &mut [Vec<f64>], i: usize, j: usize
 const MAX_PAIR_VISITS: u64 = 60_000_000;
 /// Upper bound on distinct segment pairs tested for a crossing.
 const MAX_PAIRS_TESTED: usize = 20_000_000;
+/// Upper bound on crossings kept (two split parameters each, 16 bytes): one region operation holds
+/// at most about 128 MB of them, far more than any real hatch has.
+const MAX_CROSSINGS: usize = 8_000_000;
 
 /// Split parameters of every segment (always `0` and `1`, plus every crossing), found with a uniform grid.
 fn pair_splits(segs: &[(Vec2, Vec2)], bb: &Bounds2) -> Option<Vec<Vec<f64>>> {
@@ -132,6 +144,7 @@ fn pair_splits_limited(segs: &[(Vec2, Vec2)], bb: &Bounds2, max_visits: u64, max
     }
     let mut splits: Vec<Vec<f64>> = vec![vec![0.0, 1.0]; n];
     let mut tested: usize = 0;
+    let mut crossings: usize = 0;
     let mut visits: u64 = 0;
     for (c, list) in grid.iter().enumerate() {
         let (gx, gy) = (c / cells, c % cells);
@@ -147,9 +160,9 @@ fn pair_splits_limited(segs: &[(Vec2, Vec2)], bb: &Bounds2, max_visits: u64, max
                 if i == j || gx != ix.max(jx) || gy != iy.max(jy) {
                     continue;
                 }
-                split_pair(segs, &mut splits, i, j);
+                crossings += usize::from(split_pair(segs, &mut splits, i, j));
                 tested += 1;
-                if tested > max_pairs {
+                if tested > max_pairs || crossings > MAX_CROSSINGS {
                     return None;
                 }
             }
@@ -161,9 +174,9 @@ fn pair_splits_limited(segs: &[(Vec2, Vec2)], bb: &Bounds2, max_visits: u64, max
             if i == j || (j < i && is_long.get(j).copied().unwrap_or(false)) {
                 continue;
             }
-            split_pair(segs, &mut splits, i, j);
+            crossings += usize::from(split_pair(segs, &mut splits, i, j));
             tested += 1;
-            if tested > max_pairs {
+            if tested > max_pairs || crossings > MAX_CROSSINGS {
                 return None;
             }
         }
@@ -200,7 +213,13 @@ fn build_with(polys: &[Vec<Vec2>], tol: f64, find_splits: impl Fn(&[(Vec2, Vec2)
     let bb = Bounds2::from_points(segs.iter().flat_map(|(a, b)| [*a, *b]));
     let n = segs.len();
     let splits = find_splits(&segs, &bb)?;
-    let mut g = Graph { pts: Vec::new(), adj: Vec::new(), src: Vec::with_capacity(n), edge_src: HashMap::new() };
+    let mut g = Graph {
+        pts: Vec::new(),
+        adj: Vec::new(),
+        src: Vec::with_capacity(n),
+        edge_src: HashMap::new(),
+        trace_work: std::cell::Cell::new(MAX_TRACE_WORK),
+    };
     let mut index: HashMap<(i64, i64), usize> = HashMap::new();
     let mut vid = |p: Vec2, g: &mut Graph| -> usize {
         let k = key(p, bb.min, tol);
@@ -304,7 +323,13 @@ fn trace(g: &Graph, u0: usize, v0: usize, max_steps: usize) -> Option<Vec<usize>
         let pv = *g.pts.get(v)?;
         let back = (*g.pts.get(u)? - pv).angle();
         let mut best: Option<(f64, usize)> = None;
-        for &w in g.adj.get(v)? {
+        let ns = g.adj.get(v)?;
+        let Some(left) = g.trace_work.get().checked_sub(ns.len() as u64) else {
+            g.trace_work.set(0);
+            return None;
+        };
+        g.trace_work.set(left);
+        for &w in ns {
             if w == u && g.adj.get(v).is_some_and(|l| l.len() > 1) {
                 continue;
             }
@@ -355,7 +380,12 @@ fn outer_face(g: &Graph, p: Vec2) -> Option<Vec<usize>> {
     for (_, a, b) in hits.into_iter().take(64) {
         let (pa, pb) = (*g.pts.get(a)?, *g.pts.get(b)?);
         let (u, v) = if (pb - pa).cross(p - pa) > 0.0 { (a, b) } else { (b, a) };
-        if let Some(ids) = trace(g, u, v, max_steps)
+        let traced = trace(g, u, v, max_steps);
+        // A spent budget ends the search: a cheaper, more distant face must not stand in for the nearer one.
+        if traced.is_none() && g.trace_work.get() == 0 {
+            return None;
+        }
+        if let Some(ids) = traced
             && ids.len() >= 3
             && let Some(lp) = points(g, &ids)
             && shoelace(&lp) > 0.0
@@ -409,7 +439,7 @@ fn prune(g: &Graph) -> Graph {
             }
         }
     }
-    Graph { pts: g.pts.clone(), adj, src: Vec::new(), edge_src: HashMap::new() }
+    Graph { pts: g.pts.clone(), adj, src: Vec::new(), edge_src: HashMap::new(), trace_work: std::cell::Cell::new(MAX_TRACE_WORK) }
 }
 
 /// The outline of the outside of a group of curves on cycles, traced counter-clockwise from its
@@ -806,6 +836,30 @@ mod tests {
         assert_eq!(old.as_deref().map(sorted_bits), new.as_deref().map(sorted_bits), "n={} visits<={max_visits} pairs<={max_pairs}", segs.len());
     }
 
+    #[test]
+    fn crossings_kept_are_capped() {
+        // 2,900 horizontal x 2,900 vertical lines cross 8.41 million times.
+        let n = 2_900;
+        let mut segs: Vec<(Vec2, Vec2)> = (0..n).map(|i| (Vec2::new(0.0, i as f64 + 0.5), Vec2::new(n as f64, i as f64 + 0.5))).collect();
+        segs.extend((0..n).map(|i| (Vec2::new(i as f64 + 0.5, 0.0), Vec2::new(i as f64 + 0.5, n as f64))));
+        let bb = bounds_of(&segs);
+        assert!(pair_splits_limited(&segs, &bb, u64::MAX, usize::MAX).is_none());
+    }
+
+    #[test]
+    fn face_tracing_through_a_hub_is_budgeted() {
+        // 200 spokes from one hub to a ring: each face passes the hub and scans all its neighbours.
+        let k = 200;
+        let ring: Vec<Vec2> = (0..=k).map(|i| Vec2::new(10.0, 0.0).rotate(i as f64 * crate::TAU / k as f64)).collect();
+        let mut polys = vec![ring.clone()];
+        polys.extend(ring.iter().take(k).map(|q| vec![Vec2::ZERO, *q]));
+        let p = Vec2::new(5.0, 0.01);
+        let g = build(&polys, 1e-9).unwrap();
+        assert!(outer_face(&g, p).is_some());
+        g.trace_work.set(10);
+        assert!(outer_face(&g, p).is_none());
+    }
+
     #[test]
     fn pair_search_matches_oracle_with_budgets() {
         // Small inputs in few cells: many pairs share several cells. Sweeping each budget over