From 3e32b3a33c477eeaee7722ef4eae6ff4bd81b388 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 12:41:36 -0600 Subject: [PATCH 01/20] fix(overlay): split result lobes that meet at a single point MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `trace_rings` closed a ring only when the walk returned to the node it started from. Where two lobes of the result touch at one point, the traversal passes through that point twice, so carrying on to the seed spliced the lobes into a single self-touching ring — geometry that `is_valid` rejects, and that `boost::geometry::intersection` never produces. Close the ring at any node already on the current walk instead, emit it, and carry on from there. This is the rule Boost's own traversal follows, and the one the `geometry/correct.hpp` dissolve uses. Verified against C++ Boost 1.83, which returns two valid polygons for the case now covered by `intersection_splits_lobes_meeting_at_a_point`: box (0,0)-(10,10) ring (5,-1) (6,-2) (2,-3) (1,1) (3,2) (2,6) (4,7) (5,0) (6,1) (7,-1) before: one ring through (5,0) twice, area 16.125 but self-touching after: (1.25,0)(1,1)(3,2)(2,6)(4,7)(5,0) and (5,0)(6,1)(6.5,0) Found while scoring this crate against C++ Boost on real polygon/box clips taken from a tilemaker run: it takes those from 3/10 to 6/10 matching Boost's ring set, with no change to the rest of the suite. --- .../geometry-overlay/src/operation/areal.rs | 51 ++++++++++++++----- .../geometry-overlay/tests/overlay_parity.rs | 43 ++++++++++++++++ 2 files changed, 81 insertions(+), 13 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index c6e6123..89887f9 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -353,29 +353,54 @@ where used[edge_index] = true; let edge = edges[edge_index]; node_indices.push(edge.end); + + // A node the walk has already stood on closes a ring right here, + // not only when the walk returns to the seed. Where two lobes of + // the result meet at a single point, the traversal passes through + // that point twice; carrying on to the seed splices the lobes into + // one self-touching ring, which is not a valid polygon and is not + // what `boost::geometry::intersection` returns. Cut the loop out, + // keep the path up to that node, and carry on walking. + if let Some(start) = node_indices[..node_indices.len() - 1] + .iter() + .position(|&index| index == edge.end) + { + let loop_nodes = node_indices.split_off(start); + node_indices.push(edge.end); + push_ring(&mut rings, nodes, &loop_nodes, tolerance); + } + if edge.end == first { break; } edge_index = next_edge(nodes, edges, &used, edge).ok_or(OverlayError::Unsupported)?; } debug_assert_eq!(node_indices.last().copied(), Some(first)); - let area = node_indices.windows(2).fold(0.0, |sum, pair| { - let a = nodes[pair[0]].coordinate; - let b = nodes[pair[1]].coordinate; - sum + a.x * b.y - b.x * a.y - }) * 0.5; - if area.abs() > tolerance * tolerance { - rings.push(Ring::from_vec( - node_indices - .into_iter() - .map(|index| nodes[index].point) - .collect(), - )); - } } Ok(rings) } +/// Append one traced cycle, dropping it when it encloses no area. +fn push_ring

(rings: &mut Vec>, nodes: &[Node

], node_indices: &[usize], tolerance: f64) +where + P: Point + Copy, + P::Scalar: Into, +{ + let area = node_indices.windows(2).fold(0.0, |sum, pair| { + let a = nodes[pair[0]].coordinate; + let b = nodes[pair[1]].coordinate; + sum + a.x * b.y - b.x * a.y + }) * 0.5; + if area.abs() > tolerance * tolerance { + rings.push(Ring::from_vec( + node_indices + .iter() + .map(|&index| nodes[index].point) + .collect(), + )); + } +} + fn next_edge

(nodes: &[Node

], edges: &[Edge], used: &[bool], incoming: Edge) -> Option where P: Point, diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 9744b6c..6bb928d 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -170,3 +170,46 @@ fn disjoint_all_four_ops() { close(area(&difference(&a, &b).unwrap()), 1.0); close(area(&sym_difference(&a, &b).unwrap()), 2.0); } + +// ---- Result lobes meeting at one point ------------------------------- + +/// A ring that dips out of the clip box twice and grazes its edge at a +/// single vertex in between. The intersection is two polygons that touch +/// at `(5, 0)`; splicing them into one ring through that point would be +/// an invalid self-touching polygon. +/// +/// C++ Boost (`boost::geometry::intersection`, 1.83) returns, in order: +/// `(1.25,0) (1,1) (3,2) (2,6) (4,7) (5,0) (1.25,0)` and +/// `(5,0) (6,1) (6.5,0) (5,0)`, with `is_valid` true. +#[test] +fn intersection_splits_lobes_meeting_at_a_point() { + let subject: Polygon

= polygon![[ + (5.0, -1.0), + (6.0, -2.0), + (2.0, -3.0), + (1.0, 1.0), + (3.0, 2.0), + (2.0, 6.0), + (4.0, 7.0), + (5.0, 0.0), + (6.0, 1.0), + (7.0, -1.0), + (5.0, -1.0) + ]]; + let clip = square(0.0, 0.0, 10.0); + + let result = intersection(&subject, &clip).unwrap(); + assert_eq!( + result.polygons().count(), + 2, + "the two lobes must stay separate polygons" + ); + + // 15.375 for the large lobe plus 0.75 for the small one; C++ Boost + // reports the same 16.125 for this input. + close(area(&result), 16.125); + + let mut sizes: Vec = result.polygons().map(|pg| pg.exterior().0.len()).collect(); + sizes.sort_unstable(); + assert_eq!(sizes, [4, 7], "each lobe keeps its own closed ring"); +} From f598fcc9ea48ef557848cc0bab69db22cb9946be Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 12:54:12 -0600 Subject: [PATCH 02/20] fix(overlay): report wrong orientation before self-intersection MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `is_valid_ring_with` tested orientation last, so a ring that was both wound the wrong way and self-intersecting came back `SelfIntersection`. Boost reports `failure_wrong_orientation` for it: orientation is checked first, and only a correctly wound ring goes on to be reported as self-intersecting. Spikes still precede both. The distinction is not cosmetic for callers that branch on the code — tilemaker's `buildWayGeometry` redoes its clip on `failure_self_intersections` and leaves the geometry alone on `failure_wrong_orientation`. Measured against Boost 1.83 on 2000 random rings (clean, bow-tie, figure-of-eight and meander shapes, both windings): 203 of the first 400 disagreed before, 0 of 2000 after. The order the new `orientation_is_reported_before_self_intersection` test pins: ccw + real spike failure=12 area=-16 ccw + selfint, no spike failure=22 area=-17 cw self-int, +area failure=21 area=+98 bowtie failure=22 area=+0 Two existing bow-tie assertions expected the old order; a bow-tie's lobes cancel to zero signed area, so Boost fails it on orientation. Updated with the measured value. --- crates/geometry-overlay/src/validity.rs | 96 +++++++++++++++++-- .../tests/spatial_predicates.rs | 6 +- 2 files changed, 91 insertions(+), 11 deletions(-) diff --git a/crates/geometry-overlay/src/validity.rs b/crates/geometry-overlay/src/validity.rs index 005648b..3d43fb7 100644 --- a/crates/geometry-overlay/src/validity.rs +++ b/crates/geometry-overlay/src/validity.rs @@ -540,15 +540,18 @@ where return Err(ValidityFailure::Spikes); } - // No non-adjacent edge may intersect another. - if has_self_intersection(&pts) { - return Err(ValidityFailure::SelfIntersection); - } - // Orientation — Boost's failure_wrong_orientation. The strategy area // already folds the declared PointOrder: a correctly wound exterior // is positive, a correctly wound hole negative. Zero area // (degenerate) fails either way. + // + // Checked *before* self-intersection, because that is the order Boost + // reports in: a counter-clockwise ring that also crosses itself comes + // back `failure_wrong_orientation`, and only a correctly wound ring + // goes on to be reported as `failure_self_intersections`. Spikes still + // come first — a wrongly wound ring carrying a spike is + // `failure_spikes`. Verified against Boost 1.83; see + // `orientation_is_reported_before_self_intersection`. let area = ShoelaceArea.area(ring); let zero = ::ZERO; let properly_oriented = if is_interior { @@ -560,6 +563,11 @@ where return Err(ValidityFailure::WrongOrientation); } + // No non-adjacent edge may intersect another. + if has_self_intersection(&pts) { + return Err(ValidityFailure::SelfIntersection); + } + Ok(()) } @@ -842,7 +850,10 @@ mod tests { is_valid_ring(&huge_bowtie), Err(ValidityFailure::CoordinateOutOfRange) ); - // The in-range analogue is still caught as a self-intersection. + // The in-range analogue is still caught — as WrongOrientation, + // which is what Boost 1.83 reports for it: a bow-tie has zero + // signed area, so it fails the orientation test before the + // self-intersection test is ever reached. let small_bowtie: Ring

= Ring::from_vec(vec![ P::new(0.0, 0.0), P::new(2.0, 2.0), @@ -852,7 +863,7 @@ mod tests { ]); assert_eq!( is_valid_ring(&small_bowtie), - Err(ValidityFailure::SelfIntersection) + Err(ValidityFailure::WrongOrientation) ); } @@ -869,7 +880,11 @@ mod tests { #[test] fn self_intersecting_bowtie() { - // A "bow-tie" quadrilateral whose diagonals cross. + // A "bow-tie" quadrilateral whose diagonals cross. Its two lobes + // cancel, so its signed area is zero and Boost reports the + // orientation failure rather than the crossing: + // + // bowtie (zero area) valid=0 failure=22 area=+0 let r: Ring

= Ring::from_vec(vec![ P::new(0.0, 0.0), P::new(2.0, 2.0), @@ -877,7 +892,70 @@ mod tests { P::new(0.0, 2.0), P::new(0.0, 0.0), ]); - assert_eq!(is_valid_ring(&r), Err(ValidityFailure::SelfIntersection)); + assert_eq!(is_valid_ring(&r), Err(ValidityFailure::WrongOrientation)); + } + + /// The order the ring checks report in, pinned against Boost 1.83. + /// + /// Boost runs spikes, then orientation, then self-intersection, and + /// stops at the first failure. Getting the last two the wrong way round + /// is not cosmetic: a caller that branches on the code — tilemaker's + /// `buildWayGeometry` re-clips on `failure_self_intersections` but not + /// on `failure_wrong_orientation` — takes a different path. + #[test] + fn orientation_is_reported_before_self_intersection() { + // ccw + selfint no spike valid=0 failure=22 area=-17 + let wound_wrong_and_crossing: Ring

= Ring::from_vec(vec![ + P::new(0.0, 0.0), + P::new(4.0, 0.0), + P::new(4.0, 4.0), + P::new(0.0, 4.0), + P::new(0.0, 0.0), + P::new(1.0, -1.0), + P::new(3.0, -3.0), + P::new(1.0, -3.0), + P::new(3.0, -1.0), + P::new(0.0, 0.0), + ]); + assert_eq!( + is_valid_ring(&wound_wrong_and_crossing), + Err(ValidityFailure::WrongOrientation) + ); + + // cw self-int, +area valid=0 failure=21 area=+98 + let wound_right_and_crossing: Ring

= Ring::from_vec(vec![ + P::new(0.0, 0.0), + P::new(0.0, 10.0), + P::new(10.0, 10.0), + P::new(10.0, 0.0), + P::new(0.0, 0.0), + P::new(3.0, -4.0), + P::new(7.0, -4.0), + P::new(3.0, -8.0), + P::new(7.0, -8.0), + P::new(0.0, 0.0), + ]); + assert_eq!( + is_valid_ring(&wound_right_and_crossing), + Err(ValidityFailure::SelfIntersection) + ); + + // ccw + real spike valid=0 failure=12 area=-16 + // Spikes still win over orientation. + let wound_wrong_with_spike: Ring

= Ring::from_vec(vec![ + P::new(0.0, 0.0), + P::new(2.0, 0.0), + P::new(2.0, -2.0), + P::new(2.0, 0.0), + P::new(4.0, 0.0), + P::new(4.0, 4.0), + P::new(0.0, 4.0), + P::new(0.0, 0.0), + ]); + assert_eq!( + is_valid_ring(&wound_wrong_with_spike), + Err(ValidityFailure::Spikes) + ); } #[test] diff --git a/crates/geometry-overlay/tests/spatial_predicates.rs b/crates/geometry-overlay/tests/spatial_predicates.rs index cec6427..2b610ea 100644 --- a/crates/geometry-overlay/tests/spatial_predicates.rs +++ b/crates/geometry-overlay/tests/spatial_predicates.rs @@ -86,11 +86,13 @@ fn valid_and_invalid_polygons() { let good: Polygon

= square(0.0, 0.0, 3.0); assert!(is_valid_polygon(&good).is_ok()); - // A self-intersecting "bow-tie" exterior. + // A self-intersecting "bow-tie" exterior. Its lobes cancel, so its + // signed area is zero and Boost 1.83 reports the orientation failure + // rather than the crossing — orientation is checked first. let bowtie: Polygon

= polygon![[(0.0, 0.0), (2.0, 2.0), (2.0, 0.0), (0.0, 2.0), (0.0, 0.0)]]; assert_eq!( is_valid_polygon(&bowtie), - Err(ValidityFailure::SelfIntersection) + Err(ValidityFailure::WrongOrientation) ); } From 480d277cbc70bddfb4a6943a35016fddaf66b311 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 12:59:07 -0600 Subject: [PATCH 03/20] fix(algorithm): collapse repeated vertices in remove_spikes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The predicate required the dot product to be strictly negative, so a zero-length step never qualified and every repeated vertex survived. Boost's is `point_is_spike_or_equal` — the `or_equal` half is what collapses a duplicate, and `remove_spikes` needs it to clean up after itself: peeling the apex off `(4,0) (6,0) (4,0)` leaves `(4,0) (4,0)` adjacent, which a spike-only predicate then leaves in the output. `dot <= 0` cannot over-match. Two non-zero steps that are both parallel (`cross == 0`) and perpendicular (`dot == 0`) do not exist, so the new arm fires only when one step has zero length. Boost 1.83 on a clockwise `model::polygon`, all five now covered by `repeated_vertices_are_collapsed`: consecutive dup -> (0,0) (0,4) (4,4) (4,0) (0,0) dup at start -> (0,0) (0,4) (4,4) (4,0) (0,0) triple dup -> (0,0) (0,4) (4,4) (4,0) (0,0) real spike -> (0,0) (0,4) (4,4) (4,0) (0,0) dup + spike -> (0,0) (0,4) (4,4) (4,0) (0,0) `is_valid`'s own spike predicate stays strict and keeps `dot < 0`: Boost's default validity policy accepts duplicate points, so the two answer different questions. Noted where they diverge. Measured against Boost on 2000 random rings: 1 in 400 disagreed before, 0 of 2000 after. --- .../geometry-algorithm/src/remove_spikes.rs | 136 ++++++++++++++++-- crates/geometry-overlay/src/validity.rs | 14 +- 2 files changed, 131 insertions(+), 19 deletions(-) diff --git a/crates/geometry-algorithm/src/remove_spikes.rs b/crates/geometry-algorithm/src/remove_spikes.rs index b9d8f52..37d7b44 100644 --- a/crates/geometry-algorithm/src/remove_spikes.rs +++ b/crates/geometry-algorithm/src/remove_spikes.rs @@ -1,11 +1,15 @@ //! `remove_spikes(&mut g)` — drop collinear-and-reversed vertices. //! //! Mirrors `boost::geometry::remove_spikes` from -//! `boost/geometry/algorithms/remove_spikes.hpp`. A spike is a triple -//! `(a, b, c)` where `(b-a) × (c-b) == 0` (collinear) AND -//! `(b-a) · (c-b) < 0` (reversed). The middle vertex `b` is removed; -//! the walk repeats until no spike remains, because collapsing one -//! spike can create a new one at the now-adjacent pair. +//! `boost/geometry/algorithms/remove_spikes.hpp`. The predicate is +//! Boost's `point_is_spike_or_equal`, and the `or_equal` half carries +//! its weight: a triple `(a, b, c)` qualifies when `(b-a) × (c-b) == 0` +//! (collinear) and `(b-a) · (c-b) <= 0`, which covers both a reversal +//! and a zero-length step — that is, a repeated vertex. The middle +//! vertex `b` is removed; the walk repeats until nothing qualifies, +//! because collapsing one spike can create a new one at the +//! now-adjacent pair, and peeling a spike off a ring routinely leaves a +//! repeated vertex behind. //! //! Per-kind: //! * `Linestring`, `Ring` → spike-walk the backing `Vec

` @@ -34,9 +38,19 @@ pub trait RemoveSpikes { fn remove_spikes(&mut self); } -/// True iff `b` is a spike between `a` and `c` (2D cross `== 0` AND -/// dot `< 0`). -fn is_spike_2d(a: &P, b: &P, c: &P) -> bool { +/// True iff `b` is a spike between `a` and `c`, **or** duplicates one of +/// them: 2D cross `== 0` and dot `<= 0`. +/// +/// Mirrors `detail::point_is_spike_or_equal` +/// (`algorithms/detail/point_is_spike_or_equal.hpp`). Requiring `dot < 0` +/// instead would leave every repeated vertex in place, including the ones +/// this function creates: removing the apex of `(4,0) (6,0) (4,0)` leaves +/// `(4,0) (4,0)` adjacent, and Boost collapses that. +/// +/// `dot <= 0` cannot over-match. Two non-zero vectors that are both +/// parallel (`cross == 0`) and perpendicular (`dot == 0`) do not exist, so +/// the equality arm fires only when one of the steps has zero length. +fn is_spike_or_equal_2d(a: &P, b: &P, c: &P) -> bool { let ux = b.get::<0>() - a.get::<0>(); let uy = b.get::<1>() - a.get::<1>(); let vx = c.get::<0>() - b.get::<0>(); @@ -44,7 +58,7 @@ fn is_spike_2d(a: &P, b: &P, c: &P) -> bool { let cross = ux * vy - uy * vx; let dot = ux * vx + uy * vy; let zero = ::ZERO; - cross == zero && dot < zero + cross == zero && dot <= zero } fn walk_spikes(pts: &mut alloc::vec::Vec

) { @@ -53,7 +67,7 @@ fn walk_spikes(pts: &mut alloc::vec::Vec

) { changed = false; let mut i = 1; while i + 1 < pts.len() { - if is_spike_2d(&pts[i - 1], &pts[i], &pts[i + 1]) { + if is_spike_or_equal_2d(&pts[i - 1], &pts[i], &pts[i + 1]) { pts.remove(i); changed = true; // Do not advance `i`: the new `pts[i]` (was `pts[i+1]`) @@ -105,12 +119,14 @@ fn walk_ring_spikes(pts: &mut alloc::vec::Vec

, closed: while found { found = false; // Spike at the first point: (prev = back-1, back, front). - while pts.len() >= 3 && is_spike_2d(&pts[pts.len() - 2], &pts[pts.len() - 1], &pts[0]) { + while pts.len() >= 3 + && is_spike_or_equal_2d(&pts[pts.len() - 2], &pts[pts.len() - 1], &pts[0]) + { pts.pop(); found = true; } // Spike at the second point: (back, front, front+1). - while pts.len() >= 3 && is_spike_2d(&pts[pts.len() - 1], &pts[0], &pts[1]) { + while pts.len() >= 3 && is_spike_or_equal_2d(&pts[pts.len() - 1], &pts[0], &pts[1]) { pts.remove(0); found = true; } @@ -160,8 +176,8 @@ mod tests { use super::remove_spikes; use geometry_cs::Cartesian; - use geometry_model::{Point2D, linestring}; - use geometry_trait::Linestring as _; + use geometry_model::{Point2D, Ring, linestring}; + use geometry_trait::{Linestring as _, Point as _, Ring as _}; type P = Point2D; @@ -306,4 +322,96 @@ mod tests { assert!(r.points().count() >= 4); assert_eq!(pts.first(), pts.last(), "ring must remain closed"); } + + /// Boost collapses a repeated vertex the same way it collapses a + /// spike — `point_is_spike_or_equal` covers both. Expected values from + /// `boost::geometry::remove_spikes` on a clockwise `model::polygon` + /// (Boost 1.83): + /// + /// ```text + /// consecutive dup -> (0,0) (0,4) (4,4) (4,0) (0,0) + /// dup at start -> (0,0) (0,4) (4,4) (4,0) (0,0) + /// triple dup -> (0,0) (0,4) (4,4) (4,0) (0,0) + /// real spike -> (0,0) (0,4) (4,4) (4,0) (0,0) + /// dup + spike -> (0,0) (0,4) (4,4) (4,0) (0,0) + /// ``` + /// + /// The `real spike` row is the one that shows why: removing the apex + /// of `(4,0) (6,0) (4,0)` leaves `(4,0) (4,0)` adjacent, so a + /// spike-only predicate makes duplicates out of its own output. + #[test] + fn repeated_vertices_are_collapsed() { + let square = [(0.0, 0.0), (0.0, 4.0), (4.0, 4.0), (4.0, 0.0), (0.0, 0.0)]; + + for (name, input) in [ + ( + "consecutive dup", + vec![ + (0.0, 0.0), + (0.0, 4.0), + (4.0, 4.0), + (4.0, 4.0), + (4.0, 0.0), + (0.0, 0.0), + ], + ), + ( + "dup at start", + vec![ + (0.0, 0.0), + (0.0, 0.0), + (0.0, 4.0), + (4.0, 4.0), + (4.0, 0.0), + (0.0, 0.0), + ], + ), + ( + "triple dup", + vec![ + (0.0, 0.0), + (0.0, 4.0), + (4.0, 4.0), + (4.0, 4.0), + (4.0, 4.0), + (4.0, 0.0), + (0.0, 0.0), + ], + ), + ( + "real spike", + vec![ + (0.0, 0.0), + (0.0, 4.0), + (4.0, 4.0), + (4.0, 0.0), + (6.0, 0.0), + (4.0, 0.0), + (0.0, 0.0), + ], + ), + ( + "dup + spike", + vec![ + (0.0, 0.0), + (0.0, 4.0), + (4.0, 4.0), + (4.0, 4.0), + (4.0, 0.0), + (6.0, 0.0), + (4.0, 0.0), + (0.0, 0.0), + ], + ), + ] { + let mut ring: Ring

= + Ring::from_vec(input.iter().map(|&(x, y)| P::new(x, y)).collect()); + remove_spikes(&mut ring); + let pts: Vec<(f64, f64)> = ring + .points() + .map(|p| (p.get::<0>(), p.get::<1>())) + .collect(); + assert_eq!(pts, square, "{name}"); + } + } } diff --git a/crates/geometry-overlay/src/validity.rs b/crates/geometry-overlay/src/validity.rs index 3d43fb7..818a1ab 100644 --- a/crates/geometry-overlay/src/validity.rs +++ b/crates/geometry-overlay/src/validity.rs @@ -768,11 +768,15 @@ where } /// `true` iff `b` is a spike between `a` and `c`: collinear -/// (`cross == 0`) and folding back (`dot < 0`). Same 2-D kernel as -/// `geometry_algorithm::remove_spikes::is_spike_2d` (private there); -/// duplicated locally rather than widening that crate's public -/// surface — if a third consumer appears, hoist the predicate into a -/// shared home per the aggregate-slicing rules. +/// (`cross == 0`) and folding back (`dot < 0`). +/// +/// Deliberately **stricter** than +/// `geometry_algorithm::remove_spikes::is_spike_or_equal_2d`, which also +/// fires on a zero-length step. `remove_spikes` drops a repeated vertex; +/// `is_valid` does not reject one — Boost's default policy accepts +/// duplicates (see [`ValidityOptions::BOOST_DEFAULT`]), and a ring +/// carrying one is valid until `allow_duplicates` is turned off. The two +/// predicates answer different questions, so they are not shared. fn is_spike_triple(a: &P, b: &P, c: &P) -> bool where P::Scalar: CoordinateScalar, From ad97471685c5e0dbbab31366f224d7113bc870f8 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 13:19:08 -0600 Subject: [PATCH 04/20] feat(strategy): centroid for multi-polygons, accumulated as Boost does MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `centroid` was a compile error on a multi-polygon; Boost has had one all along. Adds `CartesianMultiPolygonCentroid` and the `MultiPolygonTag` picker impl. It runs one Bashein-Detmer accumulator over every ring of every member and divides once, which is `centroid_multi` in `boost/geometry/algorithms/centroid.hpp` — not a combine of per-member centroids weighted by area. The two are not the same function. A member with zero area has zero weight and vanishes from the weighted combine, but its numerator is not zero and Boost keeps it. `CartesianPolygonCentroid` had the same shape and the same gap, so it now accumulates too. Boost 1.83 on a polygon whose exterior is a bow-tie and whose hole has real area: bowtie exterior + hole -> (10.6667, 11) area=-4 weighted combine gives -> (11, 11) and across members: zero-area MP -> (0, 0) area=0 mixed MP -> (11.3333, 11) area=4 Both are the accumulated answer, and both are now covered by `a_zero_area_part_still_moves_the_centroid`. The existing polygon tests are unchanged: for a polygon with well-formed rings the two formulations agree identically, because the per-ring area cancels against the divisor. --- crates/geometry-strategy/src/centroid.rs | 143 ++++++++++++++++++++--- crates/geometry-strategy/src/lib.rs | 4 +- 2 files changed, 131 insertions(+), 16 deletions(-) diff --git a/crates/geometry-strategy/src/centroid.rs b/crates/geometry-strategy/src/centroid.rs index aea5566..662e6e1 100644 --- a/crates/geometry-strategy/src/centroid.rs +++ b/crates/geometry-strategy/src/centroid.rs @@ -31,11 +31,13 @@ use geometry_coords::CoordinateScalar; use geometry_cs::{CartesianFamily, CoordinateSystem}; -use geometry_tag::{BoxTag, LinestringTag, MultiPointTag, PolygonTag, RingTag, SameAs, SegmentTag}; +use geometry_tag::{ + BoxTag, LinestringTag, MultiPointTag, MultiPolygonTag, PolygonTag, RingTag, SameAs, SegmentTag, +}; use geometry_trait::{ Box as BoxTrait, Geometry, Linestring as LinestringTrait, MultiPoint as MultiPointTrait, - Point as PointTrait, PointMut, Polygon as PolygonTrait, Ring as RingTrait, - Segment as SegmentTrait, box_max, box_min, segment_end, segment_start, + MultiPolygon as MultiPolygonTrait, Point as PointTrait, PointMut, Polygon as PolygonTrait, + Ring as RingTrait, Segment as SegmentTrait, box_max, box_min, segment_end, segment_start, }; use crate::area::{AreaStrategy, ShoelaceArea}; @@ -83,6 +85,14 @@ pub struct CartesianRingCentroid; #[derive(Debug, Default, Clone, Copy)] pub struct CartesianPolygonCentroid; +/// Cartesian centroid for a [`geometry_trait::MultiPolygon`] — one +/// Bashein–Detmer accumulator over every ring of every member. +/// +/// Mirrors the multi-polygon arm of +/// `boost/geometry/algorithms/centroid.hpp`. +#[derive(Debug, Default, Clone, Copy)] +pub struct CartesianMultiPolygonCentroid; + /// Cartesian centroid for a [`geometry_trait::Linestring`] — length-weighted midpoint of /// each segment, summed and divided by total length. /// @@ -260,16 +270,15 @@ where fn centroid(&self, pg: &G) -> G::Point { let zero = ::Scalar::ZERO; - let mut sum_area = zero; + let mut sum_a2 = zero; let mut sum_x = zero; let mut sum_y = zero; let mut fold_ring = |ring: &G::Ring| { - let area = ShoelaceArea.area(ring); - let c = CartesianRingCentroid.centroid(ring); - sum_area = sum_area + area; - sum_x = sum_x + area * c.get::<0>(); - sum_y = sum_y + area * c.get::<1>(); + let (a2, x, y) = bashein_detmer_sums(ring); + sum_a2 = sum_a2 + a2; + sum_x = sum_x + x; + sum_y = sum_y + y; }; fold_ring(pg.exterior()); @@ -277,10 +286,55 @@ where fold_ring(inner); } - if sum_area == zero { + if sum_a2 == zero { return pg.exterior().points().next().copied().unwrap_or_default(); } - point_2d::(sum_x / sum_area, sum_y / sum_area) + let a3 = three::<::Scalar>() * sum_a2; + point_2d::(sum_x / a3, sum_y / a3) + } +} + +// ---- MultiPolygon ---------------------------------------------------- +// +// Mirrors the multi-polygon arm of `algorithms/centroid.hpp`, which runs +// one `centroid_multi` state over every ring of every member and divides +// once. Same reason the polygon arm accumulates rather than combining +// per-part centroids: a member with zero area drops out of an +// area-weighted combine but still contributes to the running numerator, +// and Boost keeps that contribution. + +impl CentroidStrategy for CartesianMultiPolygonCentroid +where + G: MultiPolygonTrait, + G::Point: PointTrait + PointMut + Default + Copy, + <::Cs as CoordinateSystem>::Family: SameAs, +{ + type Output = G::Point; + + fn centroid(&self, mp: &G) -> G::Point { + let zero = ::Scalar::ZERO; + let mut sum_a2 = zero; + let mut sum_x = zero; + let mut sum_y = zero; + let mut first_point = None; + + for polygon in mp.polygons() { + if first_point.is_none() { + first_point = polygon.exterior().points().next().copied(); + } + for ring in core::iter::once(polygon.exterior()).chain(polygon.interiors()) { + let (a2, x, y) = bashein_detmer_sums(ring); + sum_a2 = sum_a2 + a2; + sum_x = sum_x + x; + sum_y = sum_y + y; + } + } + + if sum_a2 == zero { + return first_point.unwrap_or_default(); + } + let a3 = three::<::Scalar>() * sum_a2; + point_2d::(sum_x / a3, sum_y / a3) } } @@ -459,6 +513,10 @@ impl CentroidStrategyForKind for RingTag { type S = CartesianRingCentroid; } +impl CentroidStrategyForKind for MultiPolygonTag { + type S = CartesianMultiPolygonCentroid; +} + impl CentroidStrategyForKind for PolygonTag { type S = CartesianPolygonCentroid; } @@ -491,11 +549,13 @@ mod tests { use super::{ CartesianBoxCentroid, CartesianLinestringCentroid, CartesianMultiPointCentroid, - CartesianPolygonCentroid, CartesianRingCentroid, CartesianSegmentCentroid, - CentroidStrategy, + CartesianMultiPolygonCentroid, CartesianPolygonCentroid, CartesianRingCentroid, + CartesianSegmentCentroid, CentroidStrategy, }; use geometry_cs::Cartesian; - use geometry_model::{Box, MultiPoint, Point2D, Polygon, Ring, Segment, linestring, polygon}; + use geometry_model::{ + Box, MultiPoint, MultiPolygon, Point2D, Polygon, Ring, Segment, linestring, polygon, + }; use geometry_trait::Point as _; type Pt = Point2D; @@ -664,4 +724,59 @@ mod tests { let c = CartesianMultiPointCentroid.centroid(&mp); assert!(close_pt(&c, 2.0 / 3.0, 2.0 / 3.0, 1e-9)); } + + /// A part with zero area still contributes to the running numerator. + /// + /// Combining per-part centroids weighted by area drops it — its weight is + /// zero — and lands somewhere else. Boost 1.83 on a clockwise + /// `model::polygon` / `model::multi_polygon`: + /// + /// ```text + /// bowtie exterior + hole -> (10.6667, 11) area=-4 + /// zero-area MP -> (0, 0) area=0 + /// mixed MP -> (11.3333, 11) area=4 + /// ``` + #[test] + fn a_zero_area_part_still_moves_the_centroid() { + // Exterior is a bow-tie: zero area, non-zero numerator. + let bowtie_with_hole: Polygon = polygon![ + [(0.0, 0.0), (2.0, 2.0), (2.0, 0.0), (0.0, 2.0), (0.0, 0.0)], + [ + (10.0, 10.0), + (12.0, 10.0), + (12.0, 12.0), + (10.0, 12.0), + (10.0, 10.0) + ] + ]; + let c = CartesianPolygonCentroid.centroid(&bowtie_with_hole); + assert!(close_pt(&c, 32.0 / 3.0, 11.0, 1e-9), "{c:?}"); + + let bowtie: Polygon = + polygon![[(0.0, 0.0), (2.0, 2.0), (2.0, 0.0), (0.0, 2.0), (0.0, 0.0)]]; + let other_bowtie: Polygon = polygon![[ + (10.0, 10.0), + (12.0, 12.0), + (12.0, 10.0), + (10.0, 12.0), + (10.0, 10.0) + ]]; + let square: Polygon = polygon![[ + (10.0, 10.0), + (10.0, 12.0), + (12.0, 12.0), + (12.0, 10.0), + (10.0, 10.0) + ]]; + + // Every member degenerate: the first vertex of the first member. + let all_degenerate = MultiPolygon(vec![bowtie.clone(), other_bowtie]); + let c = CartesianMultiPolygonCentroid.centroid(&all_degenerate); + assert!(close_pt(&c, 0.0, 0.0, 1e-9), "{c:?}"); + + // One degenerate member beside a real one: it still pulls the result. + let mixed = MultiPolygon(vec![bowtie, square]); + let c = CartesianMultiPolygonCentroid.centroid(&mixed); + assert!(close_pt(&c, 34.0 / 3.0, 11.0, 1e-9), "{c:?}"); + } } diff --git a/crates/geometry-strategy/src/lib.rs b/crates/geometry-strategy/src/lib.rs index f48c0a3..49a7ab4 100644 --- a/crates/geometry-strategy/src/lib.rs +++ b/crates/geometry-strategy/src/lib.rs @@ -140,8 +140,8 @@ pub use buffer::{ pub use cartesian::{ComparablePythagoras, PointToSegment, Pythagoras}; pub use centroid::{ CartesianBoxCentroid, CartesianLinestringCentroid, CartesianMultiPointCentroid, - CartesianPolygonCentroid, CartesianRingCentroid, CartesianSegmentCentroid, CentroidStrategy, - CentroidStrategyForKind, + CartesianMultiPolygonCentroid, CartesianPolygonCentroid, CartesianRingCentroid, + CartesianSegmentCentroid, CentroidStrategy, CentroidStrategyForKind, }; pub use closest_points::{CartesianClosestPoints, ClosestPointsStrategy}; pub use compare::{ALL_DIMENSIONS, EqualTo, Greater, Less, LessExact}; From ac9bcea7cee2db1c300e9fe04552e12fb5014dfb Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 13:21:35 -0600 Subject: [PATCH 05/20] fix(overlay): distinguish the two ways multi-polygon members can be wrong Every inter-member defect was reported as `IntersectingInteriors`. Boost reports `failure_self_intersections` for most of them and reserves `failure_intersecting_interiors` for genuine nesting. Boost 1.83: overlapping members failure=21 boundaries cross edge-touching members failure=21 boundaries share a curve identical members failure=21 one inside another failure=40 interiors overlap, boundaries never meet point-touching members valid disjoint members valid So a boundary that meets the other at all is 21, and only interiors that overlap while the boundaries stay apart is 40. Callers branch on the difference: tilemaker's `buildWayGeometry` throws its fast clip away and redoes it with the robust intersection on 21, and leaves the geometry alone on 40. The existing parity test asserted 40 for the overlapping and edge-touching cases; both measure 21. Updated, and a nested case added to cover the 40 arm that is now reachable on its own. --- crates/geometry-overlay/src/validity.rs | 31 ++++++++++++++++--- .../tests/validity_completion_parity.rs | 25 ++++++++++++++- 2 files changed, 51 insertions(+), 5 deletions(-) diff --git a/crates/geometry-overlay/src/validity.rs b/crates/geometry-overlay/src/validity.rs index 818a1ab..01671b3 100644 --- a/crates/geometry-overlay/src/validity.rs +++ b/crates/geometry-overlay/src/validity.rs @@ -396,10 +396,33 @@ where for second in (first + 1)..polygons.len() { let matrix = crate::relate::relate(polygons[first], polygons[second]) .map_err(|_| ValidityFailure::SelfIntersection)?; - if matrix.interior_interior() == crate::relate::Dimension::Area - || matrix.boundary_boundary() == crate::relate::Dimension::Curve - { - return Err(ValidityFailure::IntersectingInteriors); + let interiors_overlap = + matrix.interior_interior() == crate::relate::Dimension::Area; + let boundaries = matrix.boundary_boundary(); + + // Boost distinguishes two ways members can be wrong, and + // tilemaker's `buildWayGeometry` branches on which: + // + // overlapping members failure=21 boundaries cross + // edge-touching members failure=21 boundaries share a curve + // identical members failure=21 + // one inside another failure=40 interiors overlap, boundaries do not meet + // point-touching members valid + // disjoint members valid + // + // So a boundary that meets the other at all — sharing a curve, + // or crossing it while the interiors overlap — is + // `failure_self_intersections`, and only genuine nesting is + // `failure_intersecting_interiors`. + if boundaries == crate::relate::Dimension::Curve { + return Err(ValidityFailure::SelfIntersection); + } + if interiors_overlap { + return Err(if boundaries == crate::relate::Dimension::Empty { + ValidityFailure::IntersectingInteriors + } else { + ValidityFailure::SelfIntersection + }); } } } diff --git a/crates/geometry/tests/validity_completion_parity.rs b/crates/geometry/tests/validity_completion_parity.rs index d5130d8..b8007a1 100644 --- a/crates/geometry/tests/validity_completion_parity.rs +++ b/crates/geometry/tests/validity_completion_parity.rs @@ -127,6 +127,18 @@ fn polygon_detects_disconnected_interior() { /// `test/algorithms/is_valid.cpp:929-970` — multi-polygons may touch at an /// isolated point or place a member in another member's hole, but filled /// interiors may not overlap/share an edge. +/// +/// Which *failure* they report is not the same for every way of being wrong, +/// and callers branch on it. Boost 1.83: +/// +/// ```text +/// overlapping members failure=21 boundaries cross +/// edge-touching members failure=21 boundaries share a curve +/// identical members failure=21 +/// one inside another failure=40 interiors overlap, boundaries never meet +/// point-touching members valid +/// disjoint members valid +/// ``` #[test] fn multipolygon_checks_inter_member_topology() { let invalid_member: MultiPolygon> = MultiPolygon::from_vec(vec![Polygon::new( @@ -138,13 +150,24 @@ fn multipolygon_checks_inter_member_topology() { MultiPolygon::from_vec(vec![square(0.0, 0.0, 4.0, 4.0), square(2.0, 2.0, 6.0, 6.0)]); assert_eq!( is_valid(&overlapping), - Err(ValidityFailure::IntersectingInteriors) + Err(ValidityFailure::SelfIntersection) ); let shared_edge = MultiPolygon::from_vec(vec![square(0.0, 0.0, 2.0, 2.0), square(2.0, 0.0, 4.0, 2.0)]); assert_eq!( is_valid(&shared_edge), + Err(ValidityFailure::SelfIntersection) + ); + + // Only genuine nesting reaches `IntersectingInteriors`: the interiors + // overlap and the boundaries never meet. + let nested = MultiPolygon::from_vec(vec![ + square(0.0, 0.0, 10.0, 10.0), + square(2.0, 2.0, 4.0, 4.0), + ]); + assert_eq!( + is_valid(&nested), Err(ValidityFailure::IntersectingInteriors) ); From 764b1544467946b35da2db3c677177ec9eb947eb Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 13:25:03 -0600 Subject: [PATCH 06/20] feat(overlay): multi-polygon operands for the Boolean operations MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `intersection`, `union`, `difference` and `sym_difference` took a pair of polygons only, so a caller holding multi-polygons had to decompose them, run the pairs, and re-combine — which is not the same function, and quietly changes the answer when the pieces meet. Boost dispatches every areal Boolean through one overlay whatever the operand arity: `bg::intersection(mp, box, out)` and `bg::difference(mp1, mp2, out)` are the same algorithm as the polygon pair. So the arrangement kernel is now shared, and `*_multi` runs it over the union of both operands' rings. A single polygon is the one-member case. No behaviour change for the existing entry points: `overlay` builds the same `Shape` and segment list it always did, and the split-edge core is byte-identical to what it was. Covered by `multi_polygon_operands`, where the operands overlap across members: two disjoint unit squares against one that straddles the first. --- crates/geometry-overlay/src/lib.rs | 5 +- crates/geometry-overlay/src/operation.rs | 5 +- .../geometry-overlay/src/operation/areal.rs | 88 +++++++++++++++++-- .../geometry-overlay/src/operation/boolean.rs | 86 +++++++++++++++++- .../geometry-overlay/tests/overlay_parity.rs | 33 ++++++- 5 files changed, 206 insertions(+), 11 deletions(-) diff --git a/crates/geometry-overlay/src/lib.rs b/crates/geometry-overlay/src/lib.rs index aa11aca..c99d8ae 100644 --- a/crates/geometry-overlay/src/lib.rs +++ b/crates/geometry-overlay/src/lib.rs @@ -54,7 +54,10 @@ pub use line_intersection::{LineIntersection, line_intersection}; // feature-group: Mutation & assembly pub use merge::{merge_elements, merge_multipolygon, merge_polygons, stitch_triangles}; // feature-group: Boolean operations -pub use operation::{OverlayError, difference, intersection, sym_difference, r#union, union_poly}; +pub use operation::{ + OverlayError, difference, difference_multi, intersection, intersection_multi, sym_difference, + sym_difference_multi, r#union, union_multi, union_poly, +}; // feature-group: Spatial predicates pub use relate::{ De9im, Dimension, RelateError, contains_properly, crosses, overlaps, relate as relate_matrix, diff --git a/crates/geometry-overlay/src/operation.rs b/crates/geometry-overlay/src/operation.rs index 3710ee8..39661aa 100644 --- a/crates/geometry-overlay/src/operation.rs +++ b/crates/geometry-overlay/src/operation.rs @@ -8,4 +8,7 @@ mod areal; mod boolean; -pub use boolean::{OverlayError, difference, intersection, sym_difference, r#union, union_poly}; +pub use boolean::{ + OverlayError, difference, difference_multi, intersection, intersection_multi, sym_difference, + sym_difference_multi, r#union, union_multi, union_poly, +}; diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 89887f9..e73a9ec 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -20,7 +20,9 @@ use geometry_coords::{ use geometry_cs::{CartesianFamily, CoordinateSystem}; use geometry_model::{MultiPolygon, Polygon, Ring, Segment}; use geometry_tag::SameAs; -use geometry_trait::{Point, PointMut, Polygon as PolygonTrait, Ring as RingTrait}; +use geometry_trait::{ + MultiPolygon as MultiPolygonTrait, Point, PointMut, Polygon as PolygonTrait, Ring as RingTrait, +}; use crate::assemble::assemble_multipolygon; use crate::operation::OverlayError; @@ -82,6 +84,20 @@ impl Shape { Self { rings } } + fn from_multi_polygon(multi_polygon: &G) -> Self + where + G: MultiPolygonTrait, + P: Point, + P::Scalar: Into, + { + let mut rings = Vec::new(); + for polygon in multi_polygon.polygons() { + rings.push(ring_coordinates(polygon.exterior())); + rings.extend(polygon.interiors().map(ring_coordinates)); + } + Self { rings } + } + fn contains(&self, point: Coordinate) -> bool { self.rings .iter() @@ -146,14 +162,58 @@ where P::Scalar: CoordinateScalar + Into, ::Family: SameAs, { - let first_shape = Shape::from_polygon(first); - let second_shape = Shape::from_polygon(second); - let scale = coordinate_scale(&first_shape, &second_shape); + overlay_arrangement( + &Shape::from_polygon(first), + &Shape::from_polygon(second), + polygon_segments(first), + polygon_segments(second), + operation, + ) +} + +/// The same operation over multi-polygons. +/// +/// Boost dispatches every areal Boolean through one overlay regardless of how +/// many polygons each operand holds, so this is the same kernel over the union +/// of every operand's rings rather than a second algorithm. A single polygon +/// is the one-member case. +pub(crate) fn overlay_multi( + first: &G1, + second: &G2, + operation: ArealOp, +) -> Result>, OverlayError> +where + G1: MultiPolygonTrait, + G2: MultiPolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + overlay_arrangement( + &Shape::from_multi_polygon(first), + &Shape::from_multi_polygon(second), + multi_polygon_segments(first), + multi_polygon_segments(second), + operation, + ) +} + +fn overlay_arrangement

( + first_shape: &Shape, + second_shape: &Shape, + mut first_segments: Vec>, + mut second_segments: Vec>, + operation: ArealOp, +) -> Result>, OverlayError> +where + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + let scale = coordinate_scale(first_shape, second_shape); let snap_tolerance = scale * 1e-10; let parameter_tolerance = 1e-10; - let mut first_segments = polygon_segments(first); - let mut second_segments = polygon_segments(second); for first_segment in &mut first_segments { for second_segment in &mut second_segments { let first_model = Segment::new(first_segment.start, first_segment.end); @@ -251,6 +311,22 @@ where coordinates } +fn multi_polygon_segments(multi_polygon: &G) -> Vec> +where + G: MultiPolygonTrait, + P: Point + Copy, + P::Scalar: Into, +{ + let mut segments = Vec::new(); + for polygon in multi_polygon.polygons() { + append_ring_segments(polygon.exterior(), &mut segments); + for ring in polygon.interiors() { + append_ring_segments(ring, &mut segments); + } + } + segments +} + fn polygon_segments(polygon: &G) -> Vec> where G: PolygonTrait, diff --git a/crates/geometry-overlay/src/operation/boolean.rs b/crates/geometry-overlay/src/operation/boolean.rs index 2969cfb..346e182 100644 --- a/crates/geometry-overlay/src/operation/boolean.rs +++ b/crates/geometry-overlay/src/operation/boolean.rs @@ -25,11 +25,11 @@ use geometry_coords::CoordinateScalar; use geometry_cs::{CartesianFamily, CoordinateSystem}; use geometry_model::{MultiPolygon, Polygon}; use geometry_tag::SameAs; -use geometry_trait::{PointMut, Polygon as PolygonTrait}; +use geometry_trait::{MultiPolygon as MultiPolygonTrait, PointMut, Polygon as PolygonTrait}; use crate::traverse::TraversalError; -use super::areal::{ArealOp, overlay as areal_overlay}; +use super::areal::{ArealOp, overlay as areal_overlay, overlay_multi as areal_overlay_multi}; /// Failure of a boolean overlay operation. #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -219,6 +219,88 @@ where areal_overlay(g1, g2, ArealOp::SymDifference) } +// ---- multi-polygon operands ------------------------------------------ +// +// Boost dispatches every areal Boolean through one overlay whatever the +// operand arity — `bg::intersection(mp, box, out)` and +// `bg::difference(mp1, mp2, out)` are the same algorithm as the polygon +// pair. These are that same kernel with both operands' rings, so a caller +// holding multi-polygons does not have to decompose them and re-combine +// the pieces itself, which is not the same function. + +/// Intersection of two multi-polygons — the region inside **both**. +/// +/// # Errors +/// +/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range. +pub fn intersection_multi( + g1: &G1, + g2: &G2, +) -> Result>, OverlayError> +where + G1: MultiPolygonTrait, + G2: MultiPolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + areal_overlay_multi(g1, g2, ArealOp::Intersection) +} + +/// Union of two multi-polygons — the region inside **either**. +/// +/// # Errors +/// +/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range. +pub fn union_multi(g1: &G1, g2: &G2) -> Result>, OverlayError> +where + G1: MultiPolygonTrait, + G2: MultiPolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + areal_overlay_multi(g1, g2, ArealOp::Union) +} + +/// Difference of two multi-polygons — the region inside `g1` but not `g2`. +/// +/// # Errors +/// +/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range. +pub fn difference_multi( + g1: &G1, + g2: &G2, +) -> Result>, OverlayError> +where + G1: MultiPolygonTrait, + G2: MultiPolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + areal_overlay_multi(g1, g2, ArealOp::Difference) +} + +/// Symmetric difference of two multi-polygons — inside exactly one of them. +/// +/// # Errors +/// +/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range. +pub fn sym_difference_multi( + g1: &G1, + g2: &G2, +) -> Result>, OverlayError> +where + G1: MultiPolygonTrait, + G2: MultiPolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + areal_overlay_multi(g1, g2, ArealOp::SymDifference) +} + #[cfg(test)] mod tests { use super::{OverlayError, intersection, union_poly}; diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 6bb928d..bb67979 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -15,7 +15,10 @@ use geometry_algorithm::ring_area; use geometry_cs::Cartesian; use geometry_model::{MultiPolygon, Point2D, Polygon, polygon}; -use geometry_overlay::{difference, intersection, sym_difference, union_poly}; +use geometry_overlay::{ + difference, difference_multi, intersection, intersection_multi, sym_difference, union_multi, + union_poly, +}; use geometry_trait::{MultiPolygon as _, Polygon as _}; type P = Point2D; @@ -213,3 +216,31 @@ fn intersection_splits_lobes_meeting_at_a_point() { sizes.sort_unstable(); assert_eq!(sizes, [4, 7], "each lobe keeps its own closed ring"); } + +// ---- Multi-polygon operands ------------------------------------------ + +/// The multi-polygon entry points are the same overlay over both operands' +/// rings, not a decomposition into per-member pairs. Two disjoint unit +/// squares against a third that overlaps one of them: +/// +/// ```text +/// A = {(0,0)-(1,1), (4,0)-(5,1)} area 2 +/// B = {(0.5,0)-(1.5,1)} area 1 +/// A ∪ B area 2.5 A ∩ B area 0.5 A − B area 1.5 +/// ``` +#[test] +fn multi_polygon_operands() { + let a: MultiPolygon> = + MultiPolygon::from_vec(vec![square(0.0, 0.0, 1.0), square(4.0, 0.0, 1.0)]); + let b: MultiPolygon> = MultiPolygon::from_vec(vec![polygon![[ + (0.5, 0.0), + (0.5, 1.0), + (1.5, 1.0), + (1.5, 0.0), + (0.5, 0.0) + ]]]); + + close(area(&union_multi(&a, &b).unwrap()), 2.5); + close(area(&intersection_multi(&a, &b).unwrap()), 0.5); + close(area(&difference_multi(&a, &b).unwrap()), 1.5); +} From 01cfdde1f3948f2cccd1f538ecc4eee0ad541411 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 13:42:40 -0600 Subject: [PATCH 07/20] fix(overlay): a hole sharing an edge with the exterior is a self-intersection MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `is_valid_polygon_with` reported `DisconnectedInterior` for any hole that overlapped the exterior boundary. Boost reports `failure_self_intersections` when they share a *curve*, and reserves `failure_disconnected_interior` for isolated contacts that cut the interior in two. Boost 1.83, exterior (0,0)-(10,10) clockwise: hole fully interior valid hole touching the exterior at one point valid hole sharing a segment of it failure=21 hole touching it at two isolated points failure=32 The inner-vs-inner loop a few lines below already made that distinction; the exterior-vs-inner one did not. Found in tilemaker: a building clipped to a tile edge acquires a hole that runs along the clip boundary, so it hits this on real data. The two codes take different branches there — `buildWayGeometry` redoes the clip with the robust intersection on 21 and only repairs in place on 32 — and the tile came out with two vertices the reference does not have. The existing `polygon_detects_disconnected_interior` case asserted 32 for a hole sharing the left edge; it measures 21. Replaced with the four shapes above, so both arms are covered by a case that reaches them on its own. --- crates/geometry-overlay/src/validity.rs | 14 +++++- .../tests/validity_completion_parity.rs | 50 +++++++++++++++++-- 2 files changed, 57 insertions(+), 7 deletions(-) diff --git a/crates/geometry-overlay/src/validity.rs b/crates/geometry-overlay/src/validity.rs index 01671b3..35d8a71 100644 --- a/crates/geometry-overlay/src/validity.rs +++ b/crates/geometry-overlay/src/validity.rs @@ -665,10 +665,20 @@ where return Err(ValidityFailure::InteriorRingOutside); } let interaction = ring_pair_interaction(polygon.exterior(), *inner); - if interaction.proper_crossing { + // A hole that shares a *curve* with the exterior is a + // self-intersection, the same way two holes sharing one is (below). + // Only isolated contacts that cut the interior in two are + // `failure_disconnected_interior`. Boost 1.83, exterior + // (0,0)-(10,10) clockwise: + // + // hole fully interior valid + // hole touching the exterior at one point valid + // hole touching it at two points failure=32 + // hole sharing a segment of it failure=21 + if interaction.proper_crossing || interaction.overlap { return Err(ValidityFailure::SelfIntersection); } - if interaction.overlap || interaction.contacts.len() > 1 { + if interaction.contacts.len() > 1 { return Err(ValidityFailure::DisconnectedInterior); } } diff --git a/crates/geometry/tests/validity_completion_parity.rs b/crates/geometry/tests/validity_completion_parity.rs index b8007a1..c909635 100644 --- a/crates/geometry/tests/validity_completion_parity.rs +++ b/crates/geometry/tests/validity_completion_parity.rs @@ -61,12 +61,22 @@ fn polygon_rejects_crossing_and_nested_interior_rings() { ); } -/// `test/algorithms/is_valid.cpp:612-628,680-687` — an exterior-edge contact -/// can disconnect the polygon interior while isolated boundary points remain -/// admissible. +/// How a hole may and may not meet the exterior, and which failure each way +/// of being wrong reports. Boost 1.83, exterior `(0,0)-(10,10)` clockwise: +/// +/// ```text +/// hole fully interior valid +/// hole touching the exterior at one point valid +/// hole sharing a segment of it failure=21 +/// hole touching it at two isolated points failure=32 +/// ``` +/// +/// A shared *curve* is a self-intersection, the same as two holes sharing +/// one. Only isolated contacts that cut the interior in two are +/// `failure_disconnected_interior`. #[test] fn polygon_detects_disconnected_interior() { - let edge_touch: Polygon

= polygon![ + let shared_edge: Polygon

= polygon![ [ (0.0, 0.0), (0.0, 10.0), @@ -77,10 +87,40 @@ fn polygon_detects_disconnected_interior() { [(0.0, 3.0), (3.0, 3.0), (3.0, 7.0), (0.0, 7.0), (0.0, 3.0)] ]; assert_eq!( - is_valid(&edge_touch), + is_valid(&shared_edge), + Err(ValidityFailure::SelfIntersection) + ); + + // A hole reaching from the top edge to the bottom one, meeting each at a + // single point: the interior really is cut in two. + let pinched_in_two: Polygon

= polygon![ + [ + (0.0, 0.0), + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 0.0) + ], + [(5.0, 10.0), (2.0, 5.0), (5.0, 0.0), (8.0, 5.0), (5.0, 10.0)] + ]; + assert_eq!( + is_valid(&pinched_in_two), Err(ValidityFailure::DisconnectedInterior) ); + // One isolated contact leaves the interior connected. + let touches_once: Polygon

= polygon![ + [ + (0.0, 0.0), + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 0.0) + ], + [(5.0, 10.0), (2.0, 5.0), (8.0, 5.0), (5.0, 10.0)] + ]; + assert!(is_valid(&touches_once).is_ok()); + let holes_share_edge: Polygon

= polygon![ [ (0.0, 0.0), From 46d074d308d8bd9f997e0e4e7a7cf6455938757f Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 13:59:01 -0600 Subject: [PATCH 08/20] fix(predicate): compute a segment crossing parametrically, not by determinant MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `line_cross_point` solved the two-line determinant, whose terms are products of *absolute* coordinates. On geographic input that is the wrong arithmetic: a polygon spanning 1e-4 degrees at longitude 7.4, latitude 48.7 forms terms around 362 and asks for a residue ten orders of magnitude smaller, so the crossing lands off the line it is supposed to be on. segment (7.4264305,48.7356554)-(7.4263938,48.7355668) against the horizontal clip edge y = 48.7355713 determinant y off by 9.24e-14, and 4.97e-14 off the sloped segment parametric y exact, 3.09e-20 off the sloped segment Split edges that should meet at a shared node then do not, the arrangement is left unclosed, and the traversal strands: `intersection` returned `OverlayError::Unsupported` for a building clipped to a tile edge. The parametric form touches only coordinate *differences* — the size of the geometry rather than of its position — and anchors the answer on the first segment's start, so it stays on that segment by construction. Same point in exact arithmetic, and the same cost. Measured on 17 polygon-against-clip-box cases taken from real tilemaker runs, where the expected output is C++ Boost 1.83's: 14 differed and 2 refused before, all 17 exact after. --- .../src/predicate/segment_intersection.rs | 71 +++++++++++++++++-- 1 file changed, 66 insertions(+), 5 deletions(-) diff --git a/crates/geometry-overlay/src/predicate/segment_intersection.rs b/crates/geometry-overlay/src/predicate/segment_intersection.rs index 2fbe83a..19fdca3 100644 --- a/crates/geometry-overlay/src/predicate/segment_intersection.rs +++ b/crates/geometry-overlay/src/predicate/segment_intersection.rs @@ -201,12 +201,25 @@ where let x4 = p4.get::<0>(); let y4 = p4.get::<1>(); - // Standard two-line determinant solution. + // The parametric solution, not the two-line determinant one. + // + // Both are the same point in exact arithmetic. The determinant form builds + // `x1*y2 - y1*x2`, a product of *absolute* coordinates, and then subtracts + // two such products that are nearly equal; the answer it is looking for is + // the small residue. On geographic input that residue is the whole result: + // a polygon spanning 1e-4 degrees at longitude 7.4, latitude 48.7 forms + // terms around 362 and asks for a difference ten orders of magnitude below + // them, so the crossing lands off the line it is supposed to be on. Split + // edges then fail to meet at a shared node and the arrangement cannot be + // traced. + // + // The parametric form touches only coordinate *differences*, which are the + // size of the geometry rather than of its position, and anchors the result + // on `p1` so it stays on segment `a`. let denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4); - let a = x1 * y2 - y1 * x2; - let b = x3 * y4 - y3 * x4; - let px = (a * (x3 - x4) - (x1 - x2) * b) / denom; - let py = (a * (y3 - y4) - (y1 - y2) * b) / denom; + let t = ((x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)) / denom; + let px = x1 + t * (x2 - x1); + let py = y1 + t * (y2 - y1); make_point::

(px, py) } @@ -298,6 +311,7 @@ mod tests { use super::{SegmentIntersection, segment_intersection}; use geometry_cs::Cartesian; use geometry_model::{Point2D, Segment}; + use geometry_trait::Point as _; type P = Point2D; type Seg = Segment

; @@ -398,4 +412,51 @@ mod tests { SegmentIntersection::Single(P::new(2.0, 2.0)) ); } + + /// The crossing has to lie *on* both segments, at the magnitudes + /// geographic data actually uses. + /// + /// This one is taken from a tilemaker run: a building clipped to a tile + /// edge at longitude 7.4, latitude 48.7, spanning about 1e-4 degrees. The + /// clip edge is horizontal, so the answer's y must equal the edge's y + /// exactly. The determinant form misses it by 9.24e-14 — its terms are + /// products of absolute coordinates around 362 and the residue it wants is + /// ten orders of magnitude smaller. Split edges then stop meeting at a + /// shared node and the arrangement cannot be traced. + #[test] + fn crossing_lies_on_the_segments_at_geographic_magnitude() { + let edge_y = 48.735_571_3; + let sloped = Segment::new( + P::new(7.426_430_5, 48.735_655_4), + P::new(7.426_393_8, 48.735_566_8), + ); + let clip_edge = Segment::new(P::new(7.382_592_8, edge_y), P::new(7.426_977_5, edge_y)); + + let SegmentIntersection::Single(crossing) = segment_intersection(&sloped, &clip_edge) + else { + panic!("the segments cross"); + }; + // Exact equality on purpose: the edge is horizontal, so the crossing's + // y is one of the inputs and no arithmetic should alter it. + #[expect( + clippy::float_cmp, + reason = "the crossing's y must reproduce the horizontal edge's y bit for bit" + )] + { + assert_eq!( + crossing.get::<1>(), + edge_y, + "the crossing must sit exactly on the horizontal edge" + ); + } + + // And on the sloped segment, to the last bits its own span allows. + let (x1, y1) = (7.426_430_5_f64, 48.735_655_4_f64); + let (x2, y2) = (7.426_393_8_f64, 48.735_566_8_f64); + let side = (x2 - x1) * (crossing.get::<1>() - y1) - (y2 - y1) * (crossing.get::<0>() - x1); + assert!( + side.abs() < 1e-19, + "the crossing must sit on the sloped segment, off by {side:e}" + ); + } } From cee26e9cc4f27e2d7050caacae1978efa79443a1 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 14:05:07 -0600 Subject: [PATCH 09/20] fix(overlay): emit rings, and their vertices, in Boost's order MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Which vertex a ring starts at and which polygon comes first are both observable, and both were falling out of whichever boundary edge the traversal happened to seed from. Boost's are not arbitrary: * a ring begins at a **turn**, the first one in source order; * a ring with no turn was copied whole from one operand and keeps that operand's own starting vertex — not the end the traversal reached it from, which for a hole is the other one, since a hole is walked against its stored direction; * the rings come out in source order, so the polygons do too. Nodes are numbered as the operands are walked, so "first in source order" is just the smallest node index; the tracer now records which nodes are turns (`push_split` marks them, including a crossing that lands on a vertex) and rotates each cycle onto its earliest one. This matters to any consumer that pins a ring's first vertex. tilemaker's `simplify_ring` never removes one, so a rotation of one place shows up in the tile as one extra vertex and one missing. Measured on 27 polygon-against-clip-box cases from real tilemaker runs, compared to C++ Boost 1.83 including ring start and polygon order: 0 of 27 matched before, 25 of 27 after. The two that remain are one case whose two output rings begin at the same vertex, where source order does not separate them. --- .../geometry-overlay/src/operation/areal.rs | 112 ++++++++++++++---- 1 file changed, 88 insertions(+), 24 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index e73a9ec..94deeeb 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -108,7 +108,10 @@ impl Shape { struct SourceSegment

{ start: P, end: P, - splits: Vec<(f64, P)>, + /// `(parameter, point, is_turn)`. The two endpoints are not turns; a point + /// pushed by the intersection sweep is, including one that lands on an + /// endpoint. + splits: Vec<(f64, P, bool)>, } impl

SourceSegment

@@ -120,27 +123,36 @@ where Self { start, end, - splits: alloc::vec![(0.0, start), (1.0, end)], + splits: alloc::vec![(0.0, start, false), (1.0, end, false)], } } fn push_split(&mut self, point: P, tolerance: f64) { let parameter = segment_parameter(&self.start, &self.end, &point); - if parameter >= -tolerance - && parameter <= 1.0 + tolerance - && !self - .splits - .iter() - .any(|(existing, _)| (existing - parameter).abs() <= tolerance) + if parameter < -tolerance || parameter > 1.0 + tolerance { + return; + } + if let Some(existing) = self + .splits + .iter_mut() + .find(|(at, _, _)| (at - parameter).abs() <= tolerance) { - self.splits.push((parameter.clamp(0.0, 1.0), point)); + // A crossing that lands on a vertex already split here still makes + // that vertex a turn. + existing.2 = true; + return; } + self.splits.push((parameter.clamp(0.0, 1.0), point, true)); } } struct Node

{ point: P, coordinate: Coordinate, + /// Set when any split that resolved to this node was a crossing. Boost + /// starts each output ring at a turn, so the tracer needs to know which + /// nodes are turns; see `push_ring`. + is_turn: bool, } #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -379,8 +391,8 @@ fn append_atomic_edges

( .sort_by(|left, right| left.0.total_cmp(&right.0)); for pair in segment.splits.windows(2) { debug_assert!((pair[1].0 - pair[0].0).abs() > 1e-12); - let start = canonical_node(nodes, pair[0].1, tolerance); - let end = canonical_node(nodes, pair[1].1, tolerance); + let start = canonical_node(nodes, pair[0].1, pair[0].2, tolerance); + let end = canonical_node(nodes, pair[1].1, pair[1].2, tolerance); if start != end { output.push(Edge { start, end }); } @@ -388,7 +400,7 @@ fn append_atomic_edges

( } } -fn canonical_node

(nodes: &mut Vec>, point: P, tolerance: f64) -> usize +fn canonical_node

(nodes: &mut Vec>, point: P, is_turn: bool, tolerance: f64) -> usize where P: Point + Copy, P::Scalar: Into, @@ -400,9 +412,14 @@ where node.coordinate.y - coordinate.y, ) <= tolerance }) { + nodes[index].is_turn |= is_turn; return index; } - nodes.push(Node { point, coordinate }); + nodes.push(Node { + point, + coordinate, + is_turn, + }); nodes.len() - 1 } @@ -416,7 +433,9 @@ where P::Scalar: Into, { let mut used = alloc::vec![false; edges.len()]; - let mut rings = Vec::new(); + // Each ring is kept with the node it starts at, so the whole set can be + // put back into source order below. + let mut rings: Vec<(usize, Ring

)> = Vec::new(); for seed in 0..edges.len() { if used[seed] { continue; @@ -453,12 +472,30 @@ where } debug_assert_eq!(node_indices.last().copied(), Some(first)); } - Ok(rings) + + // Which ring comes out first is observable — it decides the order of the + // polygons in the result — and Boost's is not the order the seeds happened + // to fall in. It traces from turns in source order, so ordering the rings + // by the node each one starts at reproduces it. + rings.sort_by_key(|&(start, _)| start); + Ok(rings.into_iter().map(|(_, ring)| ring).collect()) } /// Append one traced cycle, dropping it when it encloses no area. -fn push_ring

(rings: &mut Vec>, nodes: &[Node

], node_indices: &[usize], tolerance: f64) -where +/// +/// The cycle starts wherever the traversal happened to seed, which carries no +/// meaning — but which vertex a ring starts at is observable downstream, and +/// Boost's answer is not arbitrary: it begins each output ring at a *turn*, the +/// first one in source order. A ring with no turn at all was copied whole from +/// one operand and keeps that operand's own starting vertex. Reproduced here, +/// because a consumer that simplifies the ring afterwards will pin its first +/// vertex and the choice reaches the output. +fn push_ring

( + rings: &mut Vec<(usize, Ring

)>, + nodes: &[Node

], + node_indices: &[usize], + tolerance: f64, +) where P: Point + Copy, P::Scalar: Into, { @@ -467,14 +504,38 @@ where let b = nodes[pair[1]].coordinate; sum + a.x * b.y - b.x * a.y }) * 0.5; - if area.abs() > tolerance * tolerance { - rings.push(Ring::from_vec( - node_indices - .iter() - .map(|&index| nodes[index].point) - .collect(), - )); + if area.abs() <= tolerance * tolerance { + return; + } + + // The cycle is closed, so its last index repeats its first. + let cycle = &node_indices[..node_indices.len() - 1]; + // Nodes are numbered in the order the operands were walked, so the + // smallest index on the ring is the earliest in source order. + let start_at = |turns_only: bool| { + cycle + .iter() + .copied() + .enumerate() + .filter(|&(_, index)| !turns_only || nodes[index].is_turn) + .min_by_key(|&(_, index)| index) + .map(|(position, _)| position) + }; + // A ring with no turn was copied whole from one operand, and keeps that + // operand's own starting vertex rather than whichever end of it the + // traversal happened to seed from — a hole is walked against its stored + // direction, so those differ. + let first_turn = start_at(true).or_else(|| start_at(false)).unwrap_or(0); + + let mut points: Vec

= cycle[first_turn..] + .iter() + .chain(&cycle[..first_turn]) + .map(|&index| nodes[index].point) + .collect(); + if let Some(&first) = points.first() { + points.push(first); } + rings.push((cycle[first_turn], Ring::from_vec(points))); } fn next_edge

(nodes: &[Node

], edges: &[Edge], used: &[bool], incoming: Edge) -> Option @@ -576,14 +637,17 @@ mod tests { Node { point: P::new(0.0, 0.0), coordinate: Coordinate { x: 0.0, y: 0.0 }, + is_turn: false, }, Node { point: P::new(1.0, 0.0), coordinate: Coordinate { x: 1.0, y: 0.0 }, + is_turn: false, }, Node { point: P::new(2.0, 0.0), coordinate: Coordinate { x: 2.0, y: 0.0 }, + is_turn: false, }, ]; let edges = [ From e24204bcb8be434b7e119ae47e946038b3af2597 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 16:58:40 -0600 Subject: [PATCH 10/20] fix(overlay): attach a turn to the segment it terminates MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A union ring begins at a turn — the first one along the first operand's boundary — which `push_ring` decided by node index, the order the nodes happened to be created in. That order is wrong at exactly one place: a turn landing on a vertex is created when the *leaving* segment is split, so a turn on the first operand's own first vertex came out as position zero. Boost normalises the other way. `get_turns` attaches an intersection at a segment endpoint to the segment it **terminates**, so that turn is the *last* position on its ring, not the first. Nodes now carry an `arrival` order built with that normalisation — only the far end of a split counts — and the start turn is chosen by it. Measured against C++ Boost 1.83. Two quadrilaterals sharing an edge, the first given starting at one of that edge's endpoints: before D C B A F E D after A F E D C B A = Boost This is the tenth defect found by diffing the tilemaker port against the C++ it was ported from, and it is what made that port's building outlines differ: `make_valid`'s dissolve merges a self-intersecting outline's two lobes with `union_`, and the merged hole came back rotated three vertices. With this, two of the port's three remaining building differences are gone. Unions whose operands share a *collinear edge* are still not exact — Boost drops the starting turn there when the outline runs straight through it, which is its collinear turn classification and a separate piece of work. 104 of 112 such pairs differ; of the 200 pairs that meet only at crossings or points, 169 are exact and 25 of the rest differ only in the last bits of an irrational crossing. Gate: 1606 tests pass, fmt and clippy clean. --- .../geometry-overlay/src/operation/areal.rs | 54 ++++++++++++--- .../geometry-overlay/tests/overlay_parity.rs | 67 ++++++++++++++++++- 2 files changed, 111 insertions(+), 10 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 94deeeb..54fac89 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -153,6 +153,17 @@ struct Node

{ /// starts each output ring at a turn, so the tracer needs to know which /// nodes are turns; see `push_ring`. is_turn: bool, + /// Where this node sits along the operands' boundaries, first operand + /// first — and, where it lands exactly on a vertex, counted as the *end* + /// of the segment arriving there rather than the start of the one leaving. + /// + /// That normalisation is Boost's, and it is the whole difference between + /// the two: `get_turns` attaches an intersection at a segment endpoint to + /// the segment it terminates, so a turn on the first operand's *first* + /// vertex is the last position on that ring, not the first. Ordering by + /// creation instead put it first, and every union of two rings sharing an + /// edge came out started at the wrong end. See `push_ring`. + arrival: usize, } #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -249,16 +260,21 @@ where let mut nodes = Vec::new(); let mut candidates = Vec::new(); + // One counter across both operands, so every position on the first + // operand's boundary precedes every position on the second's. + let mut arrival = 0; append_atomic_edges( &mut first_segments, &mut nodes, &mut candidates, + &mut arrival, snap_tolerance, ); append_atomic_edges( &mut second_segments, &mut nodes, &mut candidates, + &mut arrival, snap_tolerance, ); @@ -380,6 +396,7 @@ fn append_atomic_edges

( segments: &mut [SourceSegment

], nodes: &mut Vec>, output: &mut Vec, + arrival: &mut usize, tolerance: f64, ) where P: Point + Copy, @@ -393,6 +410,13 @@ fn append_atomic_edges

( debug_assert!((pair[1].0 - pair[0].0).abs() > 1e-12); let start = canonical_node(nodes, pair[0].1, pair[0].2, tolerance); let end = canonical_node(nodes, pair[1].1, pair[1].2, tolerance); + // Only the far end of a split counts as an arrival, which is what + // pushes a ring's first vertex to the end of its own ring: it is + // reached as the last segment's endpoint, not the first's start. + if nodes[end].arrival == usize::MAX { + nodes[end].arrival = *arrival; + *arrival += 1; + } if start != end { output.push(Edge { start, end }); } @@ -419,6 +443,7 @@ where point, coordinate, is_turn, + arrival: usize::MAX, }); nodes.len() - 1 } @@ -510,22 +535,30 @@ fn push_ring

( // The cycle is closed, so its last index repeats its first. let cycle = &node_indices[..node_indices.len() - 1]; - // Nodes are numbered in the order the operands were walked, so the - // smallest index on the ring is the earliest in source order. - let start_at = |turns_only: bool| { + // Boost begins each output ring at the first *turn* along the first + // operand's boundary — `Node::arrival`, which is that order with Boost's + // endpoint normalisation applied. + let first_turn_by_arrival = cycle + .iter() + .copied() + .enumerate() + .filter(|&(_, index)| nodes[index].is_turn) + .min_by_key(|&(_, index)| nodes[index].arrival) + .map(|(position, _)| position); + // A ring with no turn was copied whole from one operand, and keeps that + // operand's own starting vertex rather than whichever end of it the + // traversal happened to seed from — a hole is walked against its stored + // direction, so those differ. That vertex is the one created first, which + // is node order, not arrival order. + let first_node = || { cycle .iter() .copied() .enumerate() - .filter(|&(_, index)| !turns_only || nodes[index].is_turn) .min_by_key(|&(_, index)| index) .map(|(position, _)| position) }; - // A ring with no turn was copied whole from one operand, and keeps that - // operand's own starting vertex rather than whichever end of it the - // traversal happened to seed from — a hole is walked against its stored - // direction, so those differ. - let first_turn = start_at(true).or_else(|| start_at(false)).unwrap_or(0); + let first_turn = first_turn_by_arrival.or_else(first_node).unwrap_or(0); let mut points: Vec

= cycle[first_turn..] .iter() @@ -638,16 +671,19 @@ mod tests { point: P::new(0.0, 0.0), coordinate: Coordinate { x: 0.0, y: 0.0 }, is_turn: false, + arrival: 0, }, Node { point: P::new(1.0, 0.0), coordinate: Coordinate { x: 1.0, y: 0.0 }, is_turn: false, + arrival: 1, }, Node { point: P::new(2.0, 0.0), coordinate: Coordinate { x: 2.0, y: 0.0 }, is_turn: false, + arrival: 2, }, ]; let edges = [ diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index bb67979..4540f04 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -19,7 +19,7 @@ use geometry_overlay::{ difference, difference_multi, intersection, intersection_multi, sym_difference, union_multi, union_poly, }; -use geometry_trait::{MultiPolygon as _, Polygon as _}; +use geometry_trait::{MultiPolygon as _, Point as _, Polygon as _, Ring as _}; type P = Point2D; @@ -244,3 +244,68 @@ fn multi_polygon_operands() { close(area(&intersection_multi(&a, &b).unwrap()), 0.5); close(area(&difference_multi(&a, &b).unwrap()), 1.5); } + +// ---- Where a union ring starts, when the first operand starts at a turn ---- +// +// Boost begins each output ring at a turn — the first one along the *first* +// operand's boundary. Which turn that is depends on a normalisation in +// `get_turns`: an intersection landing exactly on a vertex is attached to the +// segment it **terminates**, not the one it begins. So a turn on the first +// operand's own first vertex is the *last* position on that ring, not the +// first. +// +// Reference values from C++ Boost 1.83 on the same input, through +// `scripts/geometry-ab/cpp_ops.cpp` in the tilemaker port. + +fn vertices(mp: &MultiPolygon>) -> Vec<(f64, f64)> { + mp.polygons() + .next() + .expect("one polygon") + .exterior() + .points() + .map(|p| (p.get::<0>(), p.get::<1>())) + .collect() +} + +/// A square and a triangle sharing the square's bottom edge. Both ends of that +/// edge are corners of the union, so nothing is dropped and the only question +/// is which one the ring starts at. +/// +/// The square is given starting at `(0, 0)` — itself one of the two turns. +/// Boost starts at the *other* one, `(10, 0)`, because `(0, 0)` terminates the +/// square's last segment and so comes last. +#[test] +fn a_union_ring_starts_at_the_first_turn_along_the_first_operand() { + let square: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 0.0) + ]]; + let triangle: Polygon

= polygon![[(0.0, 0.0), (10.0, 0.0), (5.0, -8.0), (0.0, 0.0)]]; + + let expected = vec![ + (10.0, 0.0), + (5.0, -8.0), + (0.0, 0.0), + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + ]; + assert_eq!(vertices(&union_poly(&square, &triangle).unwrap()), expected); + + // Rotating the square so it no longer starts at a turn must not move the + // answer: the same turn is still the first one along its boundary. + let rotated: Polygon

= polygon![[ + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 0.0), + (0.0, 10.0) + ]]; + assert_eq!( + vertices(&union_poly(&rotated, &triangle).unwrap()), + expected + ); +} From 5cbd230ccd92dbd8d2c7b81fb91fe1423d8c650b Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 17:24:55 -0600 Subject: [PATCH 11/20] fix(overlay): clean the traversed ring, and order turns by both operands MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two more pieces of Boost's overlay, both read out of its source rather than inferred from its output. Together they take unions whose operands share a collinear edge from 8 of 112 exact against C++ Boost 1.83 to 89. **`clean_closing_dups_and_spikes`.** `traverse_with_operation` runs it over every ring the traversal produces, and it erases the ring's first point while that point is collinear with its two neighbours, repeating until it is not. A ring starts at a turn, and where two operands share an edge a turn need not be a corner — the outline can run straight through the very point the traversal seeded from. This is why Boost's union of two squares sharing an edge has five distinct vertices where the obvious answer has six. Only the *start* is cleaned, which is what makes the behaviour look arbitrary from outside: the identical straight-through vertex at the far end of the same shared edge survives. And it has to run once the ring is in its final winding, because reversing a ring changes which vertex the erasure leaves at the front — so it moved into `assemble`, next to `orient_ring`, and applies to the rings the traversal assembled rather than to one copied whole from an operand. **Turn order is a pair.** `get_turns` walks the first operand's sections in the outer loop and the second's in the inner, so its turns come out ordered by the pair. Two turns on the same stretch of the first operand are separated by where they sit on the *second* — and rotating the second operand moves the answer, which a single position counter cannot express. `Node::arrival` is now the segment index per operand, with the fraction demoted to a third key so it cannot outrank the second operand. Measured over 312 valid polygon pairs — every rotation of squares and rectangles sharing an edge in whole and in part, meeting at a corner, crossing, contained, disjoint, plus radial polygons: before after collinear edge (112) 8 exact 89 exact crossings only (200) 169 exact 173 exact, 26 of the 27 remaining differ only in the last bits of an irrational crossing The 23 that remain share a collinear edge and need what is not here: `get_turns`' sectionalisation, which groups monotone runs so that turns Boost treats as tied are not, and its collinear turn classification, which decides whether a same-direction overlap contributes a turn point at all. Gate: 1608 tests pass, fmt and clippy clean. --- crates/geometry-overlay/src/assemble.rs | 68 ++++++++++- .../geometry-overlay/src/operation/areal.rs | 99 ++++++++++------ .../geometry-overlay/tests/overlay_parity.rs | 110 ++++++++++++++++++ 3 files changed, 242 insertions(+), 35 deletions(-) diff --git a/crates/geometry-overlay/src/assemble.rs b/crates/geometry-overlay/src/assemble.rs index 50682ac..f857aa8 100644 --- a/crates/geometry-overlay/src/assemble.rs +++ b/crates/geometry-overlay/src/assemble.rs @@ -67,9 +67,25 @@ use crate::surface_point::point_on_surface; pub fn assemble_multipolygon

(rings: Vec>) -> MultiPolygon> where P: PointMut + Default + Copy, - P::Scalar: CoordinateScalar, + P::Scalar: CoordinateScalar + Into, ::Family: SameAs, { + assemble_traced(rings.into_iter().map(|ring| (ring, false)).collect()) +} + +/// As [`assemble_multipolygon`], but told which rings the traversal assembled +/// from turns. +/// +/// Boost cleans exactly those (`traverse_with_operation` calls +/// `clean_closing_dups_and_spikes` on each ring it traverses) and emits a ring +/// copied whole from one operand untouched. +pub(crate) fn assemble_traced

(rings: Vec<(Ring

, bool)>) -> MultiPolygon> +where + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, + ::Family: SameAs, +{ + let (rings, traced): (Vec>, Vec) = rings.into_iter().unzip(); // Classify each ring by containment depth rather than by signed area: // even depths are filled exteriors and odd depths are holes. let n = rings.len(); @@ -97,6 +113,9 @@ where outer_slot[i] = Some(polygons.len()); let mut outer = slots[i].take().unwrap(); orient_ring(&mut outer, true); + if traced[i] { + clean_closing_dups_and_spikes(&mut outer); + } polygons.push(Polygon::new(outer)); } } @@ -110,6 +129,9 @@ where .expect("the immediate parent of an odd-depth ring has even depth"); let mut hole = slots[i].take().unwrap(); orient_ring(&mut hole, false); + if traced[i] { + clean_closing_dups_and_spikes(&mut hole); + } polygons[slot].inners.push(hole); } } @@ -117,6 +139,50 @@ where MultiPolygon(polygons) } +/// Drop the ring's first point while the outline runs straight through it. +/// +/// C++: `clean_closing_dups_and_spikes`, which `traverse_with_operation` +/// applies to every ring the traversal produces. Boost starts a ring at a +/// turn, and a turn is not always a corner — where two operands share a +/// collinear edge, the outline can pass straight through the very point the +/// traversal seeded from. Boost erases it, re-closes, and repeats. +/// +/// This is why Boost's union of two squares sharing an edge has five distinct +/// vertices where the naive answer has six: the sixth is the start, and it sat +/// in the middle of a straight side. Note it is only ever the *start* — the +/// identical straight-through vertex at the far end of the same shared edge +/// stays. +/// +/// The test is Boost's `point_is_collinear`, which is the side test alone: a +/// spike, where the outline doubles back, is collinear too and goes the same +/// way. +fn clean_closing_dups_and_spikes

(ring: &mut Ring

) +where + P: PointTrait + Copy, + P::Scalar: Into, +{ + // C++: `minimum_ring_size`, which is four. + const MINIMUM_CLOSED_RING: usize = 4; + + let points = &mut ring.0; + while points.len() > MINIMUM_CLOSED_RING { + let at = |index: usize| -> (f64, f64) { + let point = &points[index]; + (point.get::<0>().into(), point.get::<1>().into()) + }; + let (fx, fy) = at(0); + let (sx, sy) = at(1); + let (ux, uy) = at(points.len() - 2); + let side = (fx - ux) * (sy - uy) - (fy - uy) * (sx - ux); + if side != 0.0 { + return; + } + points.remove(0); + points.pop(); + points.push(points[0]); + } +} + fn containment_depth(parents: &[Option], mut index: usize) -> usize { let mut depth = 0; while let Some(parent) = parents[index] { diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 54fac89..54a1415 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -24,7 +24,7 @@ use geometry_trait::{ MultiPolygon as MultiPolygonTrait, Point, PointMut, Polygon as PolygonTrait, Ring as RingTrait, }; -use crate::assemble::assemble_multipolygon; +use crate::assemble::assemble_traced; use crate::operation::OverlayError; use crate::predicate::segment_intersection::{SegmentIntersection, segment_intersection}; @@ -153,17 +153,26 @@ struct Node

{ /// starts each output ring at a turn, so the tracer needs to know which /// nodes are turns; see `push_ring`. is_turn: bool, - /// Where this node sits along the operands' boundaries, first operand - /// first — and, where it lands exactly on a vertex, counted as the *end* - /// of the segment arriving there rather than the start of the one leaving. + /// Where this node sits along **each** operand's boundary — and, where it + /// lands exactly on a vertex, counted as the *end* of the segment arriving + /// there rather than the start of the one leaving. /// - /// That normalisation is Boost's, and it is the whole difference between - /// the two: `get_turns` attaches an intersection at a segment endpoint to - /// the segment it terminates, so a turn on the first operand's *first* - /// vertex is the last position on that ring, not the first. Ordering by - /// creation instead put it first, and every union of two rings sharing an - /// edge came out started at the wrong end. See `push_ring`. - arrival: usize, + /// That normalisation is Boost's: `get_turns` attaches an intersection at + /// a segment endpoint to the segment it terminates, so a turn on an + /// operand's *first* vertex is the last position on that ring, not the + /// first. + /// + /// Both entries matter. Boost walks the first operand's sections in the + /// outer loop and the second operand's in the inner, so its turns come out + /// ordered by the pair, and two turns on the same stretch of the first + /// operand are separated by where they sit on the second. `usize::MAX` + /// means the operand's boundary does not pass through this node, which is + /// true of every vertex that is not a turn. + arrival: [usize; 2], + /// How far along that segment. It orders two turns only once both segment + /// indices have tied — the fraction must not outrank the second operand, + /// or two turns sharing one edge come out in the wrong order. + offset: [f64; 2], } #[derive(Debug, Clone, Copy, PartialEq, Eq)] @@ -260,21 +269,18 @@ where let mut nodes = Vec::new(); let mut candidates = Vec::new(); - // One counter across both operands, so every position on the first - // operand's boundary precedes every position on the second's. - let mut arrival = 0; append_atomic_edges( &mut first_segments, &mut nodes, &mut candidates, - &mut arrival, + 0, snap_tolerance, ); append_atomic_edges( &mut second_segments, &mut nodes, &mut candidates, - &mut arrival, + 1, snap_tolerance, ); @@ -320,7 +326,7 @@ where } let rings = trace_rings(&nodes, &boundary, snap_tolerance)?; - Ok(assemble_multipolygon(rings)) + Ok(assemble_traced(rings)) } fn ring_coordinates(ring: &R) -> Vec @@ -396,13 +402,13 @@ fn append_atomic_edges

( segments: &mut [SourceSegment

], nodes: &mut Vec>, output: &mut Vec, - arrival: &mut usize, + operand: usize, tolerance: f64, ) where P: Point + Copy, P::Scalar: Into, { - for segment in segments { + for (index, segment) in segments.iter_mut().enumerate() { segment .splits .sort_by(|left, right| left.0.total_cmp(&right.0)); @@ -413,9 +419,9 @@ fn append_atomic_edges

( // Only the far end of a split counts as an arrival, which is what // pushes a ring's first vertex to the end of its own ring: it is // reached as the last segment's endpoint, not the first's start. - if nodes[end].arrival == usize::MAX { - nodes[end].arrival = *arrival; - *arrival += 1; + if nodes[end].arrival[operand] == usize::MAX { + nodes[end].arrival[operand] = index; + nodes[end].offset[operand] = pair[1].0; } if start != end { output.push(Edge { start, end }); @@ -443,24 +449,28 @@ where point, coordinate, is_turn, - arrival: usize::MAX, + arrival: [usize::MAX; 2], + offset: [0.0; 2], }); nodes.len() - 1 } +type TracedRing

= (Ring

, bool); + fn trace_rings

( nodes: &[Node

], edges: &[Edge], tolerance: f64, -) -> Result>, OverlayError> +) -> Result>, OverlayError> where P: Point + Copy, P::Scalar: Into, { let mut used = alloc::vec![false; edges.len()]; // Each ring is kept with the node it starts at, so the whole set can be - // put back into source order below. - let mut rings: Vec<(usize, Ring

)> = Vec::new(); + // put back into source order below, and with whether the traversal + // assembled it from turns rather than copying it whole from one operand. + let mut rings: Vec<(usize, Ring

, bool)> = Vec::new(); for seed in 0..edges.len() { if used[seed] { continue; @@ -502,8 +512,11 @@ where // polygons in the result — and Boost's is not the order the seeds happened // to fall in. It traces from turns in source order, so ordering the rings // by the node each one starts at reproduces it. - rings.sort_by_key(|&(start, _)| start); - Ok(rings.into_iter().map(|(_, ring)| ring).collect()) + rings.sort_by_key(|&(start, _, _)| start); + Ok(rings + .into_iter() + .map(|(_, ring, traced)| (ring, traced)) + .collect()) } /// Append one traced cycle, dropping it when it encloses no area. @@ -516,7 +529,7 @@ where /// because a consumer that simplifies the ring afterwards will pin its first /// vertex and the choice reaches the output. fn push_ring

( - rings: &mut Vec<(usize, Ring

)>, + rings: &mut Vec<(usize, Ring

, bool)>, nodes: &[Node

], node_indices: &[usize], tolerance: f64, @@ -538,12 +551,20 @@ fn push_ring

( // Boost begins each output ring at the first *turn* along the first // operand's boundary — `Node::arrival`, which is that order with Boost's // endpoint normalisation applied. + // Boost walks the first operand's sections in the outer loop and the + // second's in the inner, so its turns are ordered by the pair. let first_turn_by_arrival = cycle .iter() .copied() .enumerate() .filter(|&(_, index)| nodes[index].is_turn) - .min_by_key(|&(_, index)| nodes[index].arrival) + .min_by(|&(_, left), &(_, right)| { + let (a, b) = (&nodes[left], &nodes[right]); + a.arrival[0] + .cmp(&b.arrival[0]) + .then(a.arrival[1].cmp(&b.arrival[1])) + .then(a.offset[0].total_cmp(&b.offset[0])) + }) .map(|(position, _)| position); // A ring with no turn was copied whole from one operand, and keeps that // operand's own starting vertex rather than whichever end of it the @@ -568,7 +589,14 @@ fn push_ring

( if let Some(&first) = points.first() { points.push(first); } - rings.push((cycle[first_turn], Ring::from_vec(points))); + // The ring is cleaned once it is in its final winding, not here: which + // vertex `clean_closing_dups_and_spikes` leaves at the front depends on + // the direction the ring runs in, and that is decided in `assemble`. + rings.push(( + cycle[first_turn], + Ring::from_vec(points), + first_turn_by_arrival.is_some(), + )); } fn next_edge

(nodes: &[Node

], edges: &[Edge], used: &[bool], incoming: Edge) -> Option @@ -671,19 +699,22 @@ mod tests { point: P::new(0.0, 0.0), coordinate: Coordinate { x: 0.0, y: 0.0 }, is_turn: false, - arrival: 0, + arrival: [0, 0], + offset: [0.0; 2], }, Node { point: P::new(1.0, 0.0), coordinate: Coordinate { x: 1.0, y: 0.0 }, is_turn: false, - arrival: 1, + arrival: [1, 1], + offset: [0.0; 2], }, Node { point: P::new(2.0, 0.0), coordinate: Coordinate { x: 2.0, y: 0.0 }, is_turn: false, - arrival: 2, + arrival: [2, 2], + offset: [0.0; 2], }, ]; let edges = [ diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 4540f04..d2df106 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -309,3 +309,113 @@ fn a_union_ring_starts_at_the_first_turn_along_the_first_operand() { expected ); } + +// ---- Unions whose operands share a collinear edge ------------------------ +// +// Two more pieces of Boost, both taken from its source rather than guessed at: +// +// * `traverse_with_operation` runs `clean_closing_dups_and_spikes` over every +// ring it traverses, which erases the ring's first point while the outline +// runs straight through it. A ring starts at a turn, and where two operands +// share an edge a turn need not be a corner. +// * `get_turns` walks the first operand's sections in the outer loop and the +// second's in the inner, so two turns on the same stretch of the first +// operand are ordered by where they sit on the *second*. +// +// Reference values from C++ Boost 1.83 on the same input. + +/// Two squares sharing a whole edge. The traversal starts at `(10, 10)` — the +/// first turn — and that point sits in the middle of the union's straight top +/// side, so Boost erases it and the ring begins at `(20, 10)`. Note the +/// identical straight-through point at the *other* end of the shared edge, +/// `(10, 0)`, survives: only the start is cleaned. +#[test] +fn a_shared_edge_loses_the_ring_start_it_ran_straight_through() { + let left: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 10.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 0.0) + ]]; + let right: Polygon

= polygon![[ + (10.0, 0.0), + (10.0, 10.0), + (20.0, 10.0), + (20.0, 0.0), + (10.0, 0.0) + ]]; + assert_eq!( + vertices(&union_poly(&left, &right).unwrap()), + vec![ + (20.0, 10.0), + (20.0, 0.0), + (10.0, 0.0), + (0.0, 0.0), + (0.0, 10.0), + (20.0, 10.0), + ] + ); +} + +/// A square and a rectangle overlapping along part of one side, so both turns +/// lie on the *same* segment of the first operand. Which of them starts the +/// ring is then decided by the second operand — and rotating it moves the +/// answer, which is why the second operand's position has to be part of the +/// ordering and the fraction along the first must not outrank it. +#[test] +fn two_turns_on_one_segment_are_ordered_by_the_second_operand() { + let square: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 0.0), + (0.0, 0.0) + ]]; + // Starting at (100, 30): the second operand's last segment ends there, + // which puts (100, 70) ahead of it. + let from_bottom: Polygon

= polygon![[ + (100.0, 30.0), + (100.0, 70.0), + (200.0, 70.0), + (200.0, 30.0), + (100.0, 30.0) + ]]; + assert_eq!( + vertices(&union_poly(&square, &from_bottom).unwrap()), + vec![ + (100.0, 70.0), + (200.0, 70.0), + (200.0, 30.0), + (100.0, 30.0), + (100.0, 0.0), + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 70.0), + ] + ); + + // Rotated, (100, 30) now ends an earlier segment and takes the start. + let from_top: Polygon

= polygon![[ + (100.0, 70.0), + (200.0, 70.0), + (200.0, 30.0), + (100.0, 30.0), + (100.0, 70.0) + ]]; + assert_eq!( + vertices(&union_poly(&square, &from_top).unwrap()), + vec![ + (100.0, 30.0), + (100.0, 0.0), + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 70.0), + (200.0, 70.0), + (200.0, 30.0), + (100.0, 30.0), + ] + ); +} From 5cd2c06c5d660ceefc42361638bb8769dac95ab7 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 17:47:32 -0600 Subject: [PATCH 12/20] fix(overlay): do not emit a point the walked operand runs straight past MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Between two turns Boost copies the vertices of the single operand it is walking, so a point that operand has no vertex at is passed over. This arrangement splits a segment wherever anything touches it, so the *other* operand's corner landing part-way along became a vertex of the traced ring that Boost never emits. Edges now carry who runs along them and nodes know whose vertex they are, which is what the rule needs: along a stretch both operands carry, the one walked is the first, so a point the first operand has no vertex at goes. Two guards keep it honest — the point must be collinear with its neighbours, because anywhere the outline actually turns the point is part of the shape whoever owns it, and a point the first operand does have a vertex at stays. Measured over the same 312 valid polygon pairs, against C++ Boost 1.83 and compared exactly, ring start included: before after shared edge, same way (32) 9 exact 20 exact shared edge, opposite (80) 80 exact 80 exact no shared edge (200) 173 exact 173 exact, 26 of the 27 remaining differing only in an irrational crossing's last bits Two other formulations were tried and dropped, both of which looked right on the cases they were reasoned from and measured worse: dropping any point that is a vertex of exactly one operand (23 differ becomes 32), and dropping the *entry* of a shared stretch (23 becomes 20, and it loses cases this rule keeps). That second one is the signpost for the 12 that remain. They are all same-direction shared edges, where `get_turn_info`'s collinear handler emits one turn for the whole overlap and puts it at the far end (`non_opposite_to_index`), so the near end is not a turn and never reaches the output. Reproducing that needs `get_turns`' own collinear turn handling, not a further rule layered over the arrangement; the limit and its current number are recorded on the function. Gate: 1609 tests pass, fmt and clippy clean. tilemaker's `scripts/verify.sh`: 16 passed, 0 failed, unchanged. --- .../geometry-overlay/src/operation/areal.rs | 146 ++++++++++++++++-- .../geometry-overlay/tests/overlay_parity.rs | 39 +++++ 2 files changed, 175 insertions(+), 10 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 54a1415..2736578 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -173,12 +173,30 @@ struct Node

{ /// indices have tied — the fraction must not outrank the second operand, /// or two turns sharing one edge come out in the wrong order. offset: [f64; 2], + /// Whether this node is a *vertex* of each operand, as against a point + /// that only lies on its boundary because the other operand touches there. + is_vertex: [bool; 2], } -#[derive(Debug, Clone, Copy, PartialEq, Eq)] +#[derive(Debug, Clone, Copy)] struct Edge { start: usize, end: usize, + /// Which operands' boundaries run along this edge. A stretch the two share + /// is carried by both. + /// + /// Boost walks one operand at a time between turns and copies *that* + /// operand's vertices, so a point sitting inside a segment of the operand + /// being walked never reaches the output, whoever else has a vertex there. + /// Reproducing that needs to know who carries each edge — see + /// `drop_points_interior_to_a_walked_segment`. + carried_by: [bool; 2], +} + +impl Edge { + fn joins(&self, other: &Self) -> bool { + self.start == other.start && self.end == other.end + } } /// Execute a polygon Boolean operation through a split-edge arrangement. @@ -285,7 +303,7 @@ where ); let sample_distance = (scale * 1e-8).max(snap_tolerance * 32.0); - let mut boundary = Vec::new(); + let mut boundary: Vec = Vec::new(); for candidate in candidates { let start = nodes[candidate.start].coordinate; let end = nodes[candidate.end].coordinate; @@ -318,10 +336,18 @@ where Edge { start: candidate.end, end: candidate.start, + carried_by: candidate.carried_by, } }; - if !boundary.contains(&edge) { - boundary.push(edge); + // The same stretch reaches here once per operand that carries it, so + // merge rather than drop the second: who carries an edge is what says + // whether a point on it is interior to a walked segment. + match boundary.iter_mut().find(|held| held.joins(&edge)) { + Some(held) => { + held.carried_by[0] |= edge.carried_by[0]; + held.carried_by[1] |= edge.carried_by[1]; + } + None => boundary.push(edge), } } @@ -423,8 +449,20 @@ fn append_atomic_edges

( nodes[end].arrival[operand] = index; nodes[end].offset[operand] = pair[1].0; } + if pair[0].0 == 0.0 || pair[0].0 == 1.0 { + nodes[start].is_vertex[operand] = true; + } + if pair[1].0 == 0.0 || pair[1].0 == 1.0 { + nodes[end].is_vertex[operand] = true; + } if start != end { - output.push(Edge { start, end }); + let mut carried_by = [false; 2]; + carried_by[operand] = true; + output.push(Edge { + start, + end, + carried_by, + }); } } } @@ -451,6 +489,7 @@ where is_turn, arrival: [usize::MAX; 2], offset: [0.0; 2], + is_vertex: [false; 2], }); nodes.len() - 1 } @@ -478,11 +517,15 @@ where let first = edges[seed].start; let mut edge_index = seed; let mut node_indices = alloc::vec![first]; + // `carried[i]` is who carries the edge from `node_indices[i]` to the + // node after it, so it always has one entry fewer. + let mut carried: alloc::vec::Vec<[bool; 2]> = alloc::vec::Vec::new(); for _ in 0..=edges.len() { debug_assert!(!used[edge_index]); used[edge_index] = true; let edge = edges[edge_index]; node_indices.push(edge.end); + carried.push(edge.carried_by); // A node the walk has already stood on closes a ring right here, // not only when the walk returns to the seed. Where two lobes of @@ -496,8 +539,9 @@ where .position(|&index| index == edge.end) { let loop_nodes = node_indices.split_off(start); + let loop_carried = carried.split_off(start); node_indices.push(edge.end); - push_ring(&mut rings, nodes, &loop_nodes, tolerance); + push_ring(&mut rings, nodes, &loop_nodes, &loop_carried, tolerance); } if edge.end == first { @@ -519,6 +563,69 @@ where .collect()) } +/// Remove a point the first operand runs straight past. +/// +/// This arrangement splits a segment wherever anything touches it, so a vertex +/// of one operand landing part-way along a segment of the other becomes a +/// point on the traced ring. Boost's traversal does not work that way: between +/// two turns it copies the vertices of the single operand it is walking, and a +/// point that operand has no vertex at is simply passed over. +/// +/// Along a stretch both operands carry, the operand walked is the first, so +/// that is the one whose vertices survive. Two guards keep the rule honest: +/// the point must be collinear with its neighbours, because anywhere the +/// outline actually turns the point is part of the shape whoever owns it; and +/// a point the first operand does have a vertex at stays, because Boost would +/// copy it. +/// +/// Not every such point is caught. Where the two operands run the same way +/// along a shared edge, Boost emits one turn for the overlap and puts it at +/// the far end, so the near end goes too — 12 of 112 shared-edge pairs still +/// differ on that, and closing them needs `get_turns`' own collinear turn +/// handling rather than a rule over the arrangement. +fn drop_points_interior_to_a_walked_segment

( + nodes: &[Node

], + cycle: &mut Vec, + carried: &mut Vec<[bool; 2]>, +) where + P: Point + Copy, + P::Scalar: Into, +{ + let mut index = 0; + while index < cycle.len() && cycle.len() > 3 { + let length = cycle.len(); + let incoming = carried[(index + length - 1) % length]; + let outgoing = carried[index]; + let node = &nodes[cycle[index]]; + // Boost walks one operand at a time between turns and copies *that* + // operand's vertices. Along a stretch both operands carry, it walks + // the first — so a point the first operand runs straight past, having + // no vertex there, never reaches the output however many vertices the + // second has at it. + let walked_over = incoming[0] && outgoing[0] && !node.is_vertex[0]; + if !walked_over { + index += 1; + continue; + } + let previous = nodes[cycle[(index + length - 1) % length]].coordinate; + let here = node.coordinate; + let next = nodes[cycle[(index + 1) % length]].coordinate; + let side = (here.x - previous.x) * (next.y - previous.y) + - (here.y - previous.y) * (next.x - previous.x); + if side != 0.0 { + index += 1; + continue; + } + // The two edges become one, carried by whoever carried both halves. + let merged = [incoming[0] && outgoing[0], incoming[1] && outgoing[1]]; + cycle.remove(index); + carried.remove(index); + let previous_edge = (index + cycle.len() - 1) % cycle.len(); + carried[previous_edge] = merged; + index = 0; + } +} + /// Append one traced cycle, dropping it when it encloses no area. /// /// The cycle starts wherever the traversal happened to seed, which carries no @@ -532,6 +639,7 @@ fn push_ring

( rings: &mut Vec<(usize, Ring

, bool)>, nodes: &[Node

], node_indices: &[usize], + carried: &[[bool; 2]], tolerance: f64, ) where P: Point + Copy, @@ -547,7 +655,10 @@ fn push_ring

( } // The cycle is closed, so its last index repeats its first. - let cycle = &node_indices[..node_indices.len() - 1]; + let mut cycle = node_indices[..node_indices.len() - 1].to_vec(); + let mut carried = carried.to_vec(); + drop_points_interior_to_a_walked_segment(nodes, &mut cycle, &mut carried); + let cycle = &cycle[..]; // Boost begins each output ring at the first *turn* along the first // operand's boundary — `Node::arrival`, which is that order with Boost's // endpoint normalisation applied. @@ -701,6 +812,7 @@ mod tests { is_turn: false, arrival: [0, 0], offset: [0.0; 2], + is_vertex: [true; 2], }, Node { point: P::new(1.0, 0.0), @@ -708,6 +820,7 @@ mod tests { is_turn: false, arrival: [1, 1], offset: [0.0; 2], + is_vertex: [true; 2], }, Node { point: P::new(2.0, 0.0), @@ -715,12 +828,25 @@ mod tests { is_turn: false, arrival: [2, 2], offset: [0.0; 2], + is_vertex: [true; 2], }, ]; let edges = [ - Edge { start: 0, end: 1 }, - Edge { start: 1, end: 2 }, - Edge { start: 2, end: 0 }, + Edge { + start: 0, + end: 1, + carried_by: [true; 2], + }, + Edge { + start: 1, + end: 2, + carried_by: [true; 2], + }, + Edge { + start: 2, + end: 0, + carried_by: [true; 2], + }, ]; assert!(trace_rings(&nodes, &edges, 1e-10).unwrap().is_empty()); diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index d2df106..2d5931f 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -419,3 +419,42 @@ fn two_turns_on_one_segment_are_ordered_by_the_second_operand() { ] ); } + +/// A square and a rectangle overlapping along part of one side, running the +/// same way round. The rectangle's corner at `(50, 100)` lands part-way along +/// the square's top edge, and this arrangement splits there — but Boost walks +/// the square across it without stopping, because the square has no vertex +/// there, so the point never reaches the output. +/// +/// The square's own corner at `(100, 100)` is a different matter: it is a +/// vertex of the operand being walked, and it survives — until the ring-start +/// cleaning takes it, which is why the answer begins at `(150, 100)`. +#[test] +fn a_point_the_walked_operand_runs_straight_past_is_not_emitted() { + let square: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 0.0), + (0.0, 0.0) + ]]; + let overlapping: Polygon

= polygon![[ + (50.0, 100.0), + (150.0, 100.0), + (150.0, 50.0), + (50.0, 50.0), + (50.0, 100.0) + ]]; + assert_eq!( + vertices(&union_poly(&square, &overlapping).unwrap()), + vec![ + (150.0, 100.0), + (150.0, 50.0), + (100.0, 50.0), + (100.0, 0.0), + (0.0, 0.0), + (0.0, 100.0), + (150.0, 100.0), + ] + ); +} From e70c42395c82613d46f12080ffbc7a1786b21fe5 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 18:24:09 -0600 Subject: [PATCH 13/20] fix(overlay): append a turn the way the traversal does, and order turns by section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous commit approximated the traversal's vertex handling with a rule over the finished cycle, and left twelve shared-edge pairs it could not reach. This replaces the rule with the mechanism, read out of Boost's source: the traversal appends a *turn* through `append_no_collinear`, which drops any point the new one now runs straight through, and copies the vertices between two turns through `copy_segments`, which does no such check. So a shared edge's near end goes not because anything decided it was interior to a walked segment, but because the far end of the overlap is a turn collinear with it and takes its place — while an operand's own straight-through vertex, merely copied, stays. Two things had to come with it. `append_no_collinear` and `clean_closing_dups_and_spikes` both judge the point *before* the one they are given, so which vertex of a straight run survives depends on the direction the ring runs in. Edges are now oriented with the filled side on the right, the winding Boost's traversal produces. And the turn order needed `sectionalize`. `get_turns` partitions each operand into sections — runs of consecutive segments heading the same way in both dimensions, capped at ten — and walks the section pairs, so the section pair outranks the segment index: two turns on one segment of the first operand can come out in either order depending on which section of the second operand each one meets. Ring starts are ordered by that pair first, and two rings starting at the same node by which operand carries the edge leaving it, since `iterate` tries operation 0 before operation 1. `Node::is_vertex` was only there for the rule this removes, and goes with it. Measured over 312 valid polygon pairs against C++ Boost 1.83, compared exactly, ring start and ring order included: before after shares a collinear edge (112) 100 exact 112 exact meets only at crossings (200) 173 exact 173 exact differing structurally 13 0 The 27 that still differ all sit on an irrational crossing and differ in its last bits alone, worst 2.84e-14 across the set. Gate: 1612 tests pass, fmt and clippy clean. tilemaker's `scripts/verify.sh`: 16 passed, 0 failed. --- .../geometry-overlay/src/operation/areal.rs | 266 +++++++++++------- .../geometry-overlay/tests/overlay_parity.rs | 144 +++++++++- 2 files changed, 305 insertions(+), 105 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 2736578..38367e4 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -112,6 +112,9 @@ struct SourceSegment

{ /// pushed by the intersection sweep is, including one that lands on an /// endpoint. splits: Vec<(f64, P, bool)>, + /// Which monotone run of its ring this segment belongs to. C++: + /// `sectionalize`, whose sections are what `get_turns` iterates over. + section: usize, } impl

SourceSegment

@@ -119,11 +122,12 @@ where P: Point + Copy, P::Scalar: Into, { - fn new(start: P, end: P) -> Self { + fn new(start: P, end: P, section: usize) -> Self { Self { start, end, splits: alloc::vec![(0.0, start, false), (1.0, end, false)], + section, } } @@ -169,13 +173,15 @@ struct Node

{ /// means the operand's boundary does not pass through this node, which is /// true of every vertex that is not a turn. arrival: [usize; 2], + /// The section of each operand that reaches this node, which outranks the + /// segment: `get_turns` partitions both operands into sections first and + /// walks the pairs, so two turns in the same pair of sections keep their + /// segment order while turns in different pairs do not. + section: [usize; 2], /// How far along that segment. It orders two turns only once both segment /// indices have tied — the fraction must not outrank the second operand, /// or two turns sharing one edge come out in the wrong order. offset: [f64; 2], - /// Whether this node is a *vertex* of each operand, as against a point - /// that only lies on its boundary because the other operand touches there. - is_vertex: [bool; 2], } #[derive(Debug, Clone, Copy)] @@ -330,7 +336,14 @@ where if left_result == right_result { continue; } - let edge = if left_result { + // Oriented so the filled side is on the right, which walks an outer + // ring clockwise and a hole counter-clockwise — the directions Boost's + // traversal produces, and the ones `append_no_collinear` and + // `clean_closing_dups_and_spikes` are written against. Both look at + // the point *before* the one they judge, so a ring traced the other + // way round drops the vertex at the far end of a straight run instead + // of the near one. + let edge = if right_result { candidate } else { Edge { @@ -378,10 +391,11 @@ where P::Scalar: Into, { let mut segments = Vec::new(); + let mut sections = Sectionizer::new(0); for polygon in multi_polygon.polygons() { - append_ring_segments(polygon.exterior(), &mut segments); + append_ring_segments(polygon.exterior(), &mut segments, &mut sections); for ring in polygon.interiors() { - append_ring_segments(ring, &mut segments); + append_ring_segments(ring, &mut segments, &mut sections); } } segments @@ -394,15 +408,82 @@ where P::Scalar: Into, { let mut segments = Vec::new(); - append_ring_segments(polygon.exterior(), &mut segments); + let mut sections = Sectionizer::new(0); + append_ring_segments(polygon.exterior(), &mut segments, &mut sections); for ring in polygon.interiors() { - append_ring_segments(ring, &mut segments); + append_ring_segments(ring, &mut segments, &mut sections); } segments } -fn append_ring_segments(ring: &R, output: &mut Vec>) -where +/// C++: `sectionalize`'s cap, "defaults to 10, this seems to give the fastest +/// results". +const MAX_SEGMENTS_PER_SECTION: usize = 10; + +/// A section is a run of consecutive segments heading the same way in both +/// dimensions. C++: `sectionalize`, which starts a new one whenever the pair +/// of signs changes, or the run grows past `max_count`. Sections do not span +/// rings. +struct Sectionizer { + next: usize, + directions: Option<(i8, i8)>, + count: usize, +} + +impl Sectionizer { + fn new(next: usize) -> Self { + Self { + next, + directions: None, + count: 0, + } + } + + fn start_ring(&mut self) { + if self.count > 0 { + self.next += 1; + } + self.directions = None; + self.count = 0; + } + + fn section_for

(&mut self, start: &P, end: &P) -> usize + where + P: Point, + P::Scalar: Into, + { + let sign = |a: f64, b: f64| -> i8 { + if b > a { + 1 + } else if b < a { + -1 + } else { + 0 + } + }; + let directions = ( + sign(start.get::<0>().into(), end.get::<0>().into()), + sign(start.get::<1>().into(), end.get::<1>().into()), + ); + if self.count > 0 + && (Some(directions) != self.directions || self.count > MAX_SEGMENTS_PER_SECTION) + { + self.next += 1; + self.count = 0; + } + if self.count == 0 { + self.directions = Some(directions); + } + self.count += 1; + self.next + } +} + +fn append_ring_segments( + ring: &R, + output: &mut Vec>, + sections: &mut Sectionizer, +) where R: RingTrait, P: Point + Copy, P::Scalar: Into, @@ -411,16 +492,17 @@ where if points.len() < 2 { return; } + sections.start_ring(); for pair in points.windows(2) { if points_differ(&pair[0], &pair[1]) { - output.push(SourceSegment::new(pair[0], pair[1])); + let section = sections.section_for(&pair[0], &pair[1]); + output.push(SourceSegment::new(pair[0], pair[1], section)); } } - if points_differ(points.last().expect("nonempty"), &points[0]) { - output.push(SourceSegment::new( - *points.last().expect("nonempty"), - points[0], - )); + let last = *points.last().expect("nonempty"); + if points_differ(&last, &points[0]) { + let section = sections.section_for(&last, &points[0]); + output.push(SourceSegment::new(last, points[0], section)); } } @@ -435,6 +517,7 @@ fn append_atomic_edges

( P::Scalar: Into, { for (index, segment) in segments.iter_mut().enumerate() { + let section = segment.section; segment .splits .sort_by(|left, right| left.0.total_cmp(&right.0)); @@ -448,12 +531,7 @@ fn append_atomic_edges

( if nodes[end].arrival[operand] == usize::MAX { nodes[end].arrival[operand] = index; nodes[end].offset[operand] = pair[1].0; - } - if pair[0].0 == 0.0 || pair[0].0 == 1.0 { - nodes[start].is_vertex[operand] = true; - } - if pair[1].0 == 0.0 || pair[1].0 == 1.0 { - nodes[end].is_vertex[operand] = true; + nodes[end].section[operand] = section; } if start != end { let mut carried_by = [false; 2]; @@ -488,8 +566,8 @@ where coordinate, is_turn, arrival: [usize::MAX; 2], + section: [usize::MAX; 2], offset: [0.0; 2], - is_vertex: [false; 2], }); nodes.len() - 1 } @@ -509,7 +587,7 @@ where // Each ring is kept with the node it starts at, so the whole set can be // put back into source order below, and with whether the traversal // assembled it from turns rather than copying it whole from one operand. - let mut rings: Vec<(usize, Ring

, bool)> = Vec::new(); + let mut rings: Vec<(usize, bool, Ring

, bool)> = Vec::new(); for seed in 0..edges.len() { if used[seed] { continue; @@ -556,73 +634,46 @@ where // polygons in the result — and Boost's is not the order the seeds happened // to fall in. It traces from turns in source order, so ordering the rings // by the node each one starts at reproduces it. - rings.sort_by_key(|&(start, _, _)| start); + rings.sort_by_key(|&(start, second_operand, _, _)| (start, second_operand)); Ok(rings .into_iter() - .map(|(_, ring, traced)| (ring, traced)) + .map(|(_, _, ring, traced)| (ring, traced)) .collect()) } -/// Remove a point the first operand runs straight past. -/// -/// This arrangement splits a segment wherever anything touches it, so a vertex -/// of one operand landing part-way along a segment of the other becomes a -/// point on the traced ring. Boost's traversal does not work that way: between -/// two turns it copies the vertices of the single operand it is walking, and a -/// point that operand has no vertex at is simply passed over. +/// Append a turn point, dropping whatever it now runs straight through. /// -/// Along a stretch both operands carry, the operand walked is the first, so -/// that is the one whose vertices survive. Two guards keep the rule honest: -/// the point must be collinear with its neighbours, because anywhere the -/// outline actually turns the point is part of the shape whoever owns it; and -/// a point the first operand does have a vertex at stays, because Boost would -/// copy it. +/// C++: `append_no_collinear`. Once the point is on, any point before it that +/// the new one continues the line of is redundant and comes off — repeatedly, +/// because removing one can leave the next in the same position. /// -/// Not every such point is caught. Where the two operands run the same way -/// along a shared edge, Boost emits one turn for the overlap and puts it at -/// the far end, so the near end goes too — 12 of 112 shared-edge pairs still -/// differ on that, and closing them needs `get_turns`' own collinear turn -/// handling rather than a rule over the arrangement. -fn drop_points_interior_to_a_walked_segment

( - nodes: &[Node

], - cycle: &mut Vec, - carried: &mut Vec<[bool; 2]>, -) where +/// Boost applies this to turn points only. The ring vertices copied between +/// two turns go on through `append_no_dups_or_spikes`, which takes out +/// duplicates and spikes but leaves a vertex that merely continues straight, +/// so an operand's own collinear vertex survives while a turn's does not. +fn append_no_collinear

(points: &mut Vec

, point: P) +where P: Point + Copy, P::Scalar: Into, { - let mut index = 0; - while index < cycle.len() && cycle.len() > 3 { - let length = cycle.len(); - let incoming = carried[(index + length - 1) % length]; - let outgoing = carried[index]; - let node = &nodes[cycle[index]]; - // Boost walks one operand at a time between turns and copies *that* - // operand's vertices. Along a stretch both operands carry, it walks - // the first — so a point the first operand runs straight past, having - // no vertex there, never reaches the output however many vertices the - // second has at it. - let walked_over = incoming[0] && outgoing[0] && !node.is_vertex[0]; - if !walked_over { - index += 1; - continue; + let at = |p: &P| (p.get::<0>().into(), p.get::<1>().into()); + let (x, y) = at(&point); + if points.len() == 1 { + let (fx, fy) = at(&points[0]); + if fx == x && fy == y { + return; } - let previous = nodes[cycle[(index + length - 1) % length]].coordinate; - let here = node.coordinate; - let next = nodes[cycle[(index + 1) % length]].coordinate; - let side = (here.x - previous.x) * (next.y - previous.y) - - (here.y - previous.y) * (next.x - previous.x); - if side != 0.0 { - index += 1; - continue; + } + points.push(point); + while points.len() >= 3 { + let (ax, ay) = at(&points[points.len() - 3]); + let (bx, by) = at(&points[points.len() - 2]); + if (bx - ax) * (y - ay) - (by - ay) * (x - ax) != 0.0 { + return; } - // The two edges become one, carried by whoever carried both halves. - let merged = [incoming[0] && outgoing[0], incoming[1] && outgoing[1]]; - cycle.remove(index); - carried.remove(index); - let previous_edge = (index + cycle.len() - 1) % cycle.len(); - carried[previous_edge] = merged; - index = 0; + let last = points.pop().expect("just pushed"); + points.pop(); + points.push(last); } } @@ -636,7 +687,7 @@ fn drop_points_interior_to_a_walked_segment

( /// because a consumer that simplifies the ring afterwards will pin its first /// vertex and the choice reaches the output. fn push_ring

( - rings: &mut Vec<(usize, Ring

, bool)>, + rings: &mut Vec<(usize, bool, Ring

, bool)>, nodes: &[Node

], node_indices: &[usize], carried: &[[bool; 2]], @@ -655,10 +706,7 @@ fn push_ring

( } // The cycle is closed, so its last index repeats its first. - let mut cycle = node_indices[..node_indices.len() - 1].to_vec(); - let mut carried = carried.to_vec(); - drop_points_interior_to_a_walked_segment(nodes, &mut cycle, &mut carried); - let cycle = &cycle[..]; + let cycle = &node_indices[..node_indices.len() - 1]; // Boost begins each output ring at the first *turn* along the first // operand's boundary — `Node::arrival`, which is that order with Boost's // endpoint normalisation applied. @@ -671,8 +719,10 @@ fn push_ring

( .filter(|&(_, index)| nodes[index].is_turn) .min_by(|&(_, left), &(_, right)| { let (a, b) = (&nodes[left], &nodes[right]); - a.arrival[0] - .cmp(&b.arrival[0]) + a.section[0] + .cmp(&b.section[0]) + .then(a.section[1].cmp(&b.section[1])) + .then(a.arrival[0].cmp(&b.arrival[0])) .then(a.arrival[1].cmp(&b.arrival[1])) .then(a.offset[0].total_cmp(&b.offset[0])) }) @@ -692,19 +742,43 @@ fn push_ring

( }; let first_turn = first_turn_by_arrival.or_else(first_node).unwrap_or(0); - let mut points: Vec

= cycle[first_turn..] - .iter() - .chain(&cycle[..first_turn]) - .map(|&index| nodes[index].point) - .collect(); + // C++: the traversal appends a *turn* point with `append_no_collinear` and + // the ring vertices between turns with `copy_segments`, which does not + // check for collinearity. So a turn that carries the outline straight on + // replaces the point before it, and a vertex of the operand being walked + // never does. + // + // This is what keeps the other operand's corner out of the result where + // the two run along the same edge: the traversal reaches that corner as a + // turn, appends it, and then the next turn — the far end of the shared + // stretch — is collinear with it and takes its place. + let mut points: Vec

= Vec::with_capacity(cycle.len() + 1); + for &index in cycle[first_turn..].iter().chain(&cycle[..first_turn]) { + let node = &nodes[index]; + if !node.is_turn { + points.push(node.point); + continue; + } + append_no_collinear(&mut points, node.point); + } + // The traversal closes a ring by arriving back at the turn it started + // from, and that arrival is an append like any other — which is exactly + // where the point before it goes, when the start carries the outline + // straight on through it. if let Some(&first) = points.first() { - points.push(first); + append_no_collinear(&mut points, first); } // The ring is cleaned once it is in its final winding, not here: which // vertex `clean_closing_dups_and_spikes` leaves at the front depends on // the direction the ring runs in, and that is decided in `assemble`. + // Two rings can begin at the same node — where the result touches itself + // at a point, both lobes start there. Boost separates them by operand: + // `iterate` tries operation 0 before operation 1 at a turn, so the lobe + // traced along the first operand is emitted first. + let leaves_along_first_operand = carried.get(first_turn).is_some_and(|edge| edge[0]); rings.push(( cycle[first_turn], + !leaves_along_first_operand, Ring::from_vec(points), first_turn_by_arrival.is_some(), )); @@ -811,24 +885,24 @@ mod tests { coordinate: Coordinate { x: 0.0, y: 0.0 }, is_turn: false, arrival: [0, 0], + section: [0, 0], offset: [0.0; 2], - is_vertex: [true; 2], }, Node { point: P::new(1.0, 0.0), coordinate: Coordinate { x: 1.0, y: 0.0 }, is_turn: false, arrival: [1, 1], + section: [1, 1], offset: [0.0; 2], - is_vertex: [true; 2], }, Node { point: P::new(2.0, 0.0), coordinate: Coordinate { x: 2.0, y: 0.0 }, is_turn: false, arrival: [2, 2], + section: [2, 2], offset: [0.0; 2], - is_vertex: [true; 2], }, ]; let edges = [ diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 2d5931f..76c5de7 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -421,16 +421,11 @@ fn two_turns_on_one_segment_are_ordered_by_the_second_operand() { } /// A square and a rectangle overlapping along part of one side, running the -/// same way round. The rectangle's corner at `(50, 100)` lands part-way along -/// the square's top edge, and this arrangement splits there — but Boost walks -/// the square across it without stopping, because the square has no vertex -/// there, so the point never reaches the output. -/// -/// The square's own corner at `(100, 100)` is a different matter: it is a -/// vertex of the operand being walked, and it survives — until the ring-start -/// cleaning takes it, which is why the answer begins at `(150, 100)`. +/// same way round. Both ends of the overlap are turns, and both carry the +/// outline straight on — so each is appended and then replaced by the next +/// turn along, which is `append_no_collinear` doing what Boost does. #[test] -fn a_point_the_walked_operand_runs_straight_past_is_not_emitted() { +fn a_turn_that_carries_the_outline_straight_on_is_replaced() { let square: Polygon

= polygon![[ (0.0, 0.0), (0.0, 100.0), @@ -458,3 +453,134 @@ fn a_point_the_walked_operand_runs_straight_past_is_not_emitted() { ] ); } + +/// The same shape the other way up: the rectangle straddles the square, and +/// the overlap runs down one side. `(100, 100)` is the square's own corner and +/// still goes, because the turn after it — the far end of the overlap — is +/// collinear with it. +#[test] +fn the_walked_operands_own_corner_goes_too_when_a_turn_follows_it_straight() { + let square: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 0.0), + (0.0, 0.0) + ]]; + let straddling: Polygon

= polygon![[ + (0.0, 50.0), + (0.0, 150.0), + (100.0, 150.0), + (100.0, 50.0), + (0.0, 50.0) + ]]; + assert_eq!( + vertices(&union_poly(&square, &straddling).unwrap()), + vec![ + (0.0, 150.0), + (100.0, 150.0), + (100.0, 50.0), + (100.0, 0.0), + (0.0, 0.0), + (0.0, 150.0), + ] + ); +} + +/// Two squares meeting at a single corner. Both output rings begin at that +/// corner, so the node they start at cannot separate them — Boost's `iterate` +/// tries operation 0 before operation 1 at a turn, which puts the lobe traced +/// along the *first* operand first. +#[test] +fn lobes_meeting_at_a_corner_are_ordered_by_operand() { + let lower: Polygon

= polygon![[ + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0), + (100.0, 0.0), + (0.0, 0.0) + ]]; + let upper: Polygon

= polygon![[ + (100.0, 100.0), + (100.0, 200.0), + (200.0, 200.0), + (200.0, 100.0), + (100.0, 100.0) + ]]; + let out = union_poly(&lower, &upper).unwrap(); + let rings: Vec> = out + .polygons() + .map(|pg| { + pg.exterior() + .points() + .map(|p| (p.get::<0>(), p.get::<1>())) + .collect() + }) + .collect(); + assert_eq!( + rings, + vec![ + vec![ + (100.0, 100.0), + (100.0, 0.0), + (0.0, 0.0), + (0.0, 100.0), + (100.0, 100.0) + ], + vec![ + (100.0, 100.0), + (100.0, 200.0), + (200.0, 200.0), + (200.0, 100.0), + (100.0, 100.0) + ], + ] + ); +} + +/// Two convex polygons crossing twice, where each crossing sits on a +/// *different* segment of the first operand but both sit in the same monotone +/// run of it. +/// +/// `get_turns` partitions each operand into sections — runs of segments +/// heading the same way in both dimensions — and walks the section pairs, so +/// two turns in one section of the first operand are ordered by the section of +/// the second, not by the first's segment index. Ordering by segment alone +/// starts this ring at the other crossing. +/// +/// The crossing coordinates are irrational, so the start is checked by +/// proximity rather than pinned digit for digit. +#[test] +fn turns_in_one_section_are_ordered_by_the_second_operands_section() { + let nine: Polygon

= polygon![[ + (181.0, 100.0), + (157.0, 43.0), + (100.0, 19.0), + (43.0, 43.0), + (19.0, 100.0), + (43.0, 157.0), + (100.0, 181.0), + (157.0, 157.0), + (181.0, 100.0) + ]]; + let ten: Polygon

= polygon![[ + (200.0, 4.0), + (188.0, -26.0), + (160.0, -43.0), + (129.0, -37.0), + (107.0, -12.0), + (107.0, 20.0), + (128.0, 45.0), + (160.0, 51.0), + (188.0, 34.0), + (200.0, 4.0) + ]]; + let start = vertices(&union_poly(&nine, &ten).unwrap())[0]; + // C++ Boost 1.83 begins here; ordering by segment would begin at the other + // crossing, near (160.293, 50.822). + assert!( + (start.0 - 109.530_944_625_407_16).abs() < 1e-9 + && (start.1 - 23.013_029_315_960_91).abs() < 1e-9, + "ring starts at {start:?}" + ); +} From e425fee69f38bc150c9d0cd6c33ea766d01dde23 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 18:24:46 -0600 Subject: [PATCH 14/20] fix(overlay): read the second operand backwards for a difference MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `difference` does not dispatch the overlay the way `union_` and `intersection` do: it passes `Reverse2 = true`, so `sectionalize` reads the second operand through a reversed view and every section and segment index `get_turns` sees for that operand counts from the other end of the ring. Since those indices are what order the turns, they decide which turn a result ring starts at — and a bite taken out of a single edge ties on everything the first operand can say, leaving the reversal as the whole decision. Ten of the 312 valid pairs started their ring at the wrong end of such a bite. The direction reaches nothing else here. The arrangement orients each edge by which side the result lies on, and the point-in-polygon test that classifies them counts crossings, so both are indifferent to how a ring was stored. The remaining three were the order of the polygons, not their contents. `add_rings` walks the selected rings in `ring_identifier` order, which puts every ring copied whole from an operand ahead of the traversed ones and orders the traversed ones by when `traverse` started them — that is, by where `get_turns` put the turn each one started from. Ordering by the node a ring starts at got close because nodes are created along the first operand's boundary, but it cannot separate two rings whose start turns differ only in the second operand. `TurnOrder` is now the key for both choices, the ring's start turn and the ring's place in the output, instead of the comparison being written out once inline and approximated once by node index. Measured over 312 valid polygon pairs against C++ Boost 1.83, compared exactly, ring start and ring order included: before after difference 36 differ, 13 structurally 25 differ, 0 structurally union 27 differ, 0 structurally unchanged intersection 27 differ, 0 structurally unchanged Every remaining difference across the three operations sits on an irrational crossing and differs in its last bits alone, worst 2.84e-14. The generator's self-intersecting inputs are unaffected either way, at 75, 94 and 161 of 600. Gate: 1614 tests pass, fmt and clippy clean. tilemaker's `scripts/verify.sh`: 16 passed, 0 failed. --- .../geometry-overlay/src/operation/areal.rs | 129 ++++++++++++++---- .../geometry-overlay/tests/overlay_parity.rs | 90 ++++++++++++ 2 files changed, 189 insertions(+), 30 deletions(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 38367e4..7a2131f 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -12,6 +12,7 @@ )] use alloc::vec::Vec; +use core::cmp::Ordering; use geometry_coords::{ CoordinateScalar, @@ -46,6 +47,18 @@ impl ArealOp { Self::SymDifference => first != second, } } + + /// Whether the second operand's rings are walked backwards. + /// + /// C++: `difference` dispatches the overlay with `Reverse2 = true`, so + /// `sectionalize` reads that operand through a reversed view and every + /// section and segment index it hands to `get_turns` counts from the other + /// end of the ring. That is what orders the turns, so it decides which one + /// a result ring starts at. Nothing else here depends on the direction: + /// the arrangement reorients each edge by which side the result lies on. + fn walks_second_operand_backwards(self) -> bool { + matches!(self, Self::Difference) + } } #[derive(Debug, Clone, Copy)] @@ -205,6 +218,60 @@ impl Edge { } } +/// Where a turn sits in `get_turns`' collection order. +/// +/// C++: `get_turns` walks the first operand's sections in the outer loop and +/// the second operand's in the inner, and inside a pair of sections the +/// segments in order — so a turn's place in `m_turns` is that nesting, and +/// `traverse` starts its rings in `m_turns` order. +#[derive(Clone, Copy)] +struct TurnOrder { + sections: [usize; 2], + arrivals: [usize; 2], + offset: f64, +} + +impl TurnOrder { + fn of

(node: &Node

) -> Self { + Self { + sections: node.section, + arrivals: node.arrival, + offset: node.offset[0], + } + } + + fn compare(&self, other: &Self) -> Ordering { + self.sections + .cmp(&other.sections) + .then_with(|| self.arrivals.cmp(&other.arrivals)) + .then_with(|| self.offset.total_cmp(&other.offset)) + } +} + +/// Where a finished ring falls in the output. +/// +/// C++: `add_rings` walks the selected rings in `ring_identifier` order. A +/// ring copied whole from an operand carries that operand's own identifier, so +/// every one of those precedes the traversed rings; a traversed ring is +/// identified by when `traverse` started it, which is where `get_turns` put +/// the turn it started from. +struct RingStart { + traversed: bool, + turn: TurnOrder, + second_operand: bool, + node: usize, +} + +impl RingStart { + fn compare(&self, other: &Self) -> Ordering { + self.traversed + .cmp(&other.traversed) + .then_with(|| self.turn.compare(&other.turn)) + .then_with(|| self.second_operand.cmp(&other.second_operand)) + .then_with(|| self.node.cmp(&other.node)) + } +} + /// Execute a polygon Boolean operation through a split-edge arrangement. pub(crate) fn overlay( first: &G1, @@ -221,8 +288,8 @@ where overlay_arrangement( &Shape::from_polygon(first), &Shape::from_polygon(second), - polygon_segments(first), - polygon_segments(second), + polygon_segments(first, false), + polygon_segments(second, operation.walks_second_operand_backwards()), operation, ) } @@ -248,8 +315,8 @@ where overlay_arrangement( &Shape::from_multi_polygon(first), &Shape::from_multi_polygon(second), - multi_polygon_segments(first), - multi_polygon_segments(second), + multi_polygon_segments(first, false), + multi_polygon_segments(second, operation.walks_second_operand_backwards()), operation, ) } @@ -384,7 +451,7 @@ where coordinates } -fn multi_polygon_segments(multi_polygon: &G) -> Vec> +fn multi_polygon_segments(multi_polygon: &G, backwards: bool) -> Vec> where G: MultiPolygonTrait, P: Point + Copy, @@ -393,15 +460,15 @@ where let mut segments = Vec::new(); let mut sections = Sectionizer::new(0); for polygon in multi_polygon.polygons() { - append_ring_segments(polygon.exterior(), &mut segments, &mut sections); + append_ring_segments(polygon.exterior(), &mut segments, &mut sections, backwards); for ring in polygon.interiors() { - append_ring_segments(ring, &mut segments, &mut sections); + append_ring_segments(ring, &mut segments, &mut sections, backwards); } } segments } -fn polygon_segments(polygon: &G) -> Vec> +fn polygon_segments(polygon: &G, backwards: bool) -> Vec> where G: PolygonTrait, P: Point + Copy, @@ -409,9 +476,9 @@ where { let mut segments = Vec::new(); let mut sections = Sectionizer::new(0); - append_ring_segments(polygon.exterior(), &mut segments, &mut sections); + append_ring_segments(polygon.exterior(), &mut segments, &mut sections, backwards); for ring in polygon.interiors() { - append_ring_segments(ring, &mut segments, &mut sections); + append_ring_segments(ring, &mut segments, &mut sections, backwards); } segments } @@ -483,12 +550,19 @@ fn append_ring_segments( ring: &R, output: &mut Vec>, sections: &mut Sectionizer, + backwards: bool, ) where R: RingTrait, P: Point + Copy, P::Scalar: Into, { - let points: Vec

= ring.points().copied().collect(); + let mut points: Vec

= ring.points().copied().collect(); + if backwards { + // C++: `reversible_view`, which reverses the closed ring — so a ring + // that stores its closing point still starts and ends on it, and one + // that does not still closes back to its own first vertex. + points.reverse(); + } if points.len() < 2 { return; } @@ -584,10 +658,9 @@ where P::Scalar: Into, { let mut used = alloc::vec![false; edges.len()]; - // Each ring is kept with the node it starts at, so the whole set can be - // put back into source order below, and with whether the traversal - // assembled it from turns rather than copying it whole from one operand. - let mut rings: Vec<(usize, bool, Ring

, bool)> = Vec::new(); + // Each ring is kept with where it starts, so the whole set can be put back + // into source order below. + let mut rings: Vec<(RingStart, Ring

)> = Vec::new(); for seed in 0..edges.len() { if used[seed] { continue; @@ -632,12 +705,11 @@ where // Which ring comes out first is observable — it decides the order of the // polygons in the result — and Boost's is not the order the seeds happened - // to fall in. It traces from turns in source order, so ordering the rings - // by the node each one starts at reproduces it. - rings.sort_by_key(|&(start, second_operand, _, _)| (start, second_operand)); + // to fall in. + rings.sort_by(|(left, _), (right, _)| left.compare(right)); Ok(rings .into_iter() - .map(|(_, _, ring, traced)| (ring, traced)) + .map(|(start, ring)| (ring, start.traversed)) .collect()) } @@ -687,7 +759,7 @@ where /// because a consumer that simplifies the ring afterwards will pin its first /// vertex and the choice reaches the output. fn push_ring

( - rings: &mut Vec<(usize, bool, Ring

, bool)>, + rings: &mut Vec<(RingStart, Ring

)>, nodes: &[Node

], node_indices: &[usize], carried: &[[bool; 2]], @@ -718,13 +790,7 @@ fn push_ring

( .enumerate() .filter(|&(_, index)| nodes[index].is_turn) .min_by(|&(_, left), &(_, right)| { - let (a, b) = (&nodes[left], &nodes[right]); - a.section[0] - .cmp(&b.section[0]) - .then(a.section[1].cmp(&b.section[1])) - .then(a.arrival[0].cmp(&b.arrival[0])) - .then(a.arrival[1].cmp(&b.arrival[1])) - .then(a.offset[0].total_cmp(&b.offset[0])) + TurnOrder::of(&nodes[left]).compare(&TurnOrder::of(&nodes[right])) }) .map(|(position, _)| position); // A ring with no turn was copied whole from one operand, and keeps that @@ -777,10 +843,13 @@ fn push_ring

( // traced along the first operand is emitted first. let leaves_along_first_operand = carried.get(first_turn).is_some_and(|edge| edge[0]); rings.push(( - cycle[first_turn], - !leaves_along_first_operand, + RingStart { + traversed: first_turn_by_arrival.is_some(), + turn: TurnOrder::of(&nodes[cycle[first_turn]]), + second_operand: !leaves_along_first_operand, + node: cycle[first_turn], + }, Ring::from_vec(points), - first_turn_by_arrival.is_some(), )); } diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 76c5de7..3f250a0 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -584,3 +584,93 @@ fn turns_in_one_section_are_ordered_by_the_second_operands_section() { "ring starts at {start:?}" ); } + +/// A pentagon with a smaller polygon cutting a bite out of one of its edges, +/// where both ends of the bite land on the *same* segment of the pentagon. +/// +/// C++: `difference` dispatches the overlay with `Reverse2 = true`, so +/// `sectionalize` reads the second operand backwards and the two turns come +/// out in the opposite order from the one their stored segments give. They tie +/// on everything the first operand can say, so that reversal is the whole +/// decision: read forwards, the ring starts at the other end of the bite. +#[test] +fn a_difference_reads_the_second_operand_backwards() { + let pentagon: Polygon

= polygon![[ + (182.0, 100.0), + (125.0, 23.0), + (34.0, 52.0), + (34.0, 148.0), + (125.0, 177.0), + (182.0, 100.0) + ]]; + let bite: Polygon

= polygon![[ + (135.0, 192.0), + (105.0, 153.0), + (60.0, 168.0), + (60.0, 216.0), + (105.0, 231.0), + (135.0, 192.0) + ]]; + let start = vertices(&difference(&pentagon, &bite).unwrap())[0]; + // C++ Boost 1.83 begins here; reading the second operand forwards would + // begin at the other end of the bite, near (122.962, 176.351). + assert!( + (start.0 - 77.966_292_134_831_46).abs() < 1e-9 + && (start.1 - 162.011_235_955_056_18).abs() < 1e-9, + "ring starts at {start:?}" + ); +} + +/// The same pentagon against a nonagon that clips three separate pieces off +/// it, so the result is three polygons and their order is what is under test. +/// +/// C++: `add_rings` emits the traversed rings in the order `traverse` started +/// them, which is where `get_turns` put each one's starting turn — not the +/// order the rings happened to be traced in. Two of these three start in the +/// same section of the first operand and are separated only by the second +/// operand's segment, so ordering by anything else swaps them. +#[test] +fn difference_pieces_come_out_in_the_order_their_turns_were_collected() { + let pentagon: Polygon

= polygon![[ + (182.0, 100.0), + (125.0, 23.0), + (34.0, 52.0), + (34.0, 148.0), + (125.0, 177.0), + (182.0, 100.0) + ]]; + let nonagon: Polygon

= polygon![[ + (161.0, 91.0), + (145.0, 49.0), + (106.0, 27.0), + (63.0, 34.0), + (33.0, 69.0), + (33.0, 113.0), + (62.0, 148.0), + (106.0, 155.0), + (145.0, 133.0), + (161.0, 91.0) + ]]; + let pieces = difference(&pentagon, &nonagon).unwrap(); + let sizes: Vec = pieces + .polygons() + .map(|pg| pg.exterior().points().count()) + .collect(); + // C++ Boost 1.83: the corner by (125, 23) first, then the body, then the + // sliver by (34, 52). Tracing order alone puts the body first. + assert_eq!(sizes, vec![4, 10, 4], "piece order"); + let corner: Vec<(f64, f64)> = pieces + .polygons() + .next() + .expect("three pieces") + .exterior() + .points() + .map(|p| (p.get::<0>(), p.get::<1>())) + .collect(); + assert!( + (corner[0].0 - 143.706_689_536_878_23).abs() < 1e-9 + && (corner[0].1 - 48.270_440_251_572_325).abs() < 1e-9, + "first piece starts at {:?}", + corner[0] + ); +} From 0fc83cabf0942f5d056ae11b45ce0f9d3927b3fb Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 18:57:24 -0600 Subject: [PATCH 15/20] fix(overlay): a ring no turn lands on keeps every vertex it has MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `traverse` never sees a ring that nothing crossed. `add_rings` copies it out of its operand through `convert_ring`, which appends nothing and drops nothing, so every vertex it was stored with survives — including one the outline runs straight through on its way back to the first. This arrangement traces such a ring like any other and closed it the way the traversal closes a *traced* ring, through `append_no_collinear`. That function exists to drop the point a new turn now runs straight past, and on the closing append it did exactly that to the last vertex of an untouched ring. Closing is now a plain append unless a turn actually started the ring. Found in tilemaker rather than in the parity set, which has no case where a whole ring survives an operation untouched: the dissolve it repairs polygons with finishes on `difference(outers, inners)`, and a repaired piece that nothing else touched came back a vertex short. Over the 767 `simplify` calls a liechtenstein run makes, the dissolve went from 7 outputs differing against C++ Boost 1.83 to 3, and the whole post-simplify tail — `correct`, the validity check and the repair — from 5 to 3; what is left there is ring order, not content. tilemaker's own gate has liechtenstein identical for the first time. The 312 valid polygon pairs are unchanged at zero structural differences for union, intersection and difference alike. Gate: 1615 tests pass, fmt and clippy clean. --- .../geometry-overlay/src/operation/areal.rs | 11 ++++++++- .../geometry-overlay/tests/overlay_parity.rs | 23 +++++++++++++++++++ 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 7a2131f..851f93c 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -831,8 +831,17 @@ fn push_ring

( // from, and that arrival is an append like any other — which is exactly // where the point before it goes, when the start carries the outline // straight on through it. + // + // Only a *traced* ring closes that way. One with no turn on it was never + // traversed at all: `add_rings` copies it out of its operand through + // `convert_ring`, which appends nothing and drops nothing, so its last + // vertex stays even where it continues the line straight into the first. if let Some(&first) = points.first() { - append_no_collinear(&mut points, first); + if first_turn_by_arrival.is_some() { + append_no_collinear(&mut points, first); + } else { + points.push(first); + } } // The ring is cleaned once it is in its final winding, not here: which // vertex `clean_closing_dups_and_spikes` leaves at the front depends on diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 3f250a0..c94692d 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -674,3 +674,26 @@ fn difference_pieces_come_out_in_the_order_their_turns_were_collected() { corner[0] ); } + +/// A polygon whose ring runs straight through its last vertex into its first, +/// differenced against something that does not touch it. +/// +/// C++: nothing traverses this ring — no turn lands on it, so `add_rings` +/// copies it out of its operand with `convert_ring`, which appends nothing and +/// drops nothing. Closing it the way the traversal closes a *traced* ring puts +/// the last vertex through `append_no_collinear`, which sees the straight run +/// into the first vertex and takes it off. +/// +/// This is what reached tilemaker: the dissolve it uses to repair a polygon +/// finishes with `difference(outers, inners)`, and a repaired piece that +/// nothing else touches came back a vertex short. +#[test] +fn an_untouched_ring_keeps_the_vertex_it_runs_straight_through() { + let sliver: Polygon

= polygon![[(3.0, 3.0), (2.0, 4.0), (3.0, 5.0), (3.0, 4.0), (3.0, 3.0)]]; + let elsewhere = square(20.0, 20.0, 4.0); + let kept = vertices(&difference(&sliver, &elsewhere).unwrap()); + assert_eq!( + kept, + vec![(3.0, 3.0), (2.0, 4.0), (3.0, 5.0), (3.0, 4.0), (3.0, 3.0)] + ); +} From 025ce640864c5d86ba2b64e5e3bf0adfbcd44ef6 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 20:04:06 -0600 Subject: [PATCH 16/20] fix(coords): two points within an epsilon are one point MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Boost's kernel does not compare coordinates exactly. `math::equals` calls two values equal when they differ by no more than one epsilon of the larger — or of `1`, so values near zero still have to agree absolutely — and `side_by_triangle` opens by calling three points collinear whenever any *two* of them are equal by that rule, before it computes any determinant at all (`side_by_triangle.hpp:150-164`). Every predicate built on the side test inherits it. `CoordinateScalar::tolerant_eq` is that rule, exact for an integer or a rational and relative for a float, and both the side predicate and `remove_spikes` now use it. The side predicate's own documentation already claimed the coincident-pair behaviour; it had the short-circuit only for points that were equal bit for bit. What it changes in practice is what counts as a spike. A hairline whose two ends are four last bits apart at a coordinate of 3540 is one point to Boost, which makes the ring a spike and collapses it to nothing; to an exact cross product it is a sliver with real area, and it survived into a monaco tile that the reference drew as nothing. Measured on tilemaker's own runs, `remove_spikes` now agrees with C++ Boost 1.83 on every call a whole extract makes: 0 of 767 under `config-example`, 0 of 5429 under OpenMapTiles, and 0 of 300 synthetic. The 312 valid polygon pairs are unchanged at zero structural differences for union, intersection and difference — the epsilon only fires where an exact determinant was already down in the last bits. Gate: 1618 tests pass, fmt and clippy clean. --- .../geometry-algorithm/src/remove_spikes.rs | 59 ++++++++++++++++++- crates/geometry-coords/src/rational.rs | 7 +++ crates/geometry-coords/src/scalar.rs | 32 ++++++++++ .../src/predicate/orientation.rs | 18 +++++- 4 files changed, 113 insertions(+), 3 deletions(-) diff --git a/crates/geometry-algorithm/src/remove_spikes.rs b/crates/geometry-algorithm/src/remove_spikes.rs index 37d7b44..f46de0a 100644 --- a/crates/geometry-algorithm/src/remove_spikes.rs +++ b/crates/geometry-algorithm/src/remove_spikes.rs @@ -58,7 +58,21 @@ fn is_spike_or_equal_2d(a: &P, b: &P, c: &P) -> bool { let cross = ux * vy - uy * vx; let dot = ux * vx + uy * vy; let zero = ::ZERO; - cross == zero && dot <= zero + // The collinearity half is Boost's `side_by_triangle`, which calls three + // points collinear whenever any *two* of them are equal by `math::equals` + // — a relative epsilon — before it looks at any determinant + // (`side_by_triangle.hpp:150-164`). A hairline whose two ends are a few + // last bits apart at a large coordinate is a spike to Boost and a genuine + // sliver to an exact cross product, which is how one survived into a tile + // that the reference drew as nothing. + let same = |ax: P::Scalar, ay: P::Scalar, bx: P::Scalar, by: P::Scalar| { + ax.tolerant_eq(bx) && ay.tolerant_eq(by) + }; + let collinear = cross == zero + || same(a.get::<0>(), a.get::<1>(), b.get::<0>(), b.get::<1>()) + || same(a.get::<0>(), a.get::<1>(), c.get::<0>(), c.get::<1>()) + || same(b.get::<0>(), b.get::<1>(), c.get::<0>(), c.get::<1>()); + collinear && dot <= zero } fn walk_spikes(pts: &mut alloc::vec::Vec

) { @@ -181,6 +195,49 @@ mod tests { type P = Point2D; + fn spike_ring(points: &[(f64, f64)]) -> Ring

{ + let mut ring = Ring::new(); + for &(x, y) in points { + ring.push(P::new(x, y)); + } + ring + } + + /// A hairline whose two ends are four last bits apart at a coordinate of + /// 3540, which is inside one epsilon of it. + /// + /// C++: `side_by_triangle` calls three points collinear when any two of + /// them are `math::equals` — a *relative* epsilon — before it computes any + /// determinant, so Boost sees a spike here and collapses the ring to a + /// single repeated point. An exact cross product sees a sliver with real + /// area and keeps it, which is how one survived into a monaco tile the + /// reference drew as nothing. + #[test] + fn a_hairline_within_an_epsilon_is_a_spike() { + let mut ring = spike_ring(&[ + (3_539.999_999_999_999_5, 482.199_999_999_999_76), + (3540.0, 482.199_999_999_999_8), + (3540.0, 479.0), + (3_539.999_999_999_999_5, 482.199_999_999_999_76), + ]); + remove_spikes(&mut ring); + assert_eq!(ring.0.len(), 2, "{:?}", ring.0); + } + + /// The same ring with its two ends far enough apart to be two points, + /// where the sliver has real area and stays. + #[test] + fn a_sliver_wider_than_an_epsilon_is_kept() { + let mut ring = spike_ring(&[ + (3_539.999_999_9, 482.199_999_9), + (3540.0, 482.2), + (3540.0, 479.0), + (3_539.999_999_9, 482.199_999_9), + ]); + remove_spikes(&mut ring); + assert_eq!(ring.0.len(), 4, "{:?}", ring.0); + } + #[test] fn out_and_back_spur_is_removed() { // (0,0) → (1,0) → (3,0) → (2,0): the tip (3,0) is a reversed diff --git a/crates/geometry-coords/src/rational.rs b/crates/geometry-coords/src/rational.rs index 6229634..9f9858d 100644 --- a/crates/geometry-coords/src/rational.rs +++ b/crates/geometry-coords/src/rational.rs @@ -425,4 +425,11 @@ impl CoordinateScalar for Rational { self } } + + /// Exact, like Boost's `math::equals` for any non-floating-point type: + /// a rational carries no rounding error for the epsilon to absorb. + #[inline] + fn tolerant_eq(self, other: Self) -> bool { + self == other + } } diff --git a/crates/geometry-coords/src/scalar.rs b/crates/geometry-coords/src/scalar.rs index bd79e5e..2f70abb 100644 --- a/crates/geometry-coords/src/scalar.rs +++ b/crates/geometry-coords/src/scalar.rs @@ -58,6 +58,21 @@ pub trait CoordinateScalar: /// (`boost/geometry/util/math.hpp`). #[must_use] fn abs(self) -> Self; + + /// Equality the way the kernel means it. + /// + /// Counterpart to `boost::geometry::math::equals` + /// (`boost/geometry/util/math.hpp`) under `equals_default_policy`: + /// exact for an integer, and for a float, equal when the difference is + /// within one epsilon of the larger magnitude — or of `1`, so that two + /// values near zero still have to agree to an absolute epsilon. + /// + /// This is not a convenience. Boost's side predicate calls three points + /// collinear when any *two* of them are equal by this rule, so a pair a + /// few last bits apart at a large coordinate is coincident to the whole + /// kernel, and every predicate built on the side test follows. + #[must_use] + fn tolerant_eq(self, other: Self) -> bool; } macro_rules! impl_scalar_float { @@ -69,6 +84,21 @@ macro_rules! impl_scalar_float { fn sqrt(self) -> Self { crate::math::sqrt(self) } #[inline] fn abs(self) -> Self { crate::math::abs(self) } + #[inline] + fn tolerant_eq(self, other: Self) -> bool { + if self == other { + return true; + } + if !self.is_finite() || !other.is_finite() { + return false; + } + // C++: `greatest(abs(a), abs(b), T(1))`, the factor + // `equals_default_policy` supplies. + let factor = crate::math::abs(self) + .max(crate::math::abs(other)) + .max(1.0); + crate::math::abs(self - other) <= <$t>::EPSILON * factor + } } )* }; } @@ -94,6 +124,8 @@ macro_rules! impl_scalar_int { } #[inline] fn abs(self) -> Self { <$t>::abs(self) } + #[inline] + fn tolerant_eq(self, other: Self) -> bool { self == other } } )* }; } diff --git a/crates/geometry-overlay/src/predicate/orientation.rs b/crates/geometry-overlay/src/predicate/orientation.rs index f70af60..ac28f70 100644 --- a/crates/geometry-overlay/src/predicate/orientation.rs +++ b/crates/geometry-overlay/src/predicate/orientation.rs @@ -22,8 +22,10 @@ //! Boost's //! `side_by_triangle` additionally treats any coincident pair among the //! three points as collinear -//! (`side_by_triangle.hpp:159-164`); this predicate does the same, -//! because a zero-length base line has no well-defined side. +//! (`side_by_triangle.hpp:150-164`); this predicate does the same, and +//! coincident means Boost's `math::equals` — a relative epsilon — not bitwise +//! equality, because a zero-length base line has no well-defined side and a +//! base line a few last bits long has none worth trusting. use geometry_coords::{CoordinateScalar, precise_math}; use geometry_trait::Point; @@ -91,6 +93,18 @@ where let rx = r.get::<0>(); let ry = r.get::<1>(); + // C++: `side_by_triangle` opens by calling the three points collinear if + // any two of them are `equals_point_point` — which is `math::equals` per + // coordinate, a *relative* epsilon (`side_by_triangle.hpp:150-164`). Two + // points a few last bits apart at a large coordinate are the same point to + // Boost, and the determinant below never gets to disagree. + let coincident = |ax: P::Scalar, ay: P::Scalar, bx: P::Scalar, by: P::Scalar| { + ax.tolerant_eq(bx) && ay.tolerant_eq(by) + }; + if coincident(px, py, qx, qy) || coincident(px, py, rx, ry) || coincident(qx, qy, rx, ry) { + return Sign::Collinear; + } + // Signed area of (p, q, r). Boost's `side_by_triangle::side_value` // computes the identical determinant // (`side_by_triangle.hpp` `side_value`). From ddfaf1059419be2f0c809f8d8a016e5e3ca00739 Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 20:04:19 -0600 Subject: [PATCH 17/20] fix(overlay): order the result rings the way `add_rings` does MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Two things decide which polygon comes out of an overlay first, and neither was right once an operand had more than one ring. **`partition`.** `get_turns` does not compare every section of one operand against every section of the other. It hands both lists to `geometry::partition`, which halves the plane recursively and visits only the pairs that can still meet — and only when *both* lists are longer than sixteen; under that it runs the plain nested loop. The turns land in `m_turns` in the order that walk finds them, `traverse` starts a ring at the first turn it has not used, and `add_rings` emits the rings in that order. Ordering the turns lexicographically by section pair is the same answer for a single small polygon and a different one as soon as the division kicks in, which is why this only ever showed up on multi-polygon operands. `section_partition` is that walk; `Node` carries the rank of its section pair instead of the pair itself. **Ring identity.** A ring no turn lands on is emitted under its own `ring_identifier` — source, then position within that operand — ahead of every traversed ring. The port ordered those by the lowest node on the cycle, which is not the same thing: this arrangement gives two rings that meet at a point the *same* node, so a piece could inherit its neighbour's position and overtake it. Edges now carry the section they came from, and since a section never spans a ring, the lowest section on a cycle names the ring it came out of. Measured against C++ Boost 1.83, compared exactly, ring order included: union / intersection / difference, 312 valid pairs 0 structural (unchanged) the same three over multi-polygon operands, 300 each 1 -> 0 structural the vendored dissolve, over the 767 `simplify` calls 3 -> 0 differ at all a liechtenstein run makes That last one is the whole of `correct(polygon, out, 1E-12)`, per-polygon simplification and repair included, and it is now exact on every call. Gate: 1618 tests pass, fmt and clippy clean. tilemaker's `scripts/verify.sh`: 16 passed, 0 failed, with liechtenstein identical. --- crates/geometry-overlay/src/operation.rs | 6 +- .../geometry-overlay/src/operation/areal.rs | 134 ++++++++-- .../src/operation/section_partition.rs | 251 ++++++++++++++++++ .../geometry-overlay/tests/overlay_parity.rs | 60 +++++ 4 files changed, 426 insertions(+), 25 deletions(-) create mode 100644 crates/geometry-overlay/src/operation/section_partition.rs diff --git a/crates/geometry-overlay/src/operation.rs b/crates/geometry-overlay/src/operation.rs index 39661aa..779102f 100644 --- a/crates/geometry-overlay/src/operation.rs +++ b/crates/geometry-overlay/src/operation.rs @@ -2,11 +2,13 @@ //! //! Mirrors the public drivers and areal machinery behind //! `boost/geometry/algorithms/{intersection,union,difference,sym_difference}.hpp`. -//! The `boolean` part owns the public entry contract and `areal` owns the -//! split-edge arrangement kernel; this root exposes only the aggregate surface. +//! The `boolean` part owns the public entry contract, `areal` owns the +//! split-edge arrangement kernel and `section_partition` the order `get_turns` +//! visits section pairs in; this root exposes only the aggregate surface. mod areal; mod boolean; +mod section_partition; pub use boolean::{ OverlayError, difference, difference_multi, intersection, intersection_multi, sym_difference, diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 851f93c..6c75dbf 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -27,6 +27,7 @@ use geometry_trait::{ use crate::assemble::assemble_traced; use crate::operation::OverlayError; +use crate::operation::section_partition::{Bounds, VisitRank}; use crate::predicate::segment_intersection::{SegmentIntersection, segment_intersection}; /// Boolean truth table applied to the two polygon interiors. @@ -191,6 +192,10 @@ struct Node

{ /// walks the pairs, so two turns in the same pair of sections keep their /// segment order while turns in different pairs do not. section: [usize; 2], + /// Where that pair of sections falls in the order `partition` visits them + /// — the whole of a turn's position in `m_turns`, above its segments. + /// `usize::MAX` until the arrangement knows both operands. + pair_rank: usize, /// How far along that segment. It orders two turns only once both segment /// indices have tied — the fraction must not outrank the second operand, /// or two turns sharing one edge come out in the wrong order. @@ -210,6 +215,11 @@ struct Edge { /// Reproducing that needs to know who carries each edge — see /// `drop_points_interior_to_a_walked_segment`. carried_by: [bool; 2], + /// Which section of each operand runs along it, `usize::MAX` for an + /// operand that does not. Sections never span a ring, so the lowest one on + /// a cycle names the ring the cycle came out of — which is the whole of + /// `ring_identifier` for a ring nothing crossed. + section: [usize; 2], } impl Edge { @@ -220,13 +230,14 @@ impl Edge { /// Where a turn sits in `get_turns`' collection order. /// -/// C++: `get_turns` walks the first operand's sections in the outer loop and -/// the second operand's in the inner, and inside a pair of sections the -/// segments in order — so a turn's place in `m_turns` is that nesting, and -/// `traverse` starts its rings in `m_turns` order. +/// C++: `partition` decides which pair of sections is looked at when, and +/// inside a pair `get_turns_in_sections` walks the first section's segments +/// outer and the second's inner — so a turn's place in `m_turns` is the pair's +/// rank and then that nesting, and `traverse` starts its rings in `m_turns` +/// order. #[derive(Clone, Copy)] struct TurnOrder { - sections: [usize; 2], + pair_rank: usize, arrivals: [usize; 2], offset: f64, } @@ -234,15 +245,15 @@ struct TurnOrder { impl TurnOrder { fn of

(node: &Node

) -> Self { Self { - sections: node.section, + pair_rank: node.pair_rank, arrivals: node.arrival, offset: node.offset[0], } } fn compare(&self, other: &Self) -> Ordering { - self.sections - .cmp(&other.sections) + self.pair_rank + .cmp(&other.pair_rank) .then_with(|| self.arrivals.cmp(&other.arrivals)) .then_with(|| self.offset.total_cmp(&other.offset)) } @@ -251,12 +262,15 @@ impl TurnOrder { /// Where a finished ring falls in the output. /// /// C++: `add_rings` walks the selected rings in `ring_identifier` order. A -/// ring copied whole from an operand carries that operand's own identifier, so -/// every one of those precedes the traversed rings; a traversed ring is +/// ring copied whole from an operand carries that operand's own identifier — +/// source, then position within it — so every one of those precedes the +/// traversed rings and they keep the operand's own order; a traversed ring is /// identified by when `traverse` started it, which is where `get_turns` put /// the turn it started from. struct RingStart { traversed: bool, + source: usize, + ring: usize, turn: TurnOrder, second_operand: bool, node: usize, @@ -264,11 +278,20 @@ struct RingStart { impl RingStart { fn compare(&self, other: &Self) -> Ordering { - self.traversed - .cmp(&other.traversed) - .then_with(|| self.turn.compare(&other.turn)) - .then_with(|| self.second_operand.cmp(&other.second_operand)) - .then_with(|| self.node.cmp(&other.node)) + self.traversed.cmp(&other.traversed).then_with(|| { + if self.traversed { + self.turn + .compare(&other.turn) + .then_with(|| self.second_operand.cmp(&other.second_operand)) + .then_with(|| self.node.cmp(&other.node)) + } else { + // Untouched rings are ordered by their identifier alone, which + // has nothing to do with where a turn fell. + self.source + .cmp(&other.source) + .then_with(|| self.ring.cmp(&other.ring)) + } + }) } } @@ -375,6 +398,16 @@ where snap_tolerance, ); + // C++: `partition` is handed the two section lists, and the order it + // visits their pairs in is the order the turns end up in. + let ranks = VisitRank::of( + §ion_bounds(&first_segments), + §ion_bounds(&second_segments), + ); + for node in &mut nodes { + node.pair_rank = ranks.rank(node.section[0], node.section[1]); + } + let sample_distance = (scale * 1e-8).max(snap_tolerance * 32.0); let mut boundary: Vec = Vec::new(); for candidate in candidates { @@ -417,6 +450,7 @@ where start: candidate.end, end: candidate.start, carried_by: candidate.carried_by, + section: candidate.section, } }; // The same stretch reaches here once per operand that carries it, so @@ -426,6 +460,8 @@ where Some(held) => { held.carried_by[0] |= edge.carried_by[0]; held.carried_by[1] |= edge.carried_by[1]; + held.section[0] = held.section[0].min(edge.section[0]); + held.section[1] = held.section[1].min(edge.section[1]); } None => boundary.push(edge), } @@ -483,6 +519,34 @@ where segments } +/// The box of every section, in section order. +/// +/// C++: each `section` carries the `bounding_box` `sectionalize` expanded over +/// its segments, and that box is the whole of what `partition` reasons about. +fn section_bounds

(segments: &[SourceSegment

]) -> Vec +where + P: Point + Copy, + P::Scalar: Into, +{ + let mut bounds: Vec = Vec::new(); + for segment in segments { + let box_ = Bounds::around( + [ + segment.start.get::<0>().into(), + segment.start.get::<1>().into(), + ], + [segment.end.get::<0>().into(), segment.end.get::<1>().into()], + ); + match bounds.get_mut(segment.section) { + Some(held) => held.expand(&box_), + // Sections are numbered from zero and in order, so a segment + // either extends the section being built or opens the next one. + None => bounds.push(box_), + } + } + bounds +} + /// C++: `sectionalize`'s cap, "defaults to 10, this seems to give the fastest /// results". const MAX_SEGMENTS_PER_SECTION: usize = 10; @@ -610,10 +674,13 @@ fn append_atomic_edges

( if start != end { let mut carried_by = [false; 2]; carried_by[operand] = true; + let mut sections = [usize::MAX; 2]; + sections[operand] = section; output.push(Edge { start, end, carried_by, + section: sections, }); } } @@ -641,6 +708,7 @@ where is_turn, arrival: [usize::MAX; 2], section: [usize::MAX; 2], + pair_rank: usize::MAX, offset: [0.0; 2], }); nodes.len() - 1 @@ -668,15 +736,15 @@ where let first = edges[seed].start; let mut edge_index = seed; let mut node_indices = alloc::vec![first]; - // `carried[i]` is who carries the edge from `node_indices[i]` to the - // node after it, so it always has one entry fewer. - let mut carried: alloc::vec::Vec<[bool; 2]> = alloc::vec::Vec::new(); + // `along[i]` is the edge from `node_indices[i]` to the node after it, + // so it always has one entry fewer. + let mut along: alloc::vec::Vec = alloc::vec::Vec::new(); for _ in 0..=edges.len() { debug_assert!(!used[edge_index]); used[edge_index] = true; let edge = edges[edge_index]; node_indices.push(edge.end); - carried.push(edge.carried_by); + along.push(edge); // A node the walk has already stood on closes a ring right here, // not only when the walk returns to the seed. Where two lobes of @@ -690,9 +758,9 @@ where .position(|&index| index == edge.end) { let loop_nodes = node_indices.split_off(start); - let loop_carried = carried.split_off(start); + let loop_along = along.split_off(start); node_indices.push(edge.end); - push_ring(&mut rings, nodes, &loop_nodes, &loop_carried, tolerance); + push_ring(&mut rings, nodes, &loop_nodes, &loop_along, tolerance); } if edge.end == first { @@ -762,7 +830,7 @@ fn push_ring

( rings: &mut Vec<(RingStart, Ring

)>, nodes: &[Node

], node_indices: &[usize], - carried: &[[bool; 2]], + along: &[Edge], tolerance: f64, ) where P: Point + Copy, @@ -850,10 +918,24 @@ fn push_ring

( // at a point, both lobes start there. Boost separates them by operand: // `iterate` tries operation 0 before operation 1 at a turn, so the lobe // traced along the first operand is emitted first. - let leaves_along_first_operand = carried.get(first_turn).is_some_and(|edge| edge[0]); + let leaves_along_first_operand = along.get(first_turn).is_some_and(|edge| edge.carried_by[0]); + // C++: a ring no turn lands on is emitted by `add_rings` under its own + // `ring_identifier` — source first, then where it sits in that operand. + // Its vertices say nothing about which: this arrangement gives two rings + // that meet at a point the same node, so the lowest node on a cycle can + // belong to a different ring altogether. The lowest section does not, + // because a section never spans a ring. + let source = usize::from(!along.iter().all(|edge| edge.carried_by[0])); + let ring = along + .iter() + .map(|edge| edge.section[source]) + .min() + .unwrap_or(usize::MAX); rings.push(( RingStart { traversed: first_turn_by_arrival.is_some(), + source, + ring, turn: TurnOrder::of(&nodes[cycle[first_turn]]), second_operand: !leaves_along_first_operand, node: cycle[first_turn], @@ -964,6 +1046,7 @@ mod tests { is_turn: false, arrival: [0, 0], section: [0, 0], + pair_rank: 0, offset: [0.0; 2], }, Node { @@ -972,6 +1055,7 @@ mod tests { is_turn: false, arrival: [1, 1], section: [1, 1], + pair_rank: 1, offset: [0.0; 2], }, Node { @@ -980,6 +1064,7 @@ mod tests { is_turn: false, arrival: [2, 2], section: [2, 2], + pair_rank: 2, offset: [0.0; 2], }, ]; @@ -988,16 +1073,19 @@ mod tests { start: 0, end: 1, carried_by: [true; 2], + section: [0; 2], }, Edge { start: 1, end: 2, carried_by: [true; 2], + section: [0; 2], }, Edge { start: 2, end: 0, carried_by: [true; 2], + section: [0; 2], }, ]; diff --git a/crates/geometry-overlay/src/operation/section_partition.rs b/crates/geometry-overlay/src/operation/section_partition.rs new file mode 100644 index 0000000..d85338b --- /dev/null +++ b/crates/geometry-overlay/src/operation/section_partition.rs @@ -0,0 +1,251 @@ +//! The order `get_turns` puts its turns in. +//! +//! `get_turns` does not compare every segment of one operand against every +//! segment of the other. It cuts each into **sections** — runs of consecutive +//! segments heading the same way — and hands the two lists to +//! `geometry::partition`, which recursively halves the plane and visits the +//! pairs of sections that can still meet. The turns land in `m_turns` in the +//! order that walk finds them, `traverse` starts a ring at the first turn it +//! has not used, and `add_rings` emits the rings in the order `traverse` +//! made them — so this order is the order of the polygons in the result. +//! +//! Under seventeen sections on either side `partition` skips the division and +//! runs the plain nested loop, which is why a single small polygon against +//! another comes out in plain section order and a multi-polygon does not. +//! +//! Mirrors `boost/geometry/algorithms/detail/partition.hpp` +//! (`partition::apply` and `partition_two_ranges::apply`). + +use alloc::vec::Vec; + +/// The `min_elements` `partition::apply` defaults to. Both collections must +/// be *larger* than this for the division to happen at all. +const MIN_ELEMENTS: usize = 16; + +/// C++: `recurse_ok`'s `level < 100`. +const MAX_LEVEL: usize = 100; + +/// An axis-aligned box, which is all `partition` knows about a section. +#[derive(Debug, Clone, Copy, PartialEq)] +pub(crate) struct Bounds { + pub min: [f64; 2], + pub max: [f64; 2], +} + +impl Bounds { + pub(crate) fn around(first: [f64; 2], second: [f64; 2]) -> Self { + Self { + min: [first[0].min(second[0]), first[1].min(second[1])], + max: [first[0].max(second[0]), first[1].max(second[1])], + } + } + + pub(crate) fn expand(&mut self, other: &Self) { + for axis in 0..2 { + self.min[axis] = self.min[axis].min(other.min[axis]); + self.max[axis] = self.max[axis].max(other.max[axis]); + } + } + + /// C++: `! disjoint_box_box`, which compares with `<` — so two boxes that + /// merely touch do overlap. + fn overlaps(&self, other: &Self) -> bool { + (0..2).all(|axis| self.max[axis] >= other.min[axis] && other.max[axis] >= self.min[axis]) + } + + /// C++: `divide_box`, splitting at the midpoint of one dimension. + #[expect( + clippy::manual_midpoint, + reason = "C++ divides the interval as `(mi + ma) / 2`, and where the two disagree the split lands on a different coordinate and the walk visits a different order" + )] + fn halves(&self, axis: usize) -> (Self, Self) { + let middle = (self.min[axis] + self.max[axis]) / 2.0; + let mut lower = *self; + let mut upper = *self; + lower.max[axis] = middle; + upper.min[axis] = middle; + (lower, upper) + } +} + +fn enclosing(bounds: &[Bounds], of: &[usize]) -> Bounds { + let mut total = bounds[of[0]]; + for &index in &of[1..] { + total.expand(&bounds[index]); + } + total +} + +/// C++: `divide_into_subsets`. A box that reaches into both halves is +/// "exceeding" and is matched against everything rather than descending. +fn split( + bounds: &[Bounds], + of: &[usize], + lower_box: &Bounds, + upper_box: &Bounds, +) -> [Vec; 3] { + let mut lower = Vec::new(); + let mut upper = Vec::new(); + let mut exceeding = Vec::new(); + for &index in of { + let in_lower = lower_box.overlaps(&bounds[index]); + let in_upper = upper_box.overlaps(&bounds[index]); + match (in_lower, in_upper) { + (true, true) => exceeding.push(index), + (true, false) => lower.push(index), + (false, true) => upper.push(index), + // C++: "Is nowhere", which the overlaps policy may allow. + (false, false) => {} + } + } + [lower, upper, exceeding] +} + +/// C++: `recurse_ok`, which also caps the depth at 100 levels. +fn deep_enough(of: &[usize], level: usize) -> bool { + of.len() >= MIN_ELEMENTS && level < MAX_LEVEL +} + +/// One run of the divide-and-conquer, over one pair of section lists. +struct Walk<'a> { + first: &'a [Bounds], + second: &'a [Bounds], + visited: Vec<(usize, usize)>, +} + +impl Walk<'_> { + /// C++: `handle_two`, the quadratic fallback — first collection outer. + fn pair_up(&mut self, ones: &[usize], twos: &[usize]) { + for &one in ones { + for &two in twos { + self.visited.push((one, two)); + } + } + } + + /// Descend if both sides are still worth dividing, else match them all. + fn narrow(&mut self, box_: &Bounds, ones: &[usize], twos: &[usize], level: usize, axis: usize) { + if deep_enough(ones, level) && deep_enough(twos, level) { + self.descend(box_, ones, twos, level + 1, 1 - axis); + } else { + self.pair_up(ones, twos); + } + } + + /// The straddlers of one side against both halves of the other. + /// + /// C++ decides this for the two halves *together* — all three lists have + /// to be big enough or none of them descends — so it is not two + /// independent `narrow` calls. + fn against_both_halves( + &mut self, + straddlers: &[usize], + halves: (&[usize], &[usize]), + straddlers_lead: bool, + level: usize, + axis: usize, + ) { + let (lower, upper) = halves; + let bounds = if straddlers_lead { + self.first + } else { + self.second + }; + if deep_enough(lower, level) && deep_enough(upper, level) && deep_enough(straddlers, level) + { + let box_ = enclosing(bounds, straddlers); + let (level, axis) = (level + 1, 1 - axis); + if straddlers_lead { + self.descend(&box_, straddlers, lower, level, axis); + self.descend(&box_, straddlers, upper, level, axis); + } else { + self.descend(&box_, lower, straddlers, level, axis); + self.descend(&box_, upper, straddlers, level, axis); + } + } else if straddlers_lead { + self.pair_up(straddlers, lower); + self.pair_up(straddlers, upper); + } else { + self.pair_up(lower, straddlers); + self.pair_up(upper, straddlers); + } + } + + /// C++: `partition_two_ranges::apply`. + fn descend( + &mut self, + box_: &Bounds, + ones: &[usize], + twos: &[usize], + level: usize, + axis: usize, + ) { + let (lower_box, upper_box) = box_.halves(axis); + let [lower1, upper1, exceeding1] = split(self.first, ones, &lower_box, &upper_box); + let [lower2, upper2, exceeding2] = split(self.second, twos, &lower_box, &upper_box); + + if !exceeding1.is_empty() { + let mut box_ = enclosing(self.first, &exceeding1); + if !exceeding2.is_empty() { + box_.expand(&enclosing(self.second, &exceeding2)); + } + self.narrow(&box_, &exceeding1, &exceeding2, level, axis); + self.against_both_halves(&exceeding1, (&lower2, &upper2), true, level, axis); + } + if !exceeding2.is_empty() { + self.against_both_halves(&exceeding2, (&lower1, &upper1), false, level, axis); + } + self.narrow(&lower_box, &lower1, &lower2, level, axis); + self.narrow(&upper_box, &upper1, &upper2, level, axis); + } +} + +/// Every pair of sections `get_turns` looks at, in the order it looks. +/// +/// C++: `geometry::partition::apply(sec1, sec2, visitor, …)`. +fn visit_order(first: &[Bounds], second: &[Bounds]) -> Vec<(usize, usize)> { + if first.is_empty() || second.is_empty() { + return Vec::new(); + } + let ones: Vec = (0..first.len()).collect(); + let twos: Vec = (0..second.len()).collect(); + let mut walk = Walk { + first, + second, + visited: Vec::new(), + }; + if first.len() > MIN_ELEMENTS && second.len() > MIN_ELEMENTS { + let mut total = enclosing(first, &ones); + total.expand(&enclosing(second, &twos)); + walk.descend(&total, &ones, &twos, 0, 0); + } else { + walk.pair_up(&ones, &twos); + } + walk.visited +} + +/// Where each pair of sections falls in that order. +/// +/// A pair the walk never reaches cannot hold a turn — the two sections' boxes +/// would have to be apart — so a lookup that misses sorts last rather than +/// claiming a position. +pub(crate) struct VisitRank(Vec<((usize, usize), usize)>); + +impl VisitRank { + pub(crate) fn of(first: &[Bounds], second: &[Bounds]) -> Self { + let mut ranked: Vec<((usize, usize), usize)> = visit_order(first, second) + .into_iter() + .enumerate() + .map(|(rank, pair)| (pair, rank)) + .collect(); + ranked.sort_unstable(); + ranked.dedup_by_key(|(pair, _)| *pair); + Self(ranked) + } + + pub(crate) fn rank(&self, first: usize, second: usize) -> usize { + self.0 + .binary_search_by_key(&(first, second), |&(pair, _)| pair) + .map_or(usize::MAX, |at| self.0[at].1) + } +} diff --git a/crates/geometry-overlay/tests/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index c94692d..126e6e5 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -697,3 +697,63 @@ fn an_untouched_ring_keeps_the_vertex_it_runs_straight_through() { vec![(3.0, 3.0), (2.0, 4.0), (3.0, 5.0), (3.0, 4.0), (3.0, 3.0)] ); } + +/// Three separate pieces, one of which a hole cuts into, and two of which no +/// turn lands on at all — with the third sharing a vertex with one of them. +/// +/// C++: `add_rings` emits the untouched rings first, under their own +/// `ring_identifier`, so they keep the operand's order; the traversed one +/// follows. The shared vertex is what makes this bite: an arrangement that +/// merges coincident points gives the second and third pieces a node in +/// common, so ordering the untouched rings by any vertex puts them the wrong +/// way round. Only the ring each cycle came out of says which is which. +/// +/// This is the shape the vendored dissolve hands to `difference` after it has +/// split a self-intersecting ring, which is where tilemaker met it. +#[test] +fn untouched_pieces_keep_their_operands_order() { + let pieces: MultiPolygon> = MultiPolygon(vec![ + polygon![[ + (3139.0, 3263.0), + (3104.0, 3325.0), + (3_103.231_759_656_652_2, 3_336.523_605_150_214_7), + (3139.0, 3263.0) + ]], + polygon![[ + (3103.0, 3344.0), + (3099.0, 3363.0), + (3103.0, 3346.0), + (3103.0, 3344.0) + ]], + polygon![[ + (3139.0, 3263.0), + (3165.0, 3210.0), + (3162.0, 3216.0), + (3139.0, 3263.0) + ]], + ]); + let bite: MultiPolygon> = MultiPolygon(vec![polygon![[ + (3_103.231_759_656_652_2, 3_336.523_605_150_214_7), + (3103.0, 3337.0), + (3103.0, 3340.0), + (3_103.231_759_656_652_2, 3_336.523_605_150_214_7) + ]]]); + let result = difference_multi(&pieces, &bite).unwrap(); + let starts: Vec<(f64, f64)> = result + .polygons() + .map(|pg| { + let first = pg.exterior().points().next().expect("a ring"); + (first.get::<0>(), first.get::<1>()) + }) + .collect(); + // C++ Boost 1.83: the two untouched pieces in operand order, then the one + // the bite ran through. + assert_eq!( + starts, + vec![ + (3103.0, 3344.0), + (3139.0, 3263.0), + (3_103.231_759_656_652_2, 3_336.523_605_150_214_7) + ] + ); +} From e60e4fe01c717a071c41bb5c23fd30b0711c643c Mon Sep 17 00:00:00 2001 From: Nick P Date: Thu, 3 Sep 2026 20:59:00 -0600 Subject: [PATCH 18/20] feat(buffer): a zero-width buffer of a polygon, which is not a no-op MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `buffer(g, 0)` is Boost's make-valid idiom and `repair_one_polygon` falls back on it when the dissolve gives up. This crate refused it outright, so an outline the dissolve could not fix came back as nothing. It is not a no-op and not a special case. `buffer_inserter` runs its whole pipeline: each side offsets — onto itself, at zero — into an **offsetted ring**, each join adds what it adds, the turns between those rings are found and checked against the original, and what survives is traversed. Two things fall out of the zero: * A **convex** corner adds nothing. `join_miter` returns early when the two offset points coincide, and at zero they always do. * A **concave** corner adds two `buffered_concave` pieces, each contributing the corner again — so that vertex is emitted three times. The winding a ring is stored with is what decides which corners are which: `closed_clockwise_view` is keyed on the ring *type*, not on how the ring is actually wound, so it is the identity for a type that is already closed and clockwise. Normalising the winding first inverts every corner, and closing the ring after reversing it rather than before starts it at the wrong vertex; both were tried and measured. Where the offsetted rings meet nothing, `discard_rings` drops nothing and the traversal has no work: the rings are the answer. That is the case the repair needs, and it is the case this implements. Where they do meet, what survives rests on the winding strategy's verdict for a turn point `get_turns` computed — a crossing part-way along two segments is decided by the last bits of a determinant. Reproducing the verdict without the arithmetic that produced the point was tried and measured (right 132 of 300, wrong 148), so it declines instead. Against C++ Boost 1.83 with the strategies tilemaker passes, over 300 synthetic polygons: **122 answered, 122 exact, 0 wrong, 178 declined** — up from 300 declined. And the outline that sent tilemaker here, an exterior collapsed to a single point around an interior that still encloses something, now comes back identical: twenty-one points, six corners tripled. `buffer_with(polygon, 0)` no longer returns `Unsupported` for a polygon whose rings meet nothing, which is a change to the public contract and to the test that pinned it. Gate: 1621 tests pass, fmt and clippy clean. tilemaker's buildings case is identical for the first time — 28,769 features, no differences. --- crates/geometry-overlay/src/buffer.rs | 40 +- crates/geometry-overlay/src/lib.rs | 1 + .../geometry-overlay/src/piece_collection.rs | 409 ++++++++++++++++++ .../geometry/tests/buffer_strategy_parity.rs | 19 +- 4 files changed, 467 insertions(+), 2 deletions(-) create mode 100644 crates/geometry-overlay/src/piece_collection.rs diff --git a/crates/geometry-overlay/src/buffer.rs b/crates/geometry-overlay/src/buffer.rs index 0e88b49..5f42184 100644 --- a/crates/geometry-overlay/src/buffer.rs +++ b/crates/geometry-overlay/src/buffer.rs @@ -338,6 +338,37 @@ where ) } +/// A polygon buffered at a distance of zero. +/// +/// C++: `buffer_inserter` builds an offsetted ring per input ring and then +/// finds the turns between them, discards those inside the original, and +/// traverses what is left. Where the offsetted rings do not meet each other +/// there are no turns, nothing is discarded and nothing is traversed, and the +/// rings themselves are the answer — which is the case +/// `repair_one_polygon` needs and the case this arm answers. +/// +/// A ring that does meet itself needs `check_turn_in_original` and the buffer +/// traversal, which are not ported; that asks for something this arm cannot +/// answer, and it says so rather than guessing. +fn zero_width_polygon_buffer( + polygon: &G, +) -> Result>, OverlayError> +where + G: PolygonTrait, + G::Point: PointMut + Default + Copy, + ::Scalar: CoordinateScalar + Into + FromF64, +{ + use crate::piece_collection::{ZeroWidthOutcome, zero_width_outcome, zero_width_rings}; + + let rings = zero_width_rings(polygon); + match zero_width_outcome(&rings) { + ZeroWidthOutcome::RingsStand => Ok(MultiPolygon( + rings.into_iter().map(Polygon::new).collect::>(), + )), + ZeroWidthOutcome::NeedsTraversal => Err(OverlayError::Unsupported), + } +} + /// Implements the point arm selected by `buffer_all` at /// `algorithms/detail/buffer/interface.hpp:269-273`. impl BufferStrategy for PointBuffer @@ -390,9 +421,16 @@ where let BufferDistanceStrategy::Symmetric(distance) = settings.distance else { return Err(OverlayError::Unsupported); }; - if !distance.is_finite() || distance == 0.0 { + if !distance.is_finite() { return Err(OverlayError::Unsupported); } + if distance == 0.0 { + // C++: a zero-width buffer is not a no-op and not a special case + // either — `buffer_inserter` runs its whole pipeline, and every + // side simply offsets onto itself. It is what `repair_one_polygon` + // falls back on, so it has to answer. + return zero_width_polygon_buffer(polygon); + } let Some(outer) = offset_ring(polygon.exterior(), distance, settings.join, true) else { return Ok(MultiPolygon(alloc::vec![])); }; diff --git a/crates/geometry-overlay/src/lib.rs b/crates/geometry-overlay/src/lib.rs index c99d8ae..d7da59c 100644 --- a/crates/geometry-overlay/src/lib.rs +++ b/crates/geometry-overlay/src/lib.rs @@ -36,6 +36,7 @@ pub mod buffer; pub mod line_intersection; pub mod merge; pub mod operation; +mod piece_collection; pub mod predicate; pub mod relate; pub mod surface_point; diff --git a/crates/geometry-overlay/src/piece_collection.rs b/crates/geometry-overlay/src/piece_collection.rs new file mode 100644 index 0000000..0abd8fd --- /dev/null +++ b/crates/geometry-overlay/src/piece_collection.rs @@ -0,0 +1,409 @@ +//! The rings Boost's buffer walks, at a distance of zero. +//! +//! `buffer_inserter` does not offset a ring and hand it back. It cuts the ring +//! into **pieces** — one per side, one per join — and appends each piece's +//! generated points to an **offsetted ring**; the turns between those rings are +//! then found, classified, and traversed. +//! +//! At a distance of zero the sides offset onto the segments themselves, so the +//! offsetted ring is the input ring with one addition: at every **concave** +//! corner Boost adds two `buffered_concave` pieces, each contributing the +//! corner again, so that vertex appears three times. A convex corner gets no +//! join at all, because `join_miter` returns early when the two offset points +//! coincide — which at zero they always do. +//! +//! That is not cosmetic. `repair_one_polygon` falls back on `buffer(0)` when +//! the dissolve gives up, and it is the only thing standing between an +//! outline with a collapsed exterior and nothing at all. +//! +//! Mirrors `boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp` +//! (`buffer_range::iterate`, `add_join`, `get_join_type`) and +//! `buffered_piece_collection.hpp` (`add_side_piece`, `add_range_to_piece`). + +use alloc::vec::Vec; + +use geometry_coords::CoordinateScalar; +use geometry_model::Ring; +use geometry_trait::{Point, PointMut, Polygon as PolygonTrait, Ring as RingTrait}; + +use crate::predicate::orientation::{Sign, orientation_2d}; + +/// What `get_join_type` makes of a corner. +/// +/// C++: `strategy::buffer::join_selector`, chosen from the side of the corner +/// (`side == -1` convex, `+1` concave) and, when the three points are +/// collinear, whether the third continues past the second or turns back. +#[derive(Debug, Clone, Copy, PartialEq, Eq)] +enum Join { + /// Nothing is added: at a distance of zero the two offset points coincide + /// and `join_miter` declines. + Convex, + /// Two `buffered_concave` pieces, each of which appends the corner again. + Concave, + /// Collinear and continuing: two consecutive sides, nothing between them. + Continue, + /// Collinear and turning back. An end cap, which a closed ring cannot ask + /// for at zero width, so nothing is added. + Spike, +} + +/// C++: `buffer_range::get_join_type`. +fn join_at

(before: &P, corner: &P, after: &P) -> Join +where + P: Point, + P::Scalar: CoordinateScalar + Into, +{ + match orientation_2d(before, corner, after) { + Sign::Negative => Join::Convex, + Sign::Positive => Join::Concave, + // C++: `same_direction`, which is `direction_code(...) == 1` — the + // perpendicular through the corner puts the third point beyond it. + Sign::Collinear => { + let dot = (corner.get::<0>().into() - before.get::<0>().into()) + * (after.get::<0>().into() - corner.get::<0>().into()) + + (corner.get::<1>().into() - before.get::<1>().into()) + * (after.get::<1>().into() - corner.get::<1>().into()); + if dot > 0.0 { + Join::Continue + } else { + Join::Spike + } + } + } +} + +fn same_point

(left: &P, right: &P) -> bool +where + P: Point, + P::Scalar: CoordinateScalar, +{ + left.get::<0>().tolerant_eq(right.get::<0>()) && left.get::<1>().tolerant_eq(right.get::<1>()) +} + +/// The ring as `closed_clockwise_view` presents it: distinct points, closed. +/// +/// The name is Boost's and so is the behaviour, which is less than the name +/// suggests: the view is keyed on the ring *type*'s declared order and +/// closure, not on how the ring is actually wound, so for a type that is +/// already closed and clockwise it is the identity. The winding a ring is +/// stored with is what the joins are judged against — which is how an interior +/// ring, stored the other way round by `correct`, gets the opposite corners. +fn clockwise_view(ring: &R) -> Vec

+where + R: RingTrait, + P: Point + Copy, + P::Scalar: CoordinateScalar + Into, +{ + let mut points: Vec

= Vec::new(); + for point in ring.points() { + if points.last().is_none_or(|last| !same_point(last, point)) { + points.push(*point); + } + } + while points.len() > 1 && same_point(&points[0], &points[points.len() - 1]) { + points.pop(); + } + if points.len() < 3 { + return points; + } + points.push(points[0]); + points +} + +/// The offsetted ring a zero-width buffer generates for one ring. +/// +/// C++: `buffer_range::iterate` over the sides, with `add_join` between each +/// consecutive pair and a closing join at the first vertex. +fn offsetted_ring

(closed: &[P]) -> Vec

+where + P: Point + Copy, + P::Scalar: CoordinateScalar + Into, +{ + let count = closed.len() - 1; + let mut out: Vec

= Vec::with_capacity(count * 2); + out.push(closed[0]); + for index in 0..count { + // The join at `closed[index]` is emitted before the side leaving it, + // and the first side has no join before it. + if index > 0 { + let before = closed[index - 1]; + let corner = closed[index]; + let after = closed[index + 1]; + if join_at(&before, &corner, &after) == Join::Concave { + out.push(corner); + out.push(corner); + } + } + out.push(closed[index + 1]); + } + // C++: `buffer_inserter_ring::iterate` adds a closing join at the first + // vertex once the sides are done. + if count >= 2 { + let before = closed[count - 1]; + let corner = closed[0]; + let after = closed[1]; + if join_at(&before, &corner, &after) == Join::Concave { + out.push(corner); + out.push(corner); + } + } + out +} + +/// Every offsetted ring a polygon generates at a distance of zero. +/// +/// A ring that collapses to a single distinct point becomes a point buffer, +/// which at zero width is a ring of coincident points enclosing nothing; it is +/// generated so the ring indices line up with Boost's, and dropped by the area +/// test downstream. `None` means the geometry asks for something this arm does +/// not implement. +pub(crate) fn zero_width_rings(polygon: &G) -> Vec> +where + G: PolygonTrait, + P: PointMut + Default + Copy, + P::Scalar: CoordinateScalar + Into, +{ + let mut rings = Vec::new(); + for ring in core::iter::once(polygon.exterior()).chain(polygon.interiors()) { + let closed = clockwise_view(ring); + if closed.len() < 4 { + // Fewer than three distinct points: nothing with area, and Boost's + // point buffer at zero width encloses nothing either. + continue; + } + rings.push(Ring::from_vec(offsetted_ring(&closed))); + } + rings +} + +/// One segment of an offsetted ring, tagged with the ring and its position. +struct Step { + ring: usize, + at: usize, + from: [f64; 2], + to: [f64; 2], +} + +/// Where two offsetted rings meet, which is the whole of what the rest of +/// Boost's pipeline reasons about. +/// +/// C++: `get_piece_turns`, which skips the pairs that are neighbours by +/// construction — a side and the join beside it always touch, and that is not +/// a turn. +fn turns

(rings: &[Ring

]) -> Vec<[f64; 2]> +where + P: Point + Copy, + P::Scalar: CoordinateScalar + Into, +{ + let mut steps: Vec = Vec::new(); + for (ring, points) in rings.iter().enumerate() { + let coords: Vec<[f64; 2]> = points + .points() + .map(|point| [point.get::<0>().into(), point.get::<1>().into()]) + .collect(); + // The concave joins put zero-length segments into the ring. They + // cannot meet anything, and leaving them in makes two segments that + // are neighbours look three apart — which is the difference between a + // ring that meets itself and one that merely turns a corner. + for pair in coords.windows(2) { + // Exact: these came out of the same vertex, so they are the same + // bits or they are a real segment. + #[expect(clippy::float_cmp, reason = "recognising a repeated vertex")] + let repeated = pair[0][0] == pair[1][0] && pair[0][1] == pair[1][1]; + if repeated { + continue; + } + steps.push(Step { + ring, + at: steps.iter().filter(|step| step.ring == ring).count(), + from: pair[0], + to: pair[1], + }); + } + } + let mut found = Vec::new(); + for (index, one) in steps.iter().enumerate() { + let length = steps.iter().filter(|step| step.ring == one.ring).count(); + for two in steps.iter().skip(index + 1) { + if one.ring == two.ring { + let apart = two.at - one.at; + if apart <= 1 || apart + 1 >= length { + continue; + } + } + if let Some(point) = meeting_point(one, two) { + found.push(point); + } + } + } + found +} + +fn meeting_point(one: &Step, two: &Step) -> Option<[f64; 2]> { + let cross = |a: [f64; 2], b: [f64; 2]| a[0] * b[1] - a[1] * b[0]; + let sub = |a: [f64; 2], b: [f64; 2]| [a[0] - b[0], a[1] - b[1]]; + let run = sub(one.to, one.from); + let other = sub(two.to, two.from); + let denominator = cross(run, other); + if denominator == 0.0 { + return None; + } + let offset = sub(two.from, one.from); + let along = cross(offset, other) / denominator; + let across = cross(offset, run) / denominator; + if !(0.0..=1.0).contains(&along) || !(0.0..=1.0).contains(&across) { + return None; + } + Some([one.from[0] + along * run[0], one.from[1] + along * run[1]]) +} + +/// Whether the offsetted rings are the answer on their own. +/// +/// C++: `discard_rings` drops every offsetted ring that has a turn, so what +/// comes out is the traversal's rings plus the offsetted rings that met +/// nothing. Where nothing met anything the traversal has no work and the rings +/// stand as they are — which is the case `repair_one_polygon` needs. +/// +/// Where they do meet, what survives rests on `check_turn_in_original`, and +/// that rests on the winding strategy's verdict for a turn point `get_turns` +/// computed: a crossing part-way along two segments is decided by the last bits +/// of a determinant. Reproducing the verdict without reproducing the arithmetic +/// that produced the point was tried and measured — right 132 times of 300 and +/// wrong 148 — so this declines instead of guessing. +pub(crate) enum ZeroWidthOutcome { + /// No ring meets another: the offsetted rings are the answer. + RingsStand, + /// Something met something, and what survives needs the buffer traversal. + NeedsTraversal, +} + +pub(crate) fn zero_width_outcome

(rings: &[Ring

]) -> ZeroWidthOutcome +where + P: Point + Copy, + P::Scalar: CoordinateScalar + Into, +{ + if turns(rings).is_empty() { + ZeroWidthOutcome::RingsStand + } else { + ZeroWidthOutcome::NeedsTraversal + } +} + +#[cfg(test)] +mod tests { + //! Checked against C++ Boost 1.83's `buffer` with + //! `distance_symmetric(0.0)`, `side_straight`, `join_miter`, + //! `end_flat` and `point_square` — the strategies tilemaker's + //! `repair_one_polygon` passes. + + use super::{ZeroWidthOutcome, zero_width_outcome, zero_width_rings}; + use geometry_cs::Cartesian; + use geometry_model::{Point2D, Polygon, Ring}; + use geometry_trait::{Point as _, Ring as _}; + + type P = Point2D; + + fn ring(points: &[(f64, f64)]) -> Ring

{ + let mut ring = Ring::new(); + for &(x, y) in points { + ring.push(P::new(x, y)); + } + ring + } + + fn points_of(ring: &Ring

) -> Vec<(f64, f64)> { + ring.points() + .map(|point| (point.get::<0>(), point.get::<1>())) + .collect() + } + + /// One concave corner in a ring that is otherwise convex. + /// + /// C++ returns the ring with that corner's vertex three times over — once + /// from the side and twice from the pair of `buffered_concave` pieces — + /// and every other vertex once. + #[test] + fn a_concave_corner_is_emitted_three_times() { + let polygon = Polygon::new(ring(&[ + (3186.0, 2762.0), + (3063.0, 2666.0), + (2953.0, 2536.0), + (2965.0, 2762.0), + (3023.0, 2920.0), + (3074.0, 2999.0), + (3186.0, 2762.0), + ])); + let rings = zero_width_rings(&polygon); + assert_eq!(rings.len(), 1); + assert_eq!( + points_of(&rings[0]), + vec![ + (3186.0, 2762.0), + (3063.0, 2666.0), + (3063.0, 2666.0), + (3063.0, 2666.0), + (2953.0, 2536.0), + (2965.0, 2762.0), + (3023.0, 2920.0), + (3074.0, 2999.0), + (3186.0, 2762.0), + ] + ); + assert!(matches!( + zero_width_outcome(&rings), + ZeroWidthOutcome::RingsStand + )); + } + + /// A ring that crosses itself needs the turns, the check against the + /// original and the traversal, none of which is ported. + #[test] + fn a_ring_that_meets_itself_is_declined() { + let polygon = Polygon::new(ring(&[ + (0.0, 0.0), + (10.0, 10.0), + (10.0, 0.0), + (0.0, 10.0), + (0.0, 0.0), + ])); + let rings = zero_width_rings(&polygon); + assert!(matches!( + zero_width_outcome(&rings), + ZeroWidthOutcome::NeedsTraversal + )); + } + + /// The outline that sent tilemaker here: an exterior collapsed to one + /// point, and an interior that still encloses something. C++ drops the + /// exterior — a point buffer of zero width holds nothing — and hands back + /// the interior, whose stored winding is what its corners are judged + /// against. + #[test] + fn a_collapsed_exterior_leaves_its_interior_standing() { + let polygon = Polygon::with_inners( + ring(&[ + (11.230_769_230_770_715, 4_095.000_000_000_000_5), + (11.230_769_230_770_715, 4_095.000_000_000_000_5), + ]), + vec![ring(&[ + (-11.0, 4091.0), + (-10.0, 4091.0), + (-10.0, 4_083.000_000_000_000_5), + (-2.0, 4087.0), + (-1.0, 4089.0), + (-4.0, 4092.0), + (-8.0, 4092.0), + (-9.0, 4093.0), + (-11.0, 4091.0), + ])], + ); + let rings = zero_width_rings(&polygon); + assert_eq!(rings.len(), 1, "the collapsed exterior encloses nothing"); + assert!(matches!( + zero_width_outcome(&rings), + ZeroWidthOutcome::RingsStand + )); + // C++ Boost 1.83: twenty-one points, six corners tripled. + assert_eq!(points_of(&rings[0]).len(), 21); + } +} diff --git a/crates/geometry/tests/buffer_strategy_parity.rs b/crates/geometry/tests/buffer_strategy_parity.rs index 77e4d39..a6d2a4e 100644 --- a/crates/geometry/tests/buffer_strategy_parity.rs +++ b/crates/geometry/tests/buffer_strategy_parity.rs @@ -286,6 +286,10 @@ fn polygon_buffer_handles_offset_topology_collapse() { /// `test/algorithms/buffer/buffer_with_strategies.cpp:88-106` — inapplicable /// distance strategies and degenerate inputs are rejected consistently. #[test] +#[expect( + clippy::too_many_lines, + reason = "one contract per geometry kind, read as a table" +)] fn public_buffer_error_and_empty_contract_is_consistent_across_kinds() { let asymmetric = BufferSettings { distance: BufferDistanceStrategy::Asymmetric { @@ -324,7 +328,20 @@ fn public_buffer_error_and_empty_contract_is_consistent_across_kinds() { buffer_with(&polygon, not_finite), Err(OverlayError::Unsupported) ); - assert_eq!(buffer_with(&polygon, zero), Err(OverlayError::Unsupported)); + // A zero-width buffer of a polygon is not an error and not a no-op: C++ + // runs the whole `buffer_inserter` pipeline, every side offsets onto + // itself, and this square comes back unchanged. It is what + // `repair_one_polygon` falls back on when the dissolve gives up. + assert_eq!( + buffer_with(&polygon, zero), + Ok(MultiPolygon(vec![polygon![[ + (0.0, 0.0), + (0.0, 2.0), + (2.0, 2.0), + (2.0, 0.0), + (0.0, 0.0) + ]]])) + ); assert!( buffer_convex_polygon(&polygon, 0.0, JoinStrategy::Miter) From 9b679bcec7228fcc1dccd39df7239442917dc81f Mon Sep 17 00:00:00 2001 From: Nick P Date: Fri, 4 Sep 2026 19:52:23 -0600 Subject: [PATCH 19/20] fix(build): keep the workspace clean under Rust 1.98 Stable moved to 1.98.1, which adds clippy::manual_midpoint, detects the unused trait imports, and rewords a trait-obligation diagnostic. None of this is a behaviour change; f64::midpoint matches (a + b) * 0.5 for in-range values. --- crates/geometry-algorithm/src/rhumb.rs | 1 - crates/geometry-overlay/src/operation/areal.rs | 4 ++-- crates/geometry-overlay/src/relate.rs | 16 ++++++++++++---- crates/geometry-rtree/src/bounds.rs | 4 ++-- crates/geometry-strategy/src/geographic/rhumb.rs | 1 - .../tests/ui/cartesian_only.stderr | 3 ++- 6 files changed, 18 insertions(+), 11 deletions(-) diff --git a/crates/geometry-algorithm/src/rhumb.rs b/crates/geometry-algorithm/src/rhumb.rs index fb39a83..020cfca 100644 --- a/crates/geometry-algorithm/src/rhumb.rs +++ b/crates/geometry-algorithm/src/rhumb.rs @@ -128,7 +128,6 @@ type PointOutput

= >::Output; mod tests { use geometry_cs::{Degree, Spherical}; use geometry_model::{Linestring, Point2D}; - use geometry_trait::Point as _; use super::*; diff --git a/crates/geometry-overlay/src/operation/areal.rs b/crates/geometry-overlay/src/operation/areal.rs index 6c75dbf..772b297 100644 --- a/crates/geometry-overlay/src/operation/areal.rs +++ b/crates/geometry-overlay/src/operation/areal.rs @@ -417,8 +417,8 @@ where let length = hypot(delta.0, delta.1); debug_assert!(length > snap_tolerance); let midpoint = Coordinate { - x: (start.x + end.x) * 0.5, - y: (start.y + end.y) * 0.5, + x: f64::midpoint(start.x, end.x), + y: f64::midpoint(start.y, end.y), }; let offset = sample_distance.min(length * 1e-4); let normal = (-delta.1 / length * offset, delta.0 / length * offset); diff --git a/crates/geometry-overlay/src/relate.rs b/crates/geometry-overlay/src/relate.rs index ec27fd5..b3fd902 100644 --- a/crates/geometry-overlay/src/relate.rs +++ b/crates/geometry-overlay/src/relate.rs @@ -882,7 +882,7 @@ where let sample = interpolate( first_segment.0, first_segment.1, - (interval[0] + interval[1]) * 0.5, + f64::midpoint(interval[0], interval[1]), ); let location = xy_location_linestring(sample, second); matrix.m[feature::INTERIOR][location.index()] = Dimension::Curve; @@ -896,7 +896,7 @@ where let sample = interpolate( second_segment.0, second_segment.1, - (interval[0] + interval[1]) * 0.5, + f64::midpoint(interval[0], interval[1]), ); let location = xy_location_linestring(sample, first); matrix.m[location.index()][feature::INTERIOR] = Dimension::Curve; @@ -1450,7 +1450,11 @@ fn record_segment_cells( let parameters = segment_parameters(segment, all_segments, &second.points); for interval in parameters.windows(2) { debug_assert!(interval[1] - interval[0] > f64::EPSILON); - let midpoint = interpolate(segment.0, segment.1, (interval[0] + interval[1]) * 0.5); + let midpoint = interpolate( + segment.0, + segment.1, + f64::midpoint(interval[0], interval[1]), + ); let first_location = topology_location(first, midpoint); let second_location = topology_location(second, midpoint); debug_assert_ne!(first_location, Location::Exterior); @@ -1534,7 +1538,11 @@ fn relate_topologies(first: &Topology, second: &Topology) -> Result f64::EPSILON); - let midpoint = interpolate(segment.0, segment.1, (interval[0] + interval[1]) * 0.5); + let midpoint = interpolate( + segment.0, + segment.1, + f64::midpoint(interval[0], interval[1]), + ); let first_location = topology_location(first, midpoint); let second_location = topology_location(second, midpoint); debug_assert_ne!(second_location, Location::Exterior); diff --git a/crates/geometry-rtree/src/bounds.rs b/crates/geometry-rtree/src/bounds.rs index 5368aec..be950e3 100644 --- a/crates/geometry-rtree/src/bounds.rs +++ b/crates/geometry-rtree/src/bounds.rs @@ -149,8 +149,8 @@ impl Bounds { #[must_use] pub fn center(&self) -> [f64; 2] { [ - (self.min[0] + self.max[0]) * 0.5, - (self.min[1] + self.max[1]) * 0.5, + f64::midpoint(self.min[0], self.max[0]), + f64::midpoint(self.min[1], self.max[1]), ] } } diff --git a/crates/geometry-strategy/src/geographic/rhumb.rs b/crates/geometry-strategy/src/geographic/rhumb.rs index 9e3f01d..019906a 100644 --- a/crates/geometry-strategy/src/geographic/rhumb.rs +++ b/crates/geometry-strategy/src/geographic/rhumb.rs @@ -192,7 +192,6 @@ fn reflect_latitude(latitude: f64) -> f64 { mod tests { use geometry_cs::{Degree, Spherical}; use geometry_model::{Linestring, Point2D}; - use geometry_trait::Point as _; use super::*; diff --git a/crates/geometry-strategy/tests/ui/cartesian_only.stderr b/crates/geometry-strategy/tests/ui/cartesian_only.stderr index d8c4380..712f962 100644 --- a/crates/geometry-strategy/tests/ui/cartesian_only.stderr +++ b/crates/geometry-strategy/tests/ui/cartesian_only.stderr @@ -9,7 +9,7 @@ error[E0277]: coordinate-system family mismatch: `GeographicFamily` is not the s = help: the trait `geometry_tag::same_as::SameAs` is not implemented for `GeographicFamily` = note: if `GeographicFamily` is `GeographicFamily` or `SphericalFamily`, you probably picked a Cartesian-only strategy (e.g. `Pythagoras`) by accident — wrap your point in `geometry_adapt::WithCs<_, Geographic>` / `WithCs<_, Spherical>` or pick a CS-appropriate strategy (`Haversine`, `Andoyer`, `Vincenty`) = note: see the silent-Cartesian discussion in the rust-port proposal §3.7 and §8 -help: the trait `DistanceStrategy` is implemented for `Pythagoras` +help: the trait `DistanceStrategy` is conditionally implemented for `Pythagoras` --> src/cartesian/distance_pythagoras.rs | | / impl DistanceStrategy for Pythagoras @@ -17,6 +17,7 @@ help: the trait `DistanceStrategy` is implemented for `Pythagoras` | | P1: Point, | | P2: Point, | | ::Family: SameAs, + | | ----------------------- unsatisfied requirement introduced here: `<, Geographic> as geometry_trait::point::Point>::Cs as CoordinateSystem>::Family: geometry_tag::same_as::SameAs` | | ::Family: SameAs, | |__________________________________________________________________^ = note: required for `Pythagoras` to implement `DistanceStrategy, Geographic>, _>` From 20f9d8e3c32afc47bad62d80ad2a6848ee900628 Mon Sep 17 00:00:00 2001 From: Nick P Date: Fri, 4 Sep 2026 19:56:25 -0600 Subject: [PATCH 20/20] docs(readme): regenerate the feature table for the multi-polygon operations The multi-polygon Boolean entry points were added without refreshing the generated table, which CI checks against the pub-use tags. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index ee8a42e..86192c7 100644 --- a/README.md +++ b/README.md @@ -370,7 +370,7 @@ plus a `_with` companion that takes an explicit strategy. | `within_dyn` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/algorithm/fn.within_dyn.html) | | **Boolean operations** — Overlay and offset of areal geometries ||| | `buffer` / `buffer_convex_polygon` / `buffer_point` / `buffer_with` / `buffer_with_strategy` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/overlay/fn.buffer.html) | -| `difference` / `intersection` / `sym_difference` / `union` / `union_poly` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/overlay/fn.difference.html) | +| `difference` / `difference_multi` / `intersection` / `intersection_multi` / `sym_difference` / `sym_difference_multi` / `union` / `union_multi` / `union_poly` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/overlay/fn.difference.html) | | `line_intersection` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/overlay/fn.line_intersection.html) | | `point_on_surface` | ✅ | [→](https://docs.rs/boost_geometry/latest/boost_geometry/overlay/fn.point_on_surface.html) | | **Construction & transformation** — Derive a new geometry from an existing one |||