diff --git a/crates/geometry-adapt-geo-types/CHANGELOG.md b/crates/geometry-adapt-geo-types/CHANGELOG.md index c4e52f0..082f8fe 100644 --- a/crates/geometry-adapt-geo-types/CHANGELOG.md +++ b/crates/geometry-adapt-geo-types/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-geo-types-v0.0.8...geometry-adapt-geo-types-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-geo-types-v0.0.7...geometry-adapt-geo-types-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-adapt-geo-types/Cargo.toml b/crates/geometry-adapt-geo-types/Cargo.toml index c602d9b..06bf5cb 100644 --- a/crates/geometry-adapt-geo-types/Cargo.toml +++ b/crates/geometry-adapt-geo-types/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-adapt-geo-types" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -12,10 +12,10 @@ readme = "README.md" [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords" } -geometry-cs = { version = "0.0.8", path = "../geometry-cs" } -geometry-trait = { version = "0.0.8", path = "../geometry-trait" } -geometry-model = { version = "0.0.8", path = "../geometry-model" } +geometry-coords = { version = "0.0.9", path = "../geometry-coords" } +geometry-cs = { version = "0.0.9", path = "../geometry-cs" } +geometry-trait = { version = "0.0.9", path = "../geometry-trait" } +geometry-model = { version = "0.0.9", path = "../geometry-model" } geo-types = "0.7" [dev-dependencies] diff --git a/crates/geometry-adapt-nalgebra/CHANGELOG.md b/crates/geometry-adapt-nalgebra/CHANGELOG.md index 6def585..e6642e9 100644 --- a/crates/geometry-adapt-nalgebra/CHANGELOG.md +++ b/crates/geometry-adapt-nalgebra/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-nalgebra-v0.0.8...geometry-adapt-nalgebra-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-nalgebra-v0.0.7...geometry-adapt-nalgebra-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-adapt-nalgebra/Cargo.toml b/crates/geometry-adapt-nalgebra/Cargo.toml index 3d6da67..3e67f65 100644 --- a/crates/geometry-adapt-nalgebra/Cargo.toml +++ b/crates/geometry-adapt-nalgebra/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-adapt-nalgebra" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -12,9 +12,9 @@ readme = "README.md" [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords" } -geometry-cs = { version = "0.0.8", path = "../geometry-cs" } -geometry-trait = { version = "0.0.8", path = "../geometry-trait" } +geometry-coords = { version = "0.0.9", path = "../geometry-coords" } +geometry-cs = { version = "0.0.9", path = "../geometry-cs" } +geometry-trait = { version = "0.0.9", path = "../geometry-trait" } nalgebra = "0.32" [dev-dependencies] diff --git a/crates/geometry-adapt/CHANGELOG.md b/crates/geometry-adapt/CHANGELOG.md index 4b48a68..dd1fab3 100644 --- a/crates/geometry-adapt/CHANGELOG.md +++ b/crates/geometry-adapt/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-v0.0.8...geometry-adapt-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-adapt-v0.0.7...geometry-adapt-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-adapt/Cargo.toml b/crates/geometry-adapt/Cargo.toml index 28c462a..5b491a9 100644 --- a/crates/geometry-adapt/Cargo.toml +++ b/crates/geometry-adapt/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-adapt" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -24,9 +24,9 @@ libm = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } [dev-dependencies] geometry-algorithm = { path = "../geometry-algorithm" } diff --git a/crates/geometry-algorithm/CHANGELOG.md b/crates/geometry-algorithm/CHANGELOG.md index 136bf16..1f3be48 100644 --- a/crates/geometry-algorithm/CHANGELOG.md +++ b/crates/geometry-algorithm/CHANGELOG.md @@ -7,6 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-algorithm-v0.0.8...geometry-algorithm-v0.0.9) - 2026-09-05 + +### Fixed + +- *(coords)* two points within an epsilon are one point +- *(algorithm)* collapse repeated vertices in remove_spikes + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-algorithm-v0.0.7...geometry-algorithm-v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry-algorithm/Cargo.toml b/crates/geometry-algorithm/Cargo.toml index c6058eb..be3f744 100644 --- a/crates/geometry-algorithm/Cargo.toml +++ b/crates/geometry-algorithm/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-algorithm" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -27,12 +27,12 @@ libm = [ ] [dependencies] -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } geometry-tag = { version = "0.0.7", path = "../geometry-tag", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-strategy = { version = "0.0.8", path = "../geometry-strategy", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-strategy = { version = "0.0.9", path = "../geometry-strategy", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } [dev-dependencies] geometry-adapt = { path = "../geometry-adapt" } diff --git a/crates/geometry-algorithm/src/remove_spikes.rs b/crates/geometry-algorithm/src/remove_spikes.rs index b9d8f52..f46de0a 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,21 @@ 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 + // 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

) { @@ -53,7 +81,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 +133,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,11 +190,54 @@ 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; + 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 @@ -306,4 +379,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-coords/CHANGELOG.md b/crates/geometry-coords/CHANGELOG.md index e04c19d..c355ce7 100644 --- a/crates/geometry-coords/CHANGELOG.md +++ b/crates/geometry-coords/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-coords-v0.0.8...geometry-coords-v0.0.9) - 2026-09-05 + +### Fixed + +- *(coords)* two points within an epsilon are one point + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-coords-v0.0.7...geometry-coords-v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry-coords/Cargo.toml b/crates/geometry-coords/Cargo.toml index 67184ca..5e1ab1f 100644 --- a/crates/geometry-coords/Cargo.toml +++ b/crates/geometry-coords/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-coords" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true 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-cs/CHANGELOG.md b/crates/geometry-cs/CHANGELOG.md index 8aee526..8d2944c 100644 --- a/crates/geometry-cs/CHANGELOG.md +++ b/crates/geometry-cs/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-cs-v0.0.8...geometry-cs-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-cs-v0.0.7...geometry-cs-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-cs/Cargo.toml b/crates/geometry-cs/Cargo.toml index f21df71..6885787 100644 --- a/crates/geometry-cs/Cargo.toml +++ b/crates/geometry-cs/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-cs" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -11,7 +11,7 @@ description = "Coordinate system types: Cartesian, Spherical, Geographic, and an readme = "README.md" [dependencies] -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } [features] default = ["std"] diff --git a/crates/geometry-io-geojson/CHANGELOG.md b/crates/geometry-io-geojson/CHANGELOG.md index 958af61..7c38c1d 100644 --- a/crates/geometry-io-geojson/CHANGELOG.md +++ b/crates/geometry-io-geojson/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-io-geojson-v0.0.8...geometry-io-geojson-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-io-geojson-v0.0.7...geometry-io-geojson-v0.0.8) - 2026-07-15 ### Fixed diff --git a/crates/geometry-io-geojson/Cargo.toml b/crates/geometry-io-geojson/Cargo.toml index 009afe3..7fc89c1 100644 --- a/crates/geometry-io-geojson/Cargo.toml +++ b/crates/geometry-io-geojson/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-io-geojson" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -21,10 +21,10 @@ std = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } itoa = "1" [lints] diff --git a/crates/geometry-io-svg/CHANGELOG.md b/crates/geometry-io-svg/CHANGELOG.md index ac62f72..5c0f769 100644 --- a/crates/geometry-io-svg/CHANGELOG.md +++ b/crates/geometry-io-svg/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-io-svg-v0.0.8...geometry-io-svg-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-io-svg-v0.0.7...geometry-io-svg-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-io-svg/Cargo.toml b/crates/geometry-io-svg/Cargo.toml index 5963159..f2fa38f 100644 --- a/crates/geometry-io-svg/Cargo.toml +++ b/crates/geometry-io-svg/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-io-svg" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -21,10 +21,10 @@ std = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } [lints] workspace = true diff --git a/crates/geometry-io-wkb/CHANGELOG.md b/crates/geometry-io-wkb/CHANGELOG.md index ad38ce7..83a182c 100644 --- a/crates/geometry-io-wkb/CHANGELOG.md +++ b/crates/geometry-io-wkb/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-io-wkb-v0.0.8...geometry-io-wkb-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-io-wkb-v0.0.7...geometry-io-wkb-v0.0.8) - 2026-07-15 ### Fixed diff --git a/crates/geometry-io-wkb/Cargo.toml b/crates/geometry-io-wkb/Cargo.toml index 6b5f6e6..8f6938c 100644 --- a/crates/geometry-io-wkb/Cargo.toml +++ b/crates/geometry-io-wkb/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-io-wkb" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -21,10 +21,10 @@ std = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } [lints] workspace = true diff --git a/crates/geometry-io-wkt/CHANGELOG.md b/crates/geometry-io-wkt/CHANGELOG.md index 4f4fd47..57e03a7 100644 --- a/crates/geometry-io-wkt/CHANGELOG.md +++ b/crates/geometry-io-wkt/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-io-wkt-v0.0.8...geometry-io-wkt-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-io-wkt-v0.0.7...geometry-io-wkt-v0.0.8) - 2026-07-15 ### Fixed diff --git a/crates/geometry-io-wkt/Cargo.toml b/crates/geometry-io-wkt/Cargo.toml index bcd4f63..7e749de 100644 --- a/crates/geometry-io-wkt/Cargo.toml +++ b/crates/geometry-io-wkt/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-io-wkt" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -21,10 +21,10 @@ std = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } itoa = "1" ryu = "1" diff --git a/crates/geometry-model/CHANGELOG.md b/crates/geometry-model/CHANGELOG.md index 368388f..e8d8097 100644 --- a/crates/geometry-model/CHANGELOG.md +++ b/crates/geometry-model/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-model-v0.0.8...geometry-model-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-model-v0.0.7...geometry-model-v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry-model/Cargo.toml b/crates/geometry-model/Cargo.toml index c4a1a2e..7b07dbc 100644 --- a/crates/geometry-model/Cargo.toml +++ b/crates/geometry-model/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-model" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -22,9 +22,9 @@ serde = ["dep:serde"] [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } serde = { version = "1", optional = true, default-features = false, features = ["derive", "alloc"] } [dev-dependencies] diff --git a/crates/geometry-overlay/CHANGELOG.md b/crates/geometry-overlay/CHANGELOG.md index da0768a..3ad222b 100644 --- a/crates/geometry-overlay/CHANGELOG.md +++ b/crates/geometry-overlay/CHANGELOG.md @@ -7,6 +7,31 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-overlay-v0.0.8...geometry-overlay-v0.0.9) - 2026-09-05 + +### Added + +- *(buffer)* a zero-width buffer of a polygon, which is not a no-op +- *(overlay)* multi-polygon operands for the Boolean operations + +### Fixed + +- *(overlay)* order the result rings the way `add_rings` does +- *(coords)* two points within an epsilon are one point +- *(overlay)* a ring no turn lands on keeps every vertex it has +- *(overlay)* read the second operand backwards for a difference +- *(overlay)* append a turn the way the traversal does, and order turns by section +- *(overlay)* do not emit a point the walked operand runs straight past +- *(overlay)* clean the traversed ring, and order turns by both operands +- *(overlay)* attach a turn to the segment it terminates +- *(overlay)* emit rings, and their vertices, in Boost's order +- *(predicate)* compute a segment crossing parametrically, not by determinant +- *(overlay)* a hole sharing an edge with the exterior is a self-intersection +- *(overlay)* distinguish the two ways multi-polygon members can be wrong +- *(algorithm)* collapse repeated vertices in remove_spikes +- *(overlay)* report wrong orientation before self-intersection +- *(overlay)* split result lobes that meet at a single point + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-overlay-v0.0.7...geometry-overlay-v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry-overlay/Cargo.toml b/crates/geometry-overlay/Cargo.toml index 633a2dd..daf7449 100644 --- a/crates/geometry-overlay/Cargo.toml +++ b/crates/geometry-overlay/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-overlay" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -31,12 +31,12 @@ libm = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } -geometry-strategy = { version = "0.0.8", path = "../geometry-strategy", default-features = false } -geometry-algorithm = { version = "0.0.8", path = "../geometry-algorithm", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } +geometry-strategy = { version = "0.0.9", path = "../geometry-strategy", default-features = false } +geometry-algorithm = { version = "0.0.9", path = "../geometry-algorithm", default-features = false } [lints] workspace = true 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/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 aa11aca..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; @@ -54,7 +55,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..779102f 100644 --- a/crates/geometry-overlay/src/operation.rs +++ b/crates/geometry-overlay/src/operation.rs @@ -2,10 +2,15 @@ //! //! 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, 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 c6e6123..6c75dbf 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, @@ -20,10 +21,13 @@ 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::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. @@ -44,6 +48,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)] @@ -82,6 +98,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() @@ -92,7 +122,13 @@ 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)>, + /// Which monotone run of its ring this segment belongs to. C++: + /// `sectionalize`, whose sections are what `get_turns` iterates over. + section: usize, } impl

SourceSegment

@@ -100,37 +136,163 @@ 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), (1.0, end)], + splits: alloc::vec![(0.0, start, false), (1.0, end, false)], + section, } } 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, + /// 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: `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], + /// 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], + /// 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. + offset: [f64; 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], + /// 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 { + fn joins(&self, other: &Self) -> bool { + self.start == other.start && self.end == other.end + } +} + +/// Where a turn sits in `get_turns`' collection 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 { + pair_rank: usize, + arrivals: [usize; 2], + offset: f64, +} + +impl TurnOrder { + fn of

(node: &Node

) -> Self { + Self { + pair_rank: node.pair_rank, + arrivals: node.arrival, + offset: node.offset[0], + } + } + + fn compare(&self, other: &Self) -> Ordering { + self.pair_rank + .cmp(&other.pair_rank) + .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 — +/// 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, +} + +impl RingStart { + fn compare(&self, other: &Self) -> Ordering { + 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)) + } + }) + } } /// Execute a polygon Boolean operation through a split-edge arrangement. @@ -146,14 +308,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, false), + polygon_segments(second, operation.walks_second_operand_backwards()), + 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, false), + multi_polygon_segments(second, operation.walks_second_operand_backwards()), + 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); @@ -181,17 +387,29 @@ where &mut first_segments, &mut nodes, &mut candidates, + 0, snap_tolerance, ); append_atomic_edges( &mut second_segments, &mut nodes, &mut candidates, + 1, 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::new(); + let mut boundary: Vec = Vec::new(); for candidate in candidates { let start = nodes[candidate.start].coordinate; let end = nodes[candidate.end].coordinate; @@ -218,21 +436,39 @@ 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 { start: candidate.end, end: candidate.start, + carried_by: candidate.carried_by, + section: candidate.section, } }; - 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]; + held.section[0] = held.section[0].min(edge.section[0]); + held.section[1] = held.section[1].min(edge.section[1]); + } + None => boundary.push(edge), } } let rings = trace_rings(&nodes, &boundary, snap_tolerance)?; - Ok(assemble_multipolygon(rings)) + Ok(assemble_traced(rings)) } fn ring_coordinates(ring: &R) -> Vec @@ -251,40 +487,160 @@ where coordinates } -fn polygon_segments(polygon: &G) -> Vec> +fn multi_polygon_segments(multi_polygon: &G, backwards: bool) -> Vec> +where + G: MultiPolygonTrait, + P: Point + Copy, + 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, &mut sections, backwards); + for ring in polygon.interiors() { + append_ring_segments(ring, &mut segments, &mut sections, backwards); + } + } + segments +} + +fn polygon_segments(polygon: &G, backwards: bool) -> Vec> where G: PolygonTrait, P: Point + Copy, 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, backwards); for ring in polygon.interiors() { - append_ring_segments(ring, &mut segments); + append_ring_segments(ring, &mut segments, &mut sections, backwards); } segments } -fn append_ring_segments(ring: &R, output: &mut Vec>) +/// 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; + +/// 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, + 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; } + 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)); } } @@ -292,27 +648,46 @@ fn append_atomic_edges

( segments: &mut [SourceSegment

], nodes: &mut Vec>, output: &mut Vec, + operand: usize, tolerance: f64, ) where P: Point + Copy, P::Scalar: Into, { - for segment in segments { + for (index, segment) in segments.iter_mut().enumerate() { + let section = segment.section; segment .splits .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); + // 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[operand] == usize::MAX { + nodes[end].arrival[operand] = index; + nodes[end].offset[operand] = pair[1].0; + nodes[end].section[operand] = section; + } if start != end { - output.push(Edge { 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, + }); } } } } -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, @@ -324,23 +699,36 @@ 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, + arrival: [usize::MAX; 2], + section: [usize::MAX; 2], + pair_rank: usize::MAX, + 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()]; - let mut rings = 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; @@ -348,32 +736,212 @@ where let first = edges[seed].start; let mut edge_index = seed; let mut node_indices = alloc::vec![first]; + // `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); + 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 + // 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); + let loop_along = along.split_off(start); + node_indices.push(edge.end); + push_ring(&mut rings, nodes, &loop_nodes, &loop_along, 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(), - )); + } + + // 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. + rings.sort_by(|(left, _), (right, _)| left.compare(right)); + Ok(rings + .into_iter() + .map(|(start, ring)| (ring, start.traversed)) + .collect()) +} + +/// Append a turn point, dropping whatever it now runs straight through. +/// +/// 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. +/// +/// 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 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; + } + } + 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; + } + let last = points.pop().expect("just pushed"); + points.pop(); + points.push(last); + } +} + +/// Append one traced cycle, dropping it when it encloses no area. +/// +/// 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<(RingStart, Ring

)>, + nodes: &[Node

], + node_indices: &[usize], + along: &[Edge], + 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 { + return; + } + + // The cycle is closed, so its last index repeats its first. + 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. + // 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(|&(_, left), &(_, right)| { + 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 + // 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() + .min_by_key(|&(_, index)| index) + .map(|(position, _)| position) + }; + let first_turn = first_turn_by_arrival.or_else(first_node).unwrap_or(0); + + // 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. + // + // 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() { + if first_turn_by_arrival.is_some() { + append_no_collinear(&mut points, first); + } else { + points.push(first); } } - Ok(rings) + // 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 = 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], + }, + Ring::from_vec(points), + )); } fn next_edge

(nodes: &[Node

], edges: &[Edge], used: &[bool], incoming: Edge) -> Option @@ -475,20 +1043,50 @@ mod tests { Node { point: P::new(0.0, 0.0), coordinate: Coordinate { x: 0.0, y: 0.0 }, + is_turn: false, + arrival: [0, 0], + section: [0, 0], + pair_rank: 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, 1], + section: [1, 1], + pair_rank: 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, 2], + section: [2, 2], + pair_rank: 2, + offset: [0.0; 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], + 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], + }, ]; assert!(trace_rings(&nodes, &edges, 1e-10).unwrap().is_empty()); 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/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/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-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`). 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}" + ); + } } diff --git a/crates/geometry-overlay/src/validity.rs b/crates/geometry-overlay/src/validity.rs index 005648b..35d8a71 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 + }); } } } @@ -540,15 +563,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 +586,11 @@ where return Err(ValidityFailure::WrongOrientation); } + // No non-adjacent edge may intersect another. + if has_self_intersection(&pts) { + return Err(ValidityFailure::SelfIntersection); + } + Ok(()) } @@ -634,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); } } @@ -760,11 +801,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, @@ -842,7 +887,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 +900,7 @@ mod tests { ]); assert_eq!( is_valid_ring(&small_bowtie), - Err(ValidityFailure::SelfIntersection) + Err(ValidityFailure::WrongOrientation) ); } @@ -869,7 +917,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 +929,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/overlay_parity.rs b/crates/geometry-overlay/tests/overlay_parity.rs index 9744b6c..126e6e5 100644 --- a/crates/geometry-overlay/tests/overlay_parity.rs +++ b/crates/geometry-overlay/tests/overlay_parity.rs @@ -15,8 +15,11 @@ 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_trait::{MultiPolygon as _, Polygon as _}; +use geometry_overlay::{ + difference, difference_multi, intersection, intersection_multi, sym_difference, union_multi, + union_poly, +}; +use geometry_trait::{MultiPolygon as _, Point as _, Polygon as _, Ring as _}; type P = Point2D; @@ -170,3 +173,587 @@ 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"); +} + +// ---- 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); +} + +// ---- 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 + ); +} + +// ---- 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), + ] + ); +} + +/// A square and a rectangle overlapping along part of one side, running the +/// 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_turn_that_carries_the_outline_straight_on_is_replaced() { + 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), + ] + ); +} + +/// 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:?}" + ); +} + +/// 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] + ); +} + +/// 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)] + ); +} + +/// 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) + ] + ); +} 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) ); } diff --git a/crates/geometry-proj/CHANGELOG.md b/crates/geometry-proj/CHANGELOG.md index 77935a1..4bb1464 100644 --- a/crates/geometry-proj/CHANGELOG.md +++ b/crates/geometry-proj/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-proj-v0.0.8...geometry-proj-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-proj-v0.0.7...geometry-proj-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-proj/Cargo.toml b/crates/geometry-proj/Cargo.toml index 71c9b63..4a33261 100644 --- a/crates/geometry-proj/Cargo.toml +++ b/crates/geometry-proj/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-proj" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -12,10 +12,10 @@ readme = "README.md" [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords" } -geometry-cs = { version = "0.0.8", path = "../geometry-cs" } -geometry-trait = { version = "0.0.8", path = "../geometry-trait" } -geometry-model = { version = "0.0.8", path = "../geometry-model" } +geometry-coords = { version = "0.0.9", path = "../geometry-coords" } +geometry-cs = { version = "0.0.9", path = "../geometry-cs" } +geometry-trait = { version = "0.0.9", path = "../geometry-trait" } +geometry-model = { version = "0.0.9", path = "../geometry-model" } proj4rs = { version = "0.1", features = ["crs-definitions"] } proj4wkt = "0.1" diff --git a/crates/geometry-rtree/CHANGELOG.md b/crates/geometry-rtree/CHANGELOG.md index 0be098b..6b9b9b3 100644 --- a/crates/geometry-rtree/CHANGELOG.md +++ b/crates/geometry-rtree/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-rtree-v0.0.8...geometry-rtree-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs, geometry-trait, geometry-model + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-rtree-v0.0.7...geometry-rtree-v0.0.8) - 2026-07-15 ### Fixed diff --git a/crates/geometry-rtree/Cargo.toml b/crates/geometry-rtree/Cargo.toml index 3d5997c..e8f52bc 100644 --- a/crates/geometry-rtree/Cargo.toml +++ b/crates/geometry-rtree/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-rtree" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -27,10 +27,10 @@ libm = [ serde = ["dep:serde", "geometry-model/serde"] [dependencies] -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false, features = ["libm"] } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false, features = ["libm"] } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } heapless = { version = "0.8", default-features = false } serde = { version = "1", optional = true, default-features = false, features = ["alloc", "derive"] } diff --git a/crates/geometry-strategy/CHANGELOG.md b/crates/geometry-strategy/CHANGELOG.md index 59340ff..d2d1c7b 100644 --- a/crates/geometry-strategy/CHANGELOG.md +++ b/crates/geometry-strategy/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-strategy-v0.0.8...geometry-strategy-v0.0.9) - 2026-09-05 + +### Added + +- *(strategy)* centroid for multi-polygons, accumulated as Boost does + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-strategy-v0.0.7...geometry-strategy-v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry-strategy/Cargo.toml b/crates/geometry-strategy/Cargo.toml index a9c0466..d6b5982 100644 --- a/crates/geometry-strategy/Cargo.toml +++ b/crates/geometry-strategy/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-strategy" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -26,10 +26,10 @@ libm = [ [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } [dev-dependencies] geometry-adapt = { path = "../geometry-adapt" } 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}; diff --git a/crates/geometry-trait/CHANGELOG.md b/crates/geometry-trait/CHANGELOG.md index 604c3be..c79a79b 100644 --- a/crates/geometry-trait/CHANGELOG.md +++ b/crates/geometry-trait/CHANGELOG.md @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/geometry-trait-v0.0.8...geometry-trait-v0.0.9) - 2026-09-05 + +### Other + +- updated the following local packages: geometry-coords, geometry-cs + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/geometry-trait-v0.0.7...geometry-trait-v0.0.8) - 2026-07-15 ### Other diff --git a/crates/geometry-trait/Cargo.toml b/crates/geometry-trait/Cargo.toml index 4977961..d215581 100644 --- a/crates/geometry-trait/Cargo.toml +++ b/crates/geometry-trait/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "geometry-trait" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -17,8 +17,8 @@ libm = ["geometry-coords/libm"] [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } [dev-dependencies] trybuild = "1" diff --git a/crates/geometry/CHANGELOG.md b/crates/geometry/CHANGELOG.md index a4b8296..053f8bd 100644 --- a/crates/geometry/CHANGELOG.md +++ b/crates/geometry/CHANGELOG.md @@ -7,6 +7,17 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +## [0.0.9](https://github.com/pentatonick/boost_geometry/compare/v0.0.8...v0.0.9) - 2026-09-05 + +### Added + +- *(buffer)* a zero-width buffer of a polygon, which is not a no-op + +### Fixed + +- *(overlay)* a hole sharing an edge with the exterior is a self-intersection +- *(overlay)* distinguish the two ways multi-polygon members can be wrong + ## [0.0.8](https://github.com/pentatonick/boost_geometry/compare/v0.0.7...v0.0.8) - 2026-07-15 ### Added diff --git a/crates/geometry/Cargo.toml b/crates/geometry/Cargo.toml index c7a26f0..2d37b05 100644 --- a/crates/geometry/Cargo.toml +++ b/crates/geometry/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "boost_geometry" -version = "0.0.8" +version = "0.0.9" edition.workspace = true rust-version.workspace = true license.workspace = true @@ -41,15 +41,15 @@ serde = ["geometry-rtree/serde"] [dependencies] geometry-tag = { version = "0.0.7", path = "../geometry-tag" } -geometry-coords = { version = "0.0.8", path = "../geometry-coords", default-features = false } -geometry-cs = { version = "0.0.8", path = "../geometry-cs", default-features = false } -geometry-trait = { version = "0.0.8", path = "../geometry-trait", default-features = false } -geometry-model = { version = "0.0.8", path = "../geometry-model", default-features = false } -geometry-strategy = { version = "0.0.8", path = "../geometry-strategy", default-features = false } -geometry-algorithm = { version = "0.0.8", path = "../geometry-algorithm", default-features = false } -geometry-adapt = { version = "0.0.8", path = "../geometry-adapt", default-features = false } -geometry-overlay = { version = "0.0.8", path = "../geometry-overlay", default-features = false } -geometry-rtree = { version = "0.0.8", path = "../geometry-rtree", default-features = false } +geometry-coords = { version = "0.0.9", path = "../geometry-coords", default-features = false } +geometry-cs = { version = "0.0.9", path = "../geometry-cs", default-features = false } +geometry-trait = { version = "0.0.9", path = "../geometry-trait", default-features = false } +geometry-model = { version = "0.0.9", path = "../geometry-model", default-features = false } +geometry-strategy = { version = "0.0.9", path = "../geometry-strategy", default-features = false } +geometry-algorithm = { version = "0.0.9", path = "../geometry-algorithm", default-features = false } +geometry-adapt = { version = "0.0.9", path = "../geometry-adapt", default-features = false } +geometry-overlay = { version = "0.0.9", path = "../geometry-overlay", default-features = false } +geometry-rtree = { version = "0.0.9", path = "../geometry-rtree", default-features = false } geometry-derive = { version = "0.0.7", path = "../geometry-derive" } [dev-dependencies] 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) diff --git a/crates/geometry/tests/validity_completion_parity.rs b/crates/geometry/tests/validity_completion_parity.rs index d5130d8..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), @@ -127,6 +167,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 +190,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) );