render: the even-odd fill cuts a slab at the crossings inside it (#927) (#1308)

Slabs ran between vertex rows and paired the active edges by their
order at the slab's middle, so two edges crossing inside a slab (a
bowtie, two loops crossing) were paired across the crossing and the
trapezoids covered the wrong side of it. A slab with such a crossing is
now cut again at that y and each part paired on its own, under a cut
budget so hostile input stays bounded. Tests check coverage against an
even-odd parity oracle for a bowtie, an off-centre bowtie, two crossing
loops and the plain cases; the three crossing cases fail on main.


Claude-Session: https://claude.ai/code/session_011K1peyrh5MdZ6ttHkVq8Gh

Co-authored-by: Claude Fable 5.1 <noreply@anthropic.com>

crates/render/src/fill.rs

@@ -4,8 +4,10 @@ use cadcraft_doc::LwPolyline;
 use cadcraft_geom::{Polyline, Vec2};
 
 /// Triangulate an even-odd region bounded by closed loops (holes allowed) by cutting it into
-/// horizontal slabs at every vertex y and emitting trapezoids. Robust for any input; output
-/// size is O(slabs × crossings), capped for hostile inputs.
+/// horizontal slabs at every vertex y and emitting trapezoids. Two edges that cross inside a slab
+/// (a bowtie, two loops crossing each other) swap order there, so such a slab is cut again at the
+/// crossing and each part is paired on its own. Robust for any input; output size is
+/// O(slabs × crossings), capped for hostile inputs.
 pub fn triangulate_evenodd(loops: &[Vec<Vec2>]) -> Vec<Vec2> {
     let mut edges: Vec<(Vec2, Vec2)> = Vec::new();
     for l in loops {
@@ -30,6 +32,10 @@ pub fn triangulate_evenodd(loops: &[Vec<Vec2>]) -> Vec<Vec2> {
     let mut out = Vec::new();
     let mut active: Vec<(Vec2, Vec2)> = Vec::new();
     let mut next = 0usize;
+    // Slabs created by cutting at crossings, over the whole region; past this the remaining
+    // crossings are left unsplit rather than letting a hostile input run away.
+    let mut cuts_left = MAX_CROSSING_CUTS;
+    let mut pending: Vec<(f64, f64)> = Vec::new();
     for w in ys.windows(2) {
         let (Some(&y0), Some(&y1)) = (w.first(), w.get(1)) else { continue };
         let ym = (y0 + y1) / 2.0;
@@ -42,25 +48,74 @@ pub fn triangulate_evenodd(loops: &[Vec<Vec2>]) -> Vec<Vec2> {
             }
         }
         active.retain(|e| e.1.y > ym);
-        let xat = |e: &(Vec2, Vec2), y: f64| {
-            let t = (y - e.0.y) / (e.1.y - e.0.y);
-            e.0.x + (e.1.x - e.0.x) * t
-        };
-        let mut xs: Vec<(f64, f64, f64)> = active.iter().map(|e| (xat(e, ym), xat(e, y0), xat(e, y1))).collect();
-        xs.sort_by(|a, b| a.0.total_cmp(&b.0));
-        for pair in xs.chunks(2) {
-            if let [l, r] = pair {
-                let (a, b, c, d) = (Vec2::new(l.1, y0), Vec2::new(r.1, y0), Vec2::new(r.2, y1), Vec2::new(l.2, y1));
-                out.extend_from_slice(&[a, b, c, a, c, d]);
+        pending.clear();
+        pending.push((y0, y1));
+        while let Some((ya, yb)) = pending.pop() {
+            let ymid = (ya + yb) / 2.0;
+            let mut xs: Vec<(f64, f64, f64)> = active.iter().map(|e| (x_at(e, ymid), x_at(e, ya), x_at(e, yb))).collect();
+            xs.sort_by(|a, b| a.0.total_cmp(&b.0));
+            if cuts_left > 0 {
+                let mut cuts: Vec<f64> = Vec::new();
+                for pair in xs.windows(2) {
+                    let (Some(l), Some(r)) = (pair.first(), pair.get(1)) else { continue };
+                    if let Some(y) = crossing_y(l, r, ya, yb) {
+                        cuts.push(y);
+                    }
+                }
+                if !cuts.is_empty() {
+                    cuts.sort_by(f64::total_cmp);
+                    cuts.dedup_by(|a, b| (*a - *b).abs() < 1e-12);
+                    cuts.truncate(cuts_left);
+                    cuts_left -= cuts.len();
+                    // Push the parts top-down so the lowest is processed first.
+                    let mut top = yb;
+                    for &c in cuts.iter().rev() {
+                        pending.push((c, top));
+                        top = c;
+                    }
+                    pending.push((ya, top));
+                    continue;
+                }
+            }
+            for pair in xs.chunks(2) {
+                if let [l, r] = pair {
+                    let (a, b, c, d) = (Vec2::new(l.1, ya), Vec2::new(r.1, ya), Vec2::new(r.2, yb), Vec2::new(l.2, yb));
+                    out.extend_from_slice(&[a, b, c, a, c, d]);
+                }
+            }
+            if out.len() > 6_000_000 {
+                return out;
             }
-        }
-        if out.len() > 6_000_000 {
-            break;
         }
     }
     out
 }
 
+/// How many extra slab boundaries crossings may add over one region.
+const MAX_CROSSING_CUTS: usize = 100_000;
+
+fn x_at(e: &(Vec2, Vec2), y: f64) -> f64 {
+    let t = (y - e.0.y) / (e.1.y - e.0.y);
+    e.0.x + (e.1.x - e.0.x) * t
+}
+
+/// Where two edges ordered `l` before `r` at the slab's middle swap sides, if that happens
+/// strictly inside the slab. Each entry is (x at the middle, x at `ya`, x at `yb`).
+fn crossing_y(l: &(f64, f64, f64), r: &(f64, f64, f64), ya: f64, yb: f64) -> Option<f64> {
+    let (da, db) = (l.1 - r.1, l.2 - r.2);
+    // Ordered at the middle, so a swap shows as `l` right of `r` at one boundary.
+    if da <= 0.0 && db <= 0.0 {
+        return None;
+    }
+    // The gap is linear in y: it is zero at this fraction of the way from ya to yb.
+    let t = da / (da - db);
+    if !t.is_finite() {
+        return None;
+    }
+    let y = ya + (yb - ya) * t;
+    (y > ya && y < yb && y - ya > 1e-12 && yb - y > 1e-12).then_some(y)
+}
+
 /// Triangles for a polyline with widths (constant or tapered per segment).
 pub fn wide_polyline(p: &LwPolyline, tol: f64) -> Vec<Vec2> {
     let pl = Polyline { vertices: p.vertices.clone(), closed: p.closed };
@@ -115,6 +170,107 @@ mod tests {
         assert!((area(&t) - 6.0).abs() < 1e-9);
     }
 
+    /// How many of the triangles contain `p` (points on an edge are avoided by the callers).
+    fn cover(tris: &[Vec2], p: Vec2) -> usize {
+        tris.chunks(3)
+            .filter(|t| {
+                let s = |a: Vec2, b: Vec2| (b - a).cross(p - a);
+                let (u, v, w) = (s(t[0], t[1]), s(t[1], t[2]), s(t[2], t[0]));
+                (u >= 0.0 && v >= 0.0 && w >= 0.0) || (u <= 0.0 && v <= 0.0 && w <= 0.0)
+            })
+            .count()
+    }
+
+    /// Even-odd parity of `p` against the loops, as an oracle independent of the triangulation.
+    fn inside(loops: &[Vec<Vec2>], p: Vec2) -> bool {
+        let mut crossings = 0;
+        for l in loops {
+            for i in 0..l.len() {
+                let (a, b) = (l[i], l[(i + 1) % l.len()]);
+                if (a.y > p.y) != (b.y > p.y) {
+                    let x = a.x + (p.y - a.y) / (b.y - a.y) * (b.x - a.x);
+                    if x > p.x {
+                        crossings += 1;
+                    }
+                }
+            }
+        }
+        crossings % 2 == 1
+    }
+
+    /// Every sample point off the edges is covered exactly once when the oracle says inside and
+    /// never when it says outside.
+    fn assert_matches_oracle(loops: &[Vec<Vec2>], lo: f64, hi: f64) {
+        let t = triangulate_evenodd(loops);
+        let n = 23;
+        for i in 0..n {
+            for j in 0..n {
+                // Irrational-ish offsets keep the samples off every vertex, edge and crossing.
+                let p = Vec2::new(lo + (hi - lo) * (i as f64 + 0.3137) / n as f64, lo + (hi - lo) * (j as f64 + 0.6180) / n as f64);
+                let want = usize::from(inside(loops, p));
+                assert_eq!(cover(&t, p), want, "coverage at {p:?}");
+            }
+        }
+    }
+
+    #[test]
+    fn bowtie_crossing_inside_a_slab() {
+        // The crossing at (1, 1) lies strictly inside the only slab, y in [0, 2].
+        let bowtie = vec![Vec2::new(0.0, 0.0), Vec2::new(2.0, 2.0), Vec2::new(2.0, 0.0), Vec2::new(0.0, 2.0)];
+        let t = triangulate_evenodd(std::slice::from_ref(&bowtie));
+        assert!((area(&t) - 2.0).abs() < 1e-9, "two lobes of area 1, got {}", area(&t));
+        assert_matches_oracle(&[bowtie], -0.5, 2.5);
+    }
+
+    #[test]
+    fn bowtie_off_centre_crossing() {
+        // Crossing at y = 1.2, not at the slab's middle, so the middle ordering alone is wrong too.
+        let bowtie = vec![Vec2::new(0.0, 0.0), Vec2::new(3.0, 3.0), Vec2::new(3.0, 0.0), Vec2::new(0.0, 2.0)];
+        assert_matches_oracle(&[bowtie], -0.5, 3.5);
+    }
+
+    #[test]
+    fn two_loops_crossing_each_other() {
+        // Two squares overlapping: even-odd leaves the overlap empty; their crossings sit inside
+        // the slabs between the vertex rows.
+        let a = vec![Vec2::new(0.0, 0.0), Vec2::new(4.0, 1.0), Vec2::new(3.0, 5.0), Vec2::new(-1.0, 4.0)];
+        let b = vec![Vec2::new(2.0, 2.0), Vec2::new(6.0, 3.0), Vec2::new(5.0, 7.0), Vec2::new(1.0, 6.0)];
+        let t = triangulate_evenodd(&[a.clone(), b.clone()]);
+        // Each parallelogram has area 17; even-odd drops the overlap twice. The overlap is
+        // sampled on a grid, so compare loosely; the oracle below is the exact check.
+        assert!((area(&t) - (2.0 * 17.0 - 2.0 * overlap_area(&a, &b))).abs() < 0.3, "area {}", area(&t));
+        assert_matches_oracle(&[a, b], -1.5, 7.5);
+    }
+
+    /// Area where both loops are inside, sampled on a fine grid; good to a few hundredths.
+    fn overlap_area(a: &[Vec2], b: &[Vec2]) -> f64 {
+        let (la, lb) = (vec![a.to_vec()], vec![b.to_vec()]);
+        let n = 400;
+        let (lo, hi) = (-2.0, 8.0);
+        let cell = (hi - lo) / n as f64;
+        let mut both = 0usize;
+        for i in 0..n {
+            for j in 0..n {
+                let p = Vec2::new(lo + (i as f64 + 0.5) * cell, lo + (j as f64 + 0.5) * cell);
+                if inside(&la, p) && inside(&lb, p) {
+                    both += 1;
+                }
+            }
+        }
+        both as f64 * cell * cell
+    }
+
+    #[test]
+    fn simple_loops_still_match_the_oracle() {
+        let outer = vec![Vec2::new(0.0, 0.0), Vec2::new(10.0, 0.0), Vec2::new(10.0, 10.0), Vec2::new(0.0, 10.0)];
+        let hole = vec![Vec2::new(4.0, 4.0), Vec2::new(6.0, 4.0), Vec2::new(6.0, 6.0), Vec2::new(4.0, 6.0)];
+        assert_matches_oracle(&[outer.clone(), hole.clone()], -1.0, 11.0);
+        let reversed: Vec<Vec2> = outer.iter().rev().copied().collect();
+        assert_matches_oracle(&[reversed, hole], -1.0, 11.0);
+        let tri = vec![Vec2::new(0.0, 0.0), Vec2::new(4.0, 0.0), Vec2::new(0.0, 3.0)];
+        assert_matches_oracle(&[tri], -1.0, 5.0);
+    }
+
     #[test]
     fn degenerate_input() {
         assert!(triangulate_evenodd(&[vec![Vec2::ZERO, Vec2::X]]).is_empty());