Merge pull request #51 from storytold/pr/pathops-near-horizontal-edges

Pathops: Booleans keep their corners on almost-horizontal edges; deterministic property tests (M6.10)

crates/pathops/src/boolean.rs

@@ -146,6 +146,79 @@ fn el_points(el: &kurbo::PathEl) -> impl Iterator<Item = Point> {
     v.into_iter().flatten()
 }
 
+/// Make edges that rise by less than the sweep tolerance `eps` exactly horizontal.
+///
+/// The sweep orders a curve against another segment over a y-span shorter than `eps` by bounding
+/// boxes alone, so an almost-horizontal edge (typically a boolean result's edge whose end, an
+/// intersection point, is a rounding error off the input's horizontal) can be mis-ordered against
+/// a curve it passes and the output loses a corner; exact horizontals are handled exactly. So
+/// on-curve y coordinates within `eps` of each other snap to one value (the lowest of each
+/// cluster, keeping distinct values at least `eps` apart), handles move with their anchors, and a
+/// handle within `eps` of its anchor's height is levelled with it (no y-extremum is left in a
+/// sub-`eps` sliver of a curve). Anchors move by less than `eps` and handles by less than `2·eps`,
+/// the order of the sweep's own tolerance.
+fn snap_horizontals(paths: &mut [&mut BezPath], eps: f64) {
+    use kurbo::PathEl::*;
+    let mut ys: Vec<f64> = paths.iter().flat_map(|p| p.elements().iter().filter_map(kurbo::PathEl::end_point)).map(|p| p.y).collect();
+    ys.sort_by(f64::total_cmp);
+    ys.dedup();
+    let mut start = f64::NEG_INFINITY;
+    let levels: Vec<f64> = ys
+        .iter()
+        .map(|&y| {
+            if y - start >= eps {
+                start = y;
+            }
+            start
+        })
+        .collect();
+    let snap = |p: &mut Point| {
+        let y = levels[ys.partition_point(|&v| v < p.y)];
+        let moved = y - p.y;
+        p.y = y;
+        moved
+    };
+    // Shift a handle with its anchor, then level it with the anchor when within `eps`.
+    let follow = |h: &mut Point, anchor: Point, moved: f64| {
+        h.y += moved;
+        if (h.y - anchor.y).abs() < eps {
+            h.y = anchor.y;
+        }
+    };
+    for p in paths.iter_mut() {
+        // The current point (snapped) and how far snapping moved it; the subpath's start.
+        let (mut cur, mut moved, mut first) = (Point::ZERO, 0.0, Point::ZERO);
+        for el in p.elements_mut() {
+            match el {
+                MoveTo(q) => {
+                    moved = snap(q);
+                    (cur, first) = (*q, *q);
+                }
+                LineTo(q) => {
+                    moved = snap(q);
+                    cur = *q;
+                }
+                QuadTo(c, q) => {
+                    let m = snap(q);
+                    follow(c, cur, moved);
+                    follow(c, *q, 0.0);
+                    (cur, moved) = (*q, m);
+                }
+                CurveTo(c1, c2, q) => {
+                    let m = snap(q);
+                    follow(c1, cur, moved);
+                    follow(c2, *q, m);
+                    (cur, moved) = (*q, m);
+                }
+                ClosePath => {
+                    // A later segment without a MoveTo starts from the subpath's (snapped) start.
+                    (cur, moved) = (first, 0.0);
+                }
+            }
+        }
+    }
+}
+
 /// An arrangement (planar subdivision) of N filled paths.
 pub(crate) struct Arrangement {
     pub(crate) top: Topology<Multi>,
@@ -194,8 +267,9 @@ impl Tidy {
 }
 
 impl Arrangement {
-    pub(crate) fn new(paths: &[(BezPath, FillRule)]) -> Result<Self, PathOpsError> {
+    pub(crate) fn new(mut paths: Vec<(BezPath, FillRule)>) -> Result<Self, PathOpsError> {
         let eps = eps_for(paths.iter().map(|p| &p.0))?;
+        snap_horizontals(&mut paths.iter_mut().map(|p| &mut p.0).collect::<Vec<_>>(), eps);
         // The sweep can panic on rare near-degenerate input; that must never take the app down,
         // so it becomes an error (callers fall back to leaving the art unchanged).
         let sweep = std::panic::catch_unwind(std::panic::AssertUnwindSafe(|| {
@@ -456,7 +530,9 @@ pub fn try_normalize(path: &PathData, rule: FillRule) -> Result<PathData, PathOp
 
 pub(crate) fn normalize_bez(bp: &BezPath, rule: FillRule) -> Result<Contours, PathOpsError> {
     let eps = eps_for([bp])?;
-    let top = Topology::<i32>::from_path(bp, eps).map_err(|_| PathOpsError::OpenPath)?;
+    let mut bp = bp.clone();
+    snap_horizontals(&mut [&mut bp], eps);
+    let top = Topology::<i32>::from_path(&bp, eps).map_err(|_| PathOpsError::OpenPath)?;
     Ok(top.contours(|w| inside(rule, *w)))
 }
 
@@ -467,7 +543,7 @@ pub fn boolean(a: &PathData, a_rule: FillRule, b: &PathData, b_rule: FillRule, o
 
 /// Boolean operation with an explicit refit `precision`; errors on non-finite input.
 pub fn try_boolean(a: &PathData, a_rule: FillRule, b: &PathData, b_rule: FillRule, op: BoolOp, precision: f64) -> Result<PathData, PathOpsError> {
-    let arr = Arrangement::new(&[(fill_bezpath(a), a_rule), (fill_bezpath(b), b_rule)])?;
+    let arr = Arrangement::new(vec![(fill_bezpath(a), a_rule), (fill_bezpath(b), b_rule)])?;
     let c = arr.contours(|m| op.apply(m[0], m[1]));
     let tidy = Tidy { precision, ..arr.tidy };
     Ok(all_contours_to_path(&c, &tidy))
@@ -477,7 +553,7 @@ pub fn try_boolean(a: &PathData, a_rule: FillRule, b: &PathData, b_rule: FillRul
 /// point is inside `paths[i]` under its fill rule.
 pub fn boolean_n(paths: &[(&PathData, FillRule)], pred: impl Fn(&[bool]) -> bool) -> PathData {
     let bps: Vec<(BezPath, FillRule)> = paths.iter().map(|(p, r)| (fill_bezpath(p), *r)).collect();
-    match Arrangement::new(&bps) {
+    match Arrangement::new(bps) {
         Ok(arr) => all_contours_to_path(&arr.contours(pred), &arr.tidy),
         Err(_) => PathData::default(),
     }

crates/pathops/src/offset.rs

@@ -90,7 +90,7 @@ pub fn offset_path(path: &PathData, delta: f64, join: Join, miter_limit: f64) ->
     } else {
         BezPath::new()
     };
-    let Ok(arr) = Arrangement::new(&[(fill, FillRule::NonZero), (ring, FillRule::NonZero), (open_bp, FillRule::NonZero)]) else {
+    let Ok(arr) = Arrangement::new(vec![(fill, FillRule::NonZero), (ring, FillRule::NonZero), (open_bp, FillRule::NonZero)]) else {
         return PathData::default();
     };
     let c = if delta >= 0.0 { arr.contours(|m| m[0] || m[1] || m[2]) } else { arr.contours(|m| m[0] && !m[1]) };

crates/pathops/src/pathfinder.rs

@@ -69,7 +69,7 @@ impl Region {
 
 fn arrangement(shapes: &[Shape]) -> Option<Arrangement> {
     let bps: Vec<(BezPath, FillRule)> = shapes.iter().map(|s| (fill_bezpath(&s.path), s.rule)).collect();
-    Arrangement::new(&bps).ok()
+    Arrangement::new(bps).ok()
 }
 
 fn topmost(m: &[bool]) -> Option<usize> {

crates/pathops/tests/common/mod.rs

@@ -0,0 +1,20 @@
+//! Shared property-test configuration for the pathops suites.
+
+use proptest::test_runner::{Config, RngSeed};
+
+/// Seed the property suites run with unless `PROPTEST_RNG_SEED` is set.
+const SEED: u64 = 0x5eed_9a7f;
+
+/// Reproducible property-test config: `cases` cases from a fixed seed, so a CI run can never fail
+/// at random. Explore further locally with `PROPTEST_CASES=5000 PROPTEST_RNG_SEED=<n>`; keep
+/// every minimal failure found that way as a test in `regressions.rs`.
+pub fn config(cases: u32) -> Config {
+    let mut c = Config { failure_persistence: None, ..Config::default() };
+    if std::env::var_os("PROPTEST_CASES").is_none() {
+        c.cases = cases;
+    }
+    if std::env::var_os("PROPTEST_RNG_SEED").is_none() {
+        c.rng_seed = RngSeed::Fixed(SEED);
+    }
+    c
+}

crates/pathops/tests/prop_pathops.rs

@@ -1,6 +1,8 @@
 //! Path-operation properties on concave polygons, ellipses and random smooth curves (complements
 //! `props.rs`, which uses rectangles and circles).
 
+mod common;
+
 use proptest::prelude::*;
 use vectorcraft_geom::{FillRule, PathData, Point, Rect};
 use vectorcraft_pathops::*;
@@ -19,7 +21,7 @@ fn near_edge(p: &PathData, q: Point, eps: f64) -> bool {
 }
 
 proptest! {
-    #![proptest_config(ProptestConfig { cases: 128, failure_persistence: None, ..ProptestConfig::default() })]
+    #![proptest_config(common::config(128))]
 
     /// Area of a simple polygon equals its shoelace area; normalize keeps it.
     #[test]

crates/pathops/tests/props.rs

@@ -1,3 +1,5 @@
+mod common;
+
 use std::f64::consts::PI;
 
 use kurbo::{ParamCurve, Point, Rect, Shape as _};
@@ -20,7 +22,7 @@ fn close(a: f64, b: f64, rel: f64) -> bool {
 }
 
 proptest! {
-    #![proptest_config(ProptestConfig { cases: 48, failure_persistence: None, ..ProptestConfig::default() })]
+    #![proptest_config(common::config(48))]
 
     #[test]
     fn inclusion_exclusion(a in arb_shape(), b in arb_shape()) {

crates/pathops/tests/regressions.rs

@@ -0,0 +1,186 @@
+//! Minimal failures found by the property suites (`props.rs`, `prop_pathops.rs`), kept as
+//! deterministic tests.
+
+use kurbo::{BezPath, Shape as _};
+use vectorcraft_geom::{Anchor, FillRule, PathData, SubPath};
+use vectorcraft_pathops::*;
+
+const NZ: FillRule = FillRule::NonZero;
+
+/// A closed path from `[x, y, in_x, in_y, out_x, out_y]` anchors.
+fn curved(anchors: &[[f64; 6]]) -> PathData {
+    let a = |v: &[f64; 6]| Anchor::with_handles((v[0], v[1]).into(), (v[2], v[3]).into(), (v[4], v[5]).into());
+    PathData::single(SubPath::new(anchors.iter().map(a).collect(), true))
+}
+
+fn poly(pts: &[(f64, f64)]) -> PathData {
+    PathData::single(SubPath::polyline(&pts.iter().map(|&p| p.into()).collect::<Vec<_>>(), true))
+}
+
+fn close(a: f64, b: f64, rel: f64) -> bool {
+    (a - b).abs() <= rel * a.abs().max(b.abs()).max(1.0)
+}
+
+/// Area identities between the four boolean operations of `a` and `b`.
+fn assert_identities(a: &PathData, b: &PathData) {
+    let op = |p: &PathData, q: &PathData, op| area(&boolean(p, NZ, q, NZ, op), NZ);
+    let (aa, ab) = (area(a, NZ), area(b, NZ));
+    let (u, i, x) = (op(a, b, BoolOp::Union), op(a, b, BoolOp::Intersect), op(a, b, BoolOp::Xor));
+    let (d1, d2) = (op(a, b, BoolOp::Difference), op(b, a, BoolOp::Difference));
+    let tol = 1e-6;
+    assert!(close(d1 + d2, x, tol), "a−b {d1} + b−a {d2} != xor {x}");
+    assert!(close(u - i, x, tol), "union {u} − intersection {i} != xor {x}");
+    assert!(close(u + i, aa + ab, tol), "union {u} + intersection {i} != {aa} + {ab}");
+    assert!(close(d1 + i, aa, tol), "a−b {d1} + intersection {i} != a {aa}");
+    assert!(close(d2 + i, ab, tol), "b−a {d2} + intersection {i} != b {ab}");
+}
+
+/// A rectangle edge cutting a circle near its top or bottom: the cut leaves an edge that is
+/// horizontal up to the intersection's rounding (a rise of ~1e-12), which the sweep used to
+/// mis-order against the circle's arc, dropping a corner (`difference_and_xor_consistent`).
+#[test]
+fn rect_edge_through_circle_cap() {
+    let cases = [
+        (
+            curved(&[
+                [19.32739846439343, 15.797214477812503, 19.32739846439343, 14.207108832129435, 19.32739846439343, 17.38732012349557],
+                [17.360120998998184, 19.882308227700324, 18.60331565324336, 18.890893573738936, 16.116926344753008, 20.873722881661713],
+                [12.939688659402073, 20.89124306319632, 14.489927036696397, 21.245074856562926, 11.389450282107749, 20.537411269829715],
+                [9.394777157707882, 18.064270474652044, 10.084698140848792, 19.49690615808836, 8.704856174566974, 16.631634791215724],
+                [9.394777157707882, 13.530158480972963, 8.704856174566974, 14.962794164409281, 10.08469814084879, 12.097522797536644],
+                [12.939688659402071, 10.703185892428685, 11.389450282107747, 11.05701768579529, 14.489927036696395, 10.34935409906208],
+                [17.360120998998184, 11.712120727924681, 16.116926344753004, 10.720706073963292, 18.60331565324336, 12.70353538188607],
+            ]),
+            poly(&[
+                (0.0, 20.17376246815288),
+                (37.21494402256868, 20.17376246815288),
+                (37.21494402256868, 21.17376246815288),
+                (0.0, 21.17376246815288),
+            ]),
+        ),
+        (
+            curved(&[
+                [23.57990847002715, 27.750243942757283, 23.57990847002715, 26.818683968734558, 23.57990847002715, 28.68180391678001],
+                [22.276172537520765, 30.00838081751153, 23.08292714017322, 29.542600830500167, 21.46941793486831, 30.474160804522892],
+                [19.668700672507995, 30.00838081751153, 20.47545527516045, 30.474160804522892, 18.86194606985554, 29.542600830500167],
+                [18.36496474000161, 27.750243942757283, 18.36496474000161, 28.68180391678001, 18.36496474000161, 26.818683968734558],
+                [19.668700672507995, 25.492107068003037, 18.86194606985554, 25.9578870550144, 20.475455275160446, 25.026327080991674],
+                [22.276172537520765, 25.492107068003037, 21.46941793486831, 25.02632708099167, 23.08292714017322, 25.9578870550144],
+            ]),
+            poly(&[
+                (0.0, 30.068083063105394),
+                (34.105297058743425, 30.068083063105394),
+                (34.105297058743425, 31.068083063105394),
+                (0.0, 31.068083063105394),
+            ]),
+        ),
+        (
+            poly(&[
+                (9.700636443853261, 23.969159526578938),
+                (26.555106262718724, 23.969159526578938),
+                (26.555106262718724, 54.82336332868054),
+                (9.700636443853261, 54.82336332868054),
+            ]),
+            curved(&[
+                [27.958532449581877, 19.508138986349657, 27.958532449581877, 17.819765427048214, 27.958532449581877, 21.1965125456511],
+                [25.869677998055128, 23.845689949270323, 27.1897016008836, 22.79300625331794, 24.549654395226657, 24.898373645222705],
+                [21.176064650340088, 24.916976567196787, 22.822107159492838, 25.292675028479692, 19.530022141187338, 24.541278105913882],
+                [17.412078989530123, 21.91529816656589, 18.144636822467053, 23.436470180886086, 16.67952115659319, 20.39412615224569],
+                [17.412078989530123, 17.100979806133427, 16.67952115659319, 18.622151820453624, 18.14463682246705, 15.579807791813227],
+                [21.176064650340088, 14.099301405502526, 19.530022141187338, 14.474999866785431, 22.822107159492838, 13.723602944219623],
+                [25.869677998055128, 15.17058802342899, 24.549654395226657, 14.117904327476609, 27.1897016008836, 16.223271719381373],
+            ]),
+        ),
+        (
+            poly(&[
+                (0.0, 36.673784808936425),
+                (36.25829503976694, 36.673784808936425),
+                (36.25829503976694, 61.69546684624716),
+                (0.0, 61.69546684624716),
+            ]),
+            curved(&[
+                [38.42179318629935, 31.567929837130343, 38.42179318629935, 29.68809407962987, 38.42179318629935, 33.447765594630816],
+                [36.096061852400304, 36.39736071465043, 37.56577662948331, 35.225302290679494, 34.6263470753173, 37.56941913862136],
+                [30.870190981402537, 37.59013163896124, 32.70289533170721, 38.00843444740472, 29.037486631097863, 37.171828830517754],
+                [26.679367927506128, 34.24806177596396, 27.49499809489732, 35.94173527024165, 25.863737760114933, 32.55438828168628],
+                [26.679367927506128, 28.887797898296725, 25.863737760114933, 30.58147139257441, 27.49499809489732, 27.194124404019036],
+                [30.870190981402533, 25.545728035299447, 29.037486631097863, 25.96403084374293, 32.70289533170721, 25.127425226855966],
+                [36.096061852400304, 26.738498959610258, 34.6263470753173, 25.566440535639327, 37.56577662948331, 27.91055738358119],
+            ]),
+        ),
+        (
+            poly(&[
+                (13.843112028864478, 38.445768817380575),
+                (40.14690269607463, 38.445768817380575),
+                (40.14690269607463, 39.445768817380575),
+                (13.843112028864478, 39.445768817380575),
+            ]),
+            curved(&[
+                [39.74299892595637, 33.15627685561371, 39.74299892595637, 31.327165913276996, 39.74299892595637, 34.98538779795043],
+                [37.48002429383932, 37.855392099105174, 38.91008081348493, 36.714960080090016, 36.04996777419372, 38.99582411812034],
+                [32.39516645174228, 39.01597771976541, 34.17841776390354, 39.42299319496389, 30.61191513958101, 38.608962244566925],
+                [28.317427066961493, 35.76408907360503, 29.111048561883386, 37.412061088590065, 27.523805572039596, 34.11611705861999],
+                [28.317427066961493, 30.548464637622395, 27.523805572039596, 32.19643665260743, 29.111048561883386, 28.90049262263736],
+                [32.39516645174227, 27.296575991462014, 30.61191513958101, 27.703591466660495, 34.17841776390354, 26.889560516263533],
+                [37.48002429383932, 28.457161612122245, 36.04996777419371, 27.316729593107084, 38.91008081348493, 29.597593631137407],
+            ]),
+        ),
+        (
+            poly(&[
+                (29.861827337105368, 18.136938945855064),
+                (66.57181423170378, 18.136938945855064),
+                (66.57181423170378, 20.881417410214578),
+                (29.861827337105368, 20.881417410214578),
+            ]),
+            curved(&[
+                [50.05371460858578, 9.878937912496232, 50.05371460858578, 6.762209599042626, 50.05371460858578, 12.995666225949838],
+                [46.19770076858542, 17.886033156736026, 48.634457086342934, 15.94278483813333, 43.76094445082791, 19.829281475338718],
+                [37.53331631845236, 19.863622367086386, 40.57190174592566, 20.55715966228417, 34.494730890979056, 19.1700850718886],
+                [30.585019470370003, 14.32254053268522, 31.937317204824815, 17.130615712816933, 29.23272173591519, 11.51446535255351],
+                [30.58501947037, 5.435335292307244, 29.23272173591519, 8.243410472438956, 31.937317204824815, 2.6272601121755335],
+                [37.53331631845236, -0.10574654209391987, 34.494730890979056, 0.5877907531038655, 40.57190174592566, -0.7992838372917053],
+                [46.19770076858542, 1.8718426682564377, 43.76094445082791, -0.07140565034625546, 48.634457086342934, 3.815090986859132],
+            ]),
+        ),
+    ];
+    for (a, b) in &cases {
+        assert_identities(a, b);
+    }
+}
+
+/// The root cause in isolation: a strip whose top edge rises by less than the sweep tolerance
+/// (1e-6 here) and ends level with the point where an arc meets a horizontal edge of a second
+/// contour. Every tilt must keep the strip's full area.
+#[test]
+fn almost_horizontal_edge_meeting_an_arc_end() {
+    for dy in [1e-14, 1e-12, 1e-10, 1e-8, 1e-7, 5e-7, 9e-7] {
+        let svg = format!("M10,10 L20,10 L20,16 C19,18 18,19 17,20 L11,20 Z M0,21 L0,20 L11,20 L13,20.9 L16,{} L37,20 L37,21 Z", 20.0 - dy);
+        let bp = BezPath::from_svg(&svg).unwrap();
+        // The contours overlap by at most a `dy`-thin sliver (and levelling the edge moves it by
+        // `dy`), so the filled area is the plain (Green's theorem) area up to `dy` × its length.
+        let want = bp.area().abs();
+        let got = area(&PathData::from_bezpath(&bp), NZ);
+        assert!((got - want).abs() <= dy * 40.0, "rise {dy}: area {got}, want {want}");
+    }
+}
+
+/// `polygon_area_matches_shoelace` once drew this "star" polygon, whose wide angular step made it
+/// cross itself (the strategy now rejects those): its shoelace area cancels across the crossing,
+/// while the filled area is that of the two loops.
+#[test]
+fn self_crossing_quad_fills_both_loops() {
+    let [p0, p1, p2, p3]: [kurbo::Point; 4] = [
+        (43.48211493626725, 43.58815768474916),
+        (39.45699865107951, 58.097908614334145),
+        (56.35982114265447, 23.14196428649071),
+        (72.23596128912416, 39.99999999999999),
+    ]
+    .map(Into::into);
+    // Where p1→p2 crosses p3→p0.
+    let (d1, d2) = (p2 - p1, p0 - p3);
+    let x = p1 + d1 * ((p3 - p1).cross(d2) / d1.cross(d2));
+    let tri = |a: kurbo::Point, b: kurbo::Point, c: kurbo::Point| (b - a).cross(c - a).abs() / 2.0;
+    let want = tri(p0, p1, x) + tri(x, p2, p3);
+    let got = area(&poly(&[p0, p1, p2, p3].map(|p| (p.x, p.y))), NZ);
+    assert!(close(got, want, 1e-9), "area {got}, want {want}");
+}

crates/testkit/src/strategies.rs

@@ -19,9 +19,14 @@ pub fn arb_rect() -> impl Strategy<Value = Rect> {
 }
 
 /// A simple (non self-intersecting) star-shaped polygon: points at increasing angles around a
-/// centre with random radii.
+/// centre with random radii. Every angular step stays below half a turn: a wider one puts the
+/// centre outside and can make the polygon cross itself.
 pub fn arb_star_polygon() -> impl Strategy<Value = Vec<Point>> {
-    (arb_point(40.0, 160.0), prop::collection::vec((0.2f64..1.0, 5.0f64..40.0), 3..12)).prop_map(|(c, v)| {
+    let steps = prop::collection::vec((0.2f64..1.0, 5.0f64..40.0), 3..12).prop_filter("an angular step of half a turn or more", |v| {
+        let total: f64 = v.iter().map(|(a, _)| a).sum();
+        v.iter().all(|(a, _)| 2.0 * a < total)
+    });
+    (arb_point(40.0, 160.0), steps).prop_map(|(c, v)| {
         let total: f64 = v.iter().map(|(a, _)| a).sum();
         let mut ang = 0.0;
         v.iter()