Merge pull request #3717 from storytold/land/3357

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

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)
     }