G11 #937: Shape Builder cuts art of thousands of segments with lines again

The planar map found each edge's source input by scanning every input
segment, and the face around each hole or merge point by scanning every
face, so it was quadratic in the segments; the Shape Builder therefore
gave up past 4096 segments (traced or type outlines) and offered whole
shapes. Index outline segments and faces by bounds on both axes
(RectIndex), walk an edge's segments lazily when finding its leaving
angle, and raise SHAPE_BUILDER_MAX_SEGMENTS to 65,536. Divide with open
paths shares the cap; its cap test now follows the constant.

Tests: a 20,000-segment blob cut by two lines, a 300-line hatch, the
index against a brute-force scan, and the existing cap test at the new
cap. Progress docs updated (G11, parity evidence, roadmap).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014gCpjL3yUhaH3jRp7ksufh

crates/pathops/src/pathfinder.rs

@@ -937,9 +937,11 @@ mod tests {
     #[test]
     fn divide_open_paths_retains_the_existing_segment_cap() {
         let mut inputs = vec![Shape::new(rect(0.0, 0.0, 100.0, 100.0), FillRule::NonZero, 7)];
-        inputs.extend((0..4093).map(|_| Shape::new(shapes::line(Point::new(-10.0, 50.0), Point::new(110.0, 50.0)), FillRule::NonZero, 99)));
+        // The square's four sides and one line past the cap.
+        let n = crate::SHAPE_BUILDER_MAX_SEGMENTS - 3;
+        inputs.extend((0..n).map(|_| Shape::new(shapes::line(Point::new(-10.0, 50.0), Point::new(110.0, 50.0)), FillRule::NonZero, 99)));
         let faces = pathfinder(PathfinderOp::Divide, &inputs);
-        assert_eq!(faces.len(), 1, "4097 segments fall back to filled areas");
+        assert_eq!(faces.len(), 1, "a segment past the cap falls back to filled areas");
         assert_eq!(faces[0].key, 7);
         assert!((crate::area(&faces[0].path, faces[0].rule) - 10000.0).abs() < 1e-6);
     }

crates/pathops/src/planar.rs

@@ -94,9 +94,65 @@ fn mid(segs: &[Seg]) -> Option<Point> {
     segs.get(segs.len() / 2).map(|s| s.c.eval(0.5))
 }
 
+/// Rectangles looked up by a point inside them, when there are too many to test one by one for
+/// every point (the map asks once per edge and once per hole, so a scan of them all made it
+/// quadratic in the input, #937). Rectangles narrow in x are sorted by their left edge, so the
+/// ones that can hold a point are one short run of them; those wide in x but narrow in y (a hatch
+/// of long horizontal lines) are sorted by their top edge the same way; the few wide both ways (a
+/// long diagonal across the art) are scanned. "Narrow" is up to twice the side of a cell holding
+/// one rectangle on average, so a run holds about the rectangles of one row or column of cells.
+struct RectIndex<T> {
+    by_x: Vec<(kurbo::Rect, T)>,
+    max_width: f64,
+    by_y: Vec<(kurbo::Rect, T)>,
+    max_height: f64,
+    wide: Vec<(kurbo::Rect, T)>,
+}
+
+impl<T> RectIndex<T> {
+    fn new(items: Vec<(kurbo::Rect, T)>) -> Self {
+        let n = items.len().max(1) as f64;
+        let all = items.iter().map(|(r, _)| *r).reduce(|a, b| a.union(b));
+        let cell = all
+            .map(|a| {
+                let side = (a.area() / n).sqrt();
+                if side > 0.0 { 2.0 * side } else { 2.0 * a.width().max(a.height()) / n }
+            })
+            .filter(|c| c.is_finite())
+            .unwrap_or(0.0)
+            .max(1e-9);
+        let (mut by_x, mut by_y, mut wide) = (Vec::new(), Vec::new(), Vec::new());
+        for it in items {
+            if it.0.width() <= cell {
+                by_x.push(it);
+            } else if it.0.height() <= cell {
+                by_y.push(it);
+            } else {
+                // A NaN extent fails both tests above, so a non-finite rectangle is never in a sorted
+                // run along that axis; `contains` rejects it wherever it lands.
+                wide.push(it);
+            }
+        }
+        by_x.sort_by(|a, b| a.0.x0.total_cmp(&b.0.x0));
+        by_y.sort_by(|a, b| a.0.y0.total_cmp(&b.0.y0));
+        let max_width = by_x.iter().map(|(r, _)| r.width()).fold(0.0, f64::max);
+        let max_height = by_y.iter().map(|(r, _)| r.height()).fold(0.0, f64::max);
+        Self { by_x, max_width, by_y, max_height, wide }
+    }
+
+    /// The items whose rectangle, grown by `tol`, holds `p`.
+    fn at(&self, p: Point, tol: f64) -> impl Iterator<Item = &T> {
+        let lo = self.by_x.partition_point(|(r, _)| r.x0 < p.x - tol - self.max_width);
+        let run_x = self.by_x.get(lo..).unwrap_or_default().iter().take_while(move |(r, _)| r.x0 <= p.x + tol);
+        let lo = self.by_y.partition_point(|(r, _)| r.y0 < p.y - tol - self.max_height);
+        let run_y = self.by_y.get(lo..).unwrap_or_default().iter().take_while(move |(r, _)| r.y0 <= p.y + tol);
+        run_x.chain(run_y).chain(&self.wide).filter(move |(r, _)| r.inflate(tol, tol).contains(p)).map(|(_, t)| t)
+    }
+}
+
 /// The inputs' segments (in the sweep's frame), to tell which inputs an edge lies on.
 struct Outlines {
-    segs: Vec<(kurbo::Rect, PathSeg, usize)>,
+    segs: RectIndex<(PathSeg, usize)>,
     tol: f64,
 }
 
@@ -105,18 +161,18 @@ impl Outlines {
         // A zero-length segment (a point) is the source of no edge, even one running through it.
         let segs = inputs
             .iter()
-            .flat_map(|(bp, i)| bp.segments().map(move |s| (s.bounding_box(), s, *i)))
+            .flat_map(|(bp, i)| bp.segments().map(move |s| (s.bounding_box(), (s, *i))))
             .filter(|(r, ..)| r.width() > 0.0 || r.height() > 0.0)
             .collect();
         // The sweep moves points by a few `eps` at most.
-        Self { segs, tol: 16.0 * eps }
+        Self { segs: RectIndex::new(segs), tol: 16.0 * eps }
     }
 
     /// The front-most input whose outline passes through `p`.
     fn source_at(&self, p: Point) -> Option<usize> {
         use kurbo::ParamCurveNearest;
         let tol = self.tol;
-        self.segs.iter().filter(|(r, s, _)| r.inflate(tol, tol).contains(p) && s.nearest(p, 1e-9).distance_sq <= tol * tol).map(|(_, _, i)| *i).max()
+        self.segs.at(p, tol).filter(|(s, _)| s.nearest(p, 1e-9).distance_sq <= tol * tol).map(|(_, i)| *i).max()
     }
 }
 
@@ -170,12 +226,24 @@ fn twin(h: usize) -> usize {
 
 /// The segments of half-edge `h`, in its direction.
 fn half_segs(edges: &[Edge], h: usize) -> Vec<Seg> {
-    let Some(e) = edges.get(h / 2) else { return Vec::new() };
-    if h.is_multiple_of(2) {
-        e.segs.clone()
-    } else {
-        e.segs.iter().rev().map(|s| Seg { c: CubicBez::new(s.c.p3, s.c.p2, s.c.p1, s.c.p0), line: s.line }).collect()
-    }
+    half_seg_iter(edges, h).collect()
+}
+
+/// The segments of half-edge `h` in its direction, without copying the edge first (an edge can
+/// hold thousands of segments, and the angle it leaves a vertex at needs the first few).
+fn half_seg_iter(edges: &[Edge], h: usize) -> impl Iterator<Item = Seg> + '_ {
+    let segs = edges.get(h / 2).map(|e| e.segs.as_slice()).unwrap_or_default();
+    let forward = h.is_multiple_of(2);
+    (0..segs.len()).filter_map(move |k| {
+        let s = segs.get(if forward { k } else { segs.len() - 1 - k })?;
+        Some(if forward { *s } else { Seg { c: CubicBez::new(s.c.p3, s.c.p2, s.c.p1, s.c.p0), line: s.line } })
+    })
+}
+
+/// Where half-edge `h` ends (the vertex at the other end of its edge).
+fn far_end(edges: &[Edge], h: usize) -> Option<Point> {
+    let e = edges.get(h / 2)?;
+    if h.is_multiple_of(2) { e.segs.last().map(|s| s.c.p3) } else { e.segs.first().map(|s| s.c.p0) }
 }
 
 /// The vertex half-edge `h` leaves from.
@@ -188,32 +256,33 @@ fn origin(edges: &[Edge], h: usize) -> Option<usize> {
 /// Edges that only meet at `v` cross a circle around it in the order they leave `v` in, however
 /// they start out (tangent curves, or the sweep's short straight stubs at a vertex).
 fn leaving_angle(edges: &[Edge], h: usize, v: Point, r: f64) -> f64 {
-    let segs = half_segs(edges, h);
-    let mut bp = BezPath::new();
-    for (k, s) in segs.iter().enumerate() {
-        if k == 0 {
-            bp.move_to(s.c.p0);
-        }
-        bp.curve_to(s.c.p1, s.c.p2, s.c.p3);
-    }
-    let mut prev = segs.first().map_or(v, |s| s.c.p0);
+    let first = half_seg_iter(edges, h).next();
+    let mut prev = first.map_or(v, |s| s.c.p0);
     let mut hit: Option<Point> = None;
-    kurbo::flatten(bp.iter(), (r * 1e-4).max(1e-12), |el| {
-        if let (None, kurbo::PathEl::LineTo(p)) = (hit, el) {
-            if p.distance(v) >= r {
-                // The point on prev → p at distance r (prev is closer).
-                let (d, f) = (p - prev, prev - v);
-                let (a, b, c) = (d.hypot2(), 2.0 * f.dot(d), f.hypot2() - r * r);
-                let t = if a > 0.0 { ((-b + (b * b - 4.0 * a * c).max(0.0).sqrt()) / (2.0 * a)).clamp(0.0, 1.0) } else { 1.0 };
-                hit = Some(prev + d * t);
+    // One segment at a time, stopping at the first that gets r away: an edge can run on for
+    // thousands of segments (a traced outline between two crossings), all of which flattening
+    // the whole edge to this tolerance would visit, at every vertex.
+    for s in half_seg_iter(edges, h) {
+        let path = [kurbo::PathEl::MoveTo(s.c.p0), kurbo::PathEl::CurveTo(s.c.p1, s.c.p2, s.c.p3)];
+        kurbo::flatten(path, (r * 1e-4).max(1e-12), |el| {
+            if let (None, kurbo::PathEl::LineTo(p)) = (hit, el) {
+                if p.distance(v) >= r {
+                    // The point on prev → p at distance r (prev is closer).
+                    let (d, f) = (p - prev, prev - v);
+                    let (a, b, c) = (d.hypot2(), 2.0 * f.dot(d), f.hypot2() - r * r);
+                    let t = if a > 0.0 { ((-b + (b * b - 4.0 * a * c).max(0.0).sqrt()) / (2.0 * a)).clamp(0.0, 1.0) } else { 1.0 };
+                    hit = Some(prev + d * t);
+                }
+                prev = p;
             }
-            prev = p;
+        });
+        if hit.is_some() {
+            break;
         }
-    });
+    }
     // An edge too short to reach r (it can't, r is at most half of every chord) leaves along its
     // tangent.
-    let far =
-        hit.unwrap_or_else(|| segs.first().map_or(v, |s| [s.c.p1, s.c.p2, s.c.p3].into_iter().find(|p| p.distance(s.c.p0) > 0.0).unwrap_or(s.c.p3)));
+    let far = hit.unwrap_or_else(|| first.map_or(v, |s| [s.c.p1, s.c.p2, s.c.p3].into_iter().find(|p| p.distance(s.c.p0) > 0.0).unwrap_or(s.c.p3)));
     (far - v).atan2()
 }
 
@@ -581,9 +650,11 @@ pub fn find_gaps(shapes: &[Shape], max_gap: f64) -> Vec<[Point; 2]> {
 }
 
 /// Most segments the Shape Builder sweeps into its planar map (open paths as edges, the gaps closed
-/// paths enclose, outline edges): the map costs far more than the filled areas' arrangement, so past
-/// this only the filled areas count.
-pub const SHAPE_BUILDER_MAX_SEGMENTS: usize = 4096;
+/// paths enclose, outline edges); past this only the filled areas count. The map costs about what
+/// the filled areas' arrangement does (both near linear in the input; the cap was 4096 while
+/// finding each edge's source scanned every segment, #937), so this only bounds the time a hover
+/// can take on art far past what a drawing holds.
+pub const SHAPE_BUILDER_MAX_SEGMENTS: usize = 65_536;
 
 /// The Shape Builder's arrangement: see [`shape_builder`].
 #[derive(Clone, Debug, Default, PartialEq)]
@@ -864,8 +935,10 @@ pub(crate) struct Map {
     edges: Vec<Edge>,
     around: Vec<Vec<usize>>,
     slot_of: Vec<usize>,
-    /// Bounded faces: their cycle of half-edges, path and area.
+    /// Bounded faces: their cycle of half-edges, path and area; and, by their bounds, which can
+    /// hold a point, so finding the face around one tests only those.
     faces: Vec<(Vec<usize>, BezPath, f64)>,
+    face_index: RectIndex<usize>,
     /// Outer boundaries of connected pieces, and the face each one is a hole in.
     holes: Vec<(Vec<usize>, BezPath)>,
     face_holes: Vec<Vec<usize>>,
@@ -1020,7 +1093,7 @@ pub(crate) fn map(shapes: &[Shape]) -> Option<Map> {
     for (v, list) in around.iter_mut().enumerate() {
         let Some(&at) = points.get(v) else { continue };
         // Half the shortest chord of the edges here: every one of them gets that far from it.
-        let r = 0.5 * list.iter().filter_map(|&h| half_segs(&edges, h).last().map(|s| s.c.p3.distance(at))).fold(f64::INFINITY, f64::min);
+        let r = 0.5 * list.iter().filter_map(|&h| far_end(&edges, h).map(|p| p.distance(at))).fold(f64::INFINITY, f64::min);
         let mut keyed: Vec<(f64, usize)> = list.iter().map(|&h| (leaving_angle(&edges, h, at, r), h)).collect();
         keyed.sort_by(|a, b| a.0.total_cmp(&b.0));
         *list = keyed.into_iter().map(|(_, h)| h).collect();
@@ -1081,14 +1154,15 @@ pub(crate) fn map(shapes: &[Shape]) -> Option<Map> {
     }
 
     // Each hole goes in the smallest face of another piece around it.
+    let face_index = RectIndex::new(faces.iter().enumerate().map(|(i, (_, bp, _))| (bp.bounding_box(), i)).collect());
     let mut face_holes: Vec<Vec<usize>> = vec![Vec::new(); faces.len()];
     for (k, (cycle, _)) in holes.iter().enumerate() {
         let Some(p) = cycle.first().and_then(|&h| origin(&edges, h)).and_then(|v| points.get(v)).copied() else { continue };
         let own = piece_of(cycle);
-        let around_it = faces
-            .iter()
-            .enumerate()
-            .filter(|(_, (fc, bp, _))| piece_of(fc) != own && bp.bounding_box().contains(p) && bp.winding(p) != 0)
+        let around_it = face_index
+            .at(p, 0.0)
+            .filter_map(|&i| faces.get(i).map(|f| (i, f)))
+            .filter(|(_, (fc, bp, _))| piece_of(fc) != own && bp.winding(p) != 0)
             .min_by(|a, b| a.1.2.total_cmp(&b.1.2))
             .map(|(i, _)| i);
         if let Some(list) = around_it.and_then(|i| face_holes.get_mut(i)) {
@@ -1100,7 +1174,7 @@ pub(crate) fn map(shapes: &[Shape]) -> Option<Map> {
     let snap = 2.0 * eps;
     let mut anchors: Vec<Point> = shapes.iter().flat_map(|s| s.path.subpaths.iter().flat_map(|sp| sp.anchors.iter().map(|a| a.p))).collect();
     anchors.sort_by(|a, b| a.x.total_cmp(&b.x).then(a.y.total_cmp(&b.y)));
-    Some(Map { edges, around, slot_of, faces, holes, face_holes, tidy, anchors, snap })
+    Some(Map { edges, around, slot_of, faces, face_index, holes, face_holes, tidy, anchors, snap })
 }
 
 impl Map {
@@ -1123,10 +1197,10 @@ impl Map {
             let q = SHEAR * p;
             // The smallest face around it: faces of a piece inside a hole lie within the face around.
             let (i, _) = self
-                .faces
-                .iter()
-                .enumerate()
-                .filter(|(_, (_, bp, _))| bp.bounding_box().contains(q) && bp.winding(q) != 0)
+                .face_index
+                .at(q, 0.0)
+                .filter_map(|&i| self.faces.get(i).map(|f| (i, f)))
+                .filter(|(_, (_, bp, _))| bp.winding(q) != 0)
                 .min_by(|a, b| a.1.2.total_cmp(&b.1.2))?;
             *picked.get_mut(i)? = true;
         }
@@ -1796,6 +1870,76 @@ mod tests {
         assert!(shape_builder(&s, false).lines.is_empty());
     }
 
+    #[test]
+    fn a_traced_outline_of_many_segments_is_cut_by_lines() {
+        // #937: art traced to thousands of segments was past the old cap (4096), so lines stopped
+        // cutting it and the Shape Builder offered the whole shapes. A blob of 20 000 segments
+        // crossed by two long lines (far wider than any of its segments: the index of the
+        // outlines keeps those apart) splits into four regions, each line in four pieces.
+        let n = 20_000;
+        let blob: Vec<Point> = (0..n)
+            .map(|k| {
+                let t = k as f64 / n as f64 * std::f64::consts::TAU;
+                let r = 100.0 + 5.0 * (t * 37.0).sin();
+                Point::new(r * t.cos(), r * t.sin())
+            })
+            .collect();
+        let s = [
+            Shape::new(PathData::single(SubPath::polyline(&blob, true)), FillRule::NonZero, 0),
+            line((-150.0, 0.0), (150.0, 3.0), 1),
+            line((0.0, -150.0), (3.0, 150.0), 2),
+        ];
+        let BuilderArrangement { regions, lines, .. } = shape_builder(&s, false);
+        assert_eq!(regions.len(), 4, "{}", regions.len());
+        assert!(regions.iter().all(|r| r.sources == [0]));
+        let whole = crate::area(&s[0].path, FillRule::NonZero);
+        assert!((regions.iter().map(area).sum::<f64>() - whole).abs() < 1e-6 * whole);
+        assert_eq!(lines.iter().filter(|l| l.key == 1).count(), 4, "cut twice by the blob and once by the other line");
+        assert_eq!(lines.iter().filter(|l| l.key == 2).count(), 4);
+    }
+
+    #[test]
+    fn a_hatch_of_parallel_lines_is_cut_line_by_line() {
+        // Hundreds of long horizontal lines over a square, all alike in x (the weak axis of an index
+        // by x alone), and one vertical line: every strip is cut in two.
+        let n: u32 = 300;
+        let mut s = vec![rect(0.0, 0.0, 100.0, 100.0, 0)];
+        for k in 1..n {
+            let y = 100.0 * f64::from(k) / f64::from(n);
+            s.push(line((-10.0, y), (110.0, y), u64::from(k)));
+        }
+        s.push(line((50.0, -10.0), (50.0, 110.0), u64::from(n)));
+        let BuilderArrangement { regions, lines, .. } = shape_builder(&s, false);
+        assert_eq!(regions.len(), 2 * n as usize, "{}", regions.len());
+        assert!(regions.iter().all(|r| r.sources == [0]));
+        assert!((regions.iter().map(area).sum::<f64>() - 10000.0).abs() < 1e-6);
+        // Each horizontal line: a loose end outside, two pieces inside, a loose end outside.
+        assert_eq!(lines.iter().filter(|l| l.key == 1).count(), 4);
+        // The vertical line: cut by the square's sides and every hatch line.
+        assert_eq!(lines.iter().filter(|l| l.key == u64::from(n)).count(), n as usize + 2);
+    }
+
+    #[test]
+    fn rect_index_finds_every_rectangle_holding_a_point() {
+        // Narrow in x, narrow in y only, and wide both ways: each lookup matches a scan of them all.
+        let mut items = vec![(kurbo::Rect::new(0.0, 0.0, 100.0, 100.0), 0usize)];
+        for k in 1..50 {
+            let a = f64::from(k as u32) * 2.0;
+            items.push((kurbo::Rect::new(a, 0.0, a + 1.0, 100.0), k));
+            items.push((kurbo::Rect::new(0.0, a, 100.0, a + 1.0), 100 + k));
+        }
+        items.push((kurbo::Rect::new(f64::NAN, 0.0, 1.0, 1.0), 999));
+        let index = RectIndex::new(items.clone());
+        for p in [Point::new(4.5, 6.5), Point::new(50.0, 50.0), Point::new(-1.0, 3.0), Point::new(99.5, 0.5)] {
+            let mut got: Vec<usize> = index.at(p, 0.1).copied().collect();
+            got.sort_unstable();
+            let mut want: Vec<usize> = items.iter().filter(|(r, _)| r.inflate(0.1, 0.1).contains(p)).map(|(_, k)| *k).collect();
+            want.sort_unstable();
+            assert_eq!(got, want, "{p:?}");
+        }
+        assert_eq!(RectIndex::<usize>::new(Vec::new()).at(Point::ZERO, 1.0).count(), 0);
+    }
+
     fn poly(pts: &[(f64, f64)], closed: bool) -> SubPath {
         SubPath::polyline(&pts.iter().map(|&(x, y)| Point::new(x, y)).collect::<Vec<_>>(), closed)
     }