geometry_overlay/operation/
boolean.rs

1//! OVL5 — the boolean overlay free functions.
2//!
3//! The public entries route through a planar split-edge arrangement that
4//! performs turn collection, colocation handling, boundary classification,
5//! traversal, and [`assemble`](mod@crate::assemble). Mirrors
6//! `boost/geometry/algorithms/intersection.hpp`, `union.hpp`,
7//! `difference.hpp`, and `sym_difference.hpp`.
8//!
9//! # Where these live
10//!
11//! The overlay plan (`phase_03-…-overlay.md` §OVL5) placed these in
12//! `geometry-algorithm`. That would form a dependency cycle:
13//! `geometry-overlay` already depends on `geometry-algorithm` (for
14//! `within` / `ring_area`), so `geometry-algorithm` cannot depend back
15//! on `geometry-overlay`. The functions therefore live here in
16//! `geometry-overlay` and are re-exported by the `geometry` facade.
17//! This is the same class of spec-stub cycle already corrected
18//! elsewhere in the port.
19//!
20//! Polygon × polygon → `MultiPolygon`, including interior rings, contained
21//! holes/islands, shared edges, and colocated vertices. Coordinates outside
22//! the exact-predicate range surface as [`OverlayError::Unsupported`].
23
24use geometry_coords::CoordinateScalar;
25use geometry_cs::{CartesianFamily, CoordinateSystem};
26use geometry_model::{MultiPolygon, Polygon};
27use geometry_tag::SameAs;
28use geometry_trait::{PointMut, Polygon as PolygonTrait};
29
30use crate::traverse::TraversalError;
31
32use super::areal::{ArealOp, overlay as areal_overlay};
33
34/// Failure of a boolean overlay operation.
35#[derive(Debug, Clone, Copy, PartialEq, Eq)]
36pub enum OverlayError {
37    /// Coordinates exceeded the predicate range or the input could not form
38    /// a supported result boundary. Also propagated from legacy
39    /// [`TraversalError`] callers.
40    Unsupported,
41}
42
43impl From<TraversalError> for OverlayError {
44    fn from(_: TraversalError) -> Self {
45        OverlayError::Unsupported
46    }
47}
48
49/// Intersection of two polygons — the region inside **both**.
50///
51/// Mirrors `boost::geometry::intersection` from
52/// `algorithms/detail/intersection/interface.hpp:342-372`. Returns an empty
53/// `MultiPolygon` when the polygons do not overlap.
54///
55/// # Errors
56///
57/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range.
58///
59/// # Examples
60///
61/// ```
62/// use geometry_cs::Cartesian;
63/// use geometry_model::{polygon, Point2D, Polygon};
64/// use geometry_overlay::operation::intersection;
65/// use geometry_trait::MultiPolygon as _;
66///
67/// type P = Point2D<f64, Cartesian>;
68/// let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
69/// let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
70/// let out = intersection(&a, &b).unwrap();
71/// assert_eq!(out.polygons().count(), 1);
72/// ```
73#[inline]
74#[must_use = "intersection can fail and the resulting geometry should be used"]
75pub fn intersection<G1, G2, P>(g1: &G1, g2: &G2) -> Result<MultiPolygon<Polygon<P>>, OverlayError>
76where
77    G1: PolygonTrait<Point = P>,
78    G2: PolygonTrait<Point = P>,
79    P: PointMut + Default + Copy,
80    P::Scalar: CoordinateScalar + Into<f64>,
81    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
82{
83    areal_overlay(g1, g2, ArealOp::Intersection)
84}
85
86/// Union of two polygons — the region inside **either**.
87///
88/// Mirrors `boost::geometry::union_` from `algorithms/union.hpp:851-881`; the
89/// C++ trailing underscore dodges the keyword, while this compatibility entry
90/// uses the unambiguous name `union_poly`.
91///
92/// # Errors
93///
94/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range.
95///
96/// # Examples
97///
98/// ```
99/// use geometry_cs::Cartesian;
100/// use geometry_model::{polygon, Point2D, Polygon};
101/// use geometry_overlay::operation::union_poly;
102/// use geometry_trait::MultiPolygon as _;
103///
104/// type P = Point2D<f64, Cartesian>;
105/// let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
106/// let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
107/// let out = union_poly(&a, &b).unwrap();
108/// assert_eq!(out.polygons().count(), 1);
109/// ```
110#[inline]
111#[must_use = "union can fail and the resulting geometry should be used"]
112pub fn union_poly<G1, G2, P>(g1: &G1, g2: &G2) -> Result<MultiPolygon<Polygon<P>>, OverlayError>
113where
114    G1: PolygonTrait<Point = P>,
115    G2: PolygonTrait<Point = P>,
116    P: PointMut + Default + Copy,
117    P::Scalar: CoordinateScalar + Into<f64>,
118    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
119{
120    areal_overlay(g1, g2, ArealOp::Union)
121}
122
123/// Union of two polygons — the region inside either input.
124///
125/// This is the Boost-style public spelling of [`union_poly`]. Rust reserves
126/// `union`, so callers write the raw identifier `r#union(a, b)`; the exported
127/// symbol is still named `union`.
128///
129/// Mirrors `boost::geometry::union_` from
130/// `boost/geometry/algorithms/union.hpp:866-880`.
131///
132/// # Errors
133///
134/// Propagates [`OverlayError::Unsupported`] from [`union_poly`].
135#[inline]
136#[must_use = "union can fail and the resulting geometry should be used"]
137pub fn r#union<G1, G2, P>(g1: &G1, g2: &G2) -> Result<MultiPolygon<Polygon<P>>, OverlayError>
138where
139    G1: PolygonTrait<Point = P>,
140    G2: PolygonTrait<Point = P>,
141    P: PointMut + Default + Copy,
142    P::Scalar: CoordinateScalar + Into<f64>,
143    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
144{
145    union_poly(g1, g2)
146}
147
148/// Difference of two polygons — the region inside the first but outside
149/// the second (`A − B`).
150///
151/// Mirrors `boost::geometry::difference` from
152/// `algorithms/difference.hpp:686-714`.
153///
154/// # Errors
155///
156/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range.
157///
158/// # Examples
159///
160/// ```
161/// use geometry_cs::Cartesian;
162/// use geometry_model::{polygon, Point2D, Polygon};
163/// use geometry_overlay::operation::difference;
164/// use geometry_trait::MultiPolygon as _;
165///
166/// type P = Point2D<f64, Cartesian>;
167/// let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
168/// let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
169/// let out = difference(&a, &b).unwrap();
170/// assert_eq!(out.polygons().count(), 1);
171/// ```
172#[inline]
173#[must_use = "difference can fail and the resulting geometry should be used"]
174pub fn difference<G1, G2, P>(g1: &G1, g2: &G2) -> Result<MultiPolygon<Polygon<P>>, OverlayError>
175where
176    G1: PolygonTrait<Point = P>,
177    G2: PolygonTrait<Point = P>,
178    P: PointMut + Default + Copy,
179    P::Scalar: CoordinateScalar + Into<f64>,
180    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
181{
182    areal_overlay(g1, g2, ArealOp::Difference)
183}
184
185/// Symmetric difference of two polygons — the region inside exactly one
186/// of them (`(A − B) ∪ (B − A)`).
187///
188/// Mirrors `boost::geometry::sym_difference` from
189/// `algorithms/sym_difference.hpp:795-824`.
190///
191/// # Errors
192///
193/// [`OverlayError::Unsupported`] when coordinates exceed the predicate range.
194///
195/// # Examples
196///
197/// ```
198/// use geometry_cs::Cartesian;
199/// use geometry_model::{polygon, Point2D, Polygon};
200/// use geometry_overlay::operation::sym_difference;
201/// use geometry_trait::MultiPolygon as _;
202///
203/// type P = Point2D<f64, Cartesian>;
204/// let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
205/// let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
206/// let out = sym_difference(&a, &b).unwrap();
207/// assert!(out.polygons().count() >= 1);
208/// ```
209#[inline]
210#[must_use = "symmetric difference can fail and the resulting geometry should be used"]
211pub fn sym_difference<G1, G2, P>(g1: &G1, g2: &G2) -> Result<MultiPolygon<Polygon<P>>, OverlayError>
212where
213    G1: PolygonTrait<Point = P>,
214    G2: PolygonTrait<Point = P>,
215    P: PointMut + Default + Copy,
216    P::Scalar: CoordinateScalar + Into<f64>,
217    <P::Cs as CoordinateSystem>::Family: SameAs<CartesianFamily>,
218{
219    areal_overlay(g1, g2, ArealOp::SymDifference)
220}
221
222#[cfg(test)]
223mod tests {
224    use super::{OverlayError, intersection, union_poly};
225    use geometry_algorithm::area;
226    use geometry_cs::Cartesian;
227    use geometry_model::{Point2D, Polygon, polygon};
228    use geometry_trait::{MultiPolygon as _, Polygon as _};
229
230    type P = Point2D<f64, Cartesian>;
231
232    fn close(a: f64, b: f64) -> bool {
233        (a - b).abs() <= 1e-5 * a.abs().max(b.abs()).max(1.0)
234    }
235
236    fn total_area(mp: &geometry_model::MultiPolygon<Polygon<P>>) -> f64 {
237        mp.polygons().map(|pg| area(pg).abs()).sum()
238    }
239
240    #[test]
241    fn intersection_of_offset_squares() {
242        let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
243        let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
244        let out = intersection(&a, &b).unwrap();
245        assert_eq!(out.polygons().count(), 1);
246        assert!(close(total_area(&out), 1.0), "area {}", total_area(&out));
247    }
248
249    #[test]
250    fn intersection_disjoint_is_empty() {
251        let a: Polygon<P> = polygon![[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]];
252        let b: Polygon<P> = polygon![[(5.0, 5.0), (6.0, 5.0), (6.0, 6.0), (5.0, 6.0), (5.0, 5.0)]];
253        let out = intersection(&a, &b).unwrap();
254        assert_eq!(out.polygons().count(), 0);
255    }
256
257    #[test]
258    fn intersection_contained_is_inner() {
259        let big: Polygon<P> = polygon![[
260            (0.0, 0.0),
261            (10.0, 0.0),
262            (10.0, 10.0),
263            (0.0, 10.0),
264            (0.0, 0.0)
265        ]];
266        let small: Polygon<P> =
267            polygon![[(2.0, 2.0), (4.0, 2.0), (4.0, 4.0), (2.0, 4.0), (2.0, 2.0)]];
268        let out = intersection(&big, &small).unwrap();
269        assert_eq!(out.polygons().count(), 1);
270        assert!(close(total_area(&out), 4.0), "area {}", total_area(&out));
271    }
272
273    #[test]
274    fn union_of_offset_squares_area() {
275        // |A| + |B| - |A∩B| = 4 + 4 - 1 = 7.
276        let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
277        let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
278        let out = union_poly(&a, &b).unwrap();
279        assert_eq!(out.polygons().count(), 1);
280        assert!(close(total_area(&out), 7.0), "area {}", total_area(&out));
281    }
282
283    #[test]
284    fn union_disjoint_is_two_polygons() {
285        let a: Polygon<P> = polygon![[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]];
286        let b: Polygon<P> = polygon![[(5.0, 5.0), (6.0, 5.0), (6.0, 6.0), (5.0, 6.0), (5.0, 5.0)]];
287        let out = union_poly(&a, &b).unwrap();
288        assert_eq!(out.polygons().count(), 2);
289        assert!(close(total_area(&out), 2.0), "area {}", total_area(&out));
290    }
291
292    #[test]
293    fn difference_of_offset_squares_area() {
294        // |A − B| = |A| − |A∩B| = 4 − 1 = 3.
295        let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
296        let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
297        let out = super::difference(&a, &b).unwrap();
298        assert_eq!(out.polygons().count(), 1);
299        assert!(close(total_area(&out), 3.0), "area {}", total_area(&out));
300    }
301
302    #[test]
303    fn difference_disjoint_is_first_whole() {
304        let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
305        let b: Polygon<P> = polygon![[(5.0, 5.0), (6.0, 5.0), (6.0, 6.0), (5.0, 6.0), (5.0, 5.0)]];
306        let out = super::difference(&a, &b).unwrap();
307        assert_eq!(out.polygons().count(), 1);
308        assert!(close(total_area(&out), 4.0), "area {}", total_area(&out));
309    }
310
311    #[test]
312    fn difference_a_inside_b_is_empty() {
313        let big: Polygon<P> = polygon![[
314            (0.0, 0.0),
315            (10.0, 0.0),
316            (10.0, 10.0),
317            (0.0, 10.0),
318            (0.0, 0.0)
319        ]];
320        let small: Polygon<P> =
321            polygon![[(2.0, 2.0), (4.0, 2.0), (4.0, 4.0), (2.0, 4.0), (2.0, 2.0)]];
322        let out = super::difference(&small, &big).unwrap();
323        assert_eq!(out.polygons().count(), 0);
324    }
325
326    #[test]
327    fn difference_with_contained_subtrahend_emits_a_hole() {
328        let big: Polygon<P> = polygon![[
329            (0.0, 0.0),
330            (10.0, 0.0),
331            (10.0, 10.0),
332            (0.0, 10.0),
333            (0.0, 0.0)
334        ]];
335        let small: Polygon<P> =
336            polygon![[(3.0, 3.0), (5.0, 3.0), (5.0, 5.0), (3.0, 5.0), (3.0, 3.0)]];
337        let difference = super::difference(&big, &small).unwrap();
338        assert_eq!(difference.polygons().count(), 1);
339        assert_eq!(difference.polygons().next().unwrap().interiors().count(), 1);
340        assert!(close(total_area(&difference), 96.0));
341        assert!(close(
342            total_area(&super::sym_difference(&big, &small).unwrap()),
343            96.0
344        ));
345    }
346
347    #[test]
348    fn input_with_holes_participates_in_all_operations() {
349        let donut: Polygon<P> = polygon![
350            [
351                (0.0, 0.0),
352                (10.0, 0.0),
353                (10.0, 10.0),
354                (0.0, 10.0),
355                (0.0, 0.0)
356            ],
357            [(3.0, 3.0), (7.0, 3.0), (7.0, 7.0), (3.0, 7.0), (3.0, 3.0)]
358        ];
359        let sq: Polygon<P> = polygon![[(2.0, 2.0), (8.0, 2.0), (8.0, 8.0), (2.0, 8.0), (2.0, 2.0)]];
360        assert!(close(total_area(&intersection(&donut, &sq).unwrap()), 20.0));
361        assert!(close(total_area(&union_poly(&donut, &sq).unwrap()), 100.0));
362        assert!(close(
363            total_area(&super::difference(&donut, &sq).unwrap()),
364            64.0
365        ));
366        assert!(close(
367            total_area(&super::sym_difference(&donut, &sq).unwrap()),
368            80.0
369        ));
370    }
371
372    #[test]
373    fn out_of_range_coordinates_are_refused_not_silently_wrong() {
374        // Regression: two huge overlapping squares (~1e14, past the ±2^26
375        // safe range) made the turn kernel silently drop every crossing as
376        // OutOfRange; the emptied turn graph was misread as "B inside A",
377        // over-reporting the intersection area ~4× as `Ok`. All ops must
378        // refuse rather than return a silently wrong result.
379        let a: Polygon<P> = polygon![[
380            (0.0, 0.0),
381            (2e14, 0.0),
382            (2e14, 2e14),
383            (0.0, 2e14),
384            (0.0, 0.0)
385        ]];
386        let b: Polygon<P> = polygon![[
387            (1e14, 1e14),
388            (3e14, 1e14),
389            (3e14, 3e14),
390            (1e14, 3e14),
391            (1e14, 1e14)
392        ]];
393        assert_eq!(intersection(&a, &b), Err(OverlayError::Unsupported));
394        assert_eq!(union_poly(&a, &b), Err(OverlayError::Unsupported));
395        assert_eq!(super::difference(&a, &b), Err(OverlayError::Unsupported));
396        assert_eq!(
397            super::sym_difference(&a, &b),
398            Err(OverlayError::Unsupported)
399        );
400    }
401
402    #[test]
403    fn sym_difference_of_offset_squares_area() {
404        // |A △ B| = |A| + |B| − 2|A∩B| = 4 + 4 − 2·1 = 6.
405        let a: Polygon<P> = polygon![[(0.0, 0.0), (2.0, 0.0), (2.0, 2.0), (0.0, 2.0), (0.0, 0.0)]];
406        let b: Polygon<P> = polygon![[(1.0, 1.0), (3.0, 1.0), (3.0, 3.0), (1.0, 3.0), (1.0, 1.0)]];
407        let out = super::sym_difference(&a, &b).unwrap();
408        assert!(close(total_area(&out), 6.0), "area {}", total_area(&out));
409    }
410
411    #[test]
412    fn union_contained_is_outer() {
413        let big: Polygon<P> = polygon![[
414            (0.0, 0.0),
415            (10.0, 0.0),
416            (10.0, 10.0),
417            (0.0, 10.0),
418            (0.0, 0.0)
419        ]];
420        let small: Polygon<P> =
421            polygon![[(2.0, 2.0), (4.0, 2.0), (4.0, 4.0), (2.0, 4.0), (2.0, 2.0)]];
422        let out = union_poly(&big, &small).unwrap();
423        assert_eq!(out.polygons().count(), 1);
424        assert!(close(total_area(&out), 100.0), "area {}", total_area(&out));
425    }
426}
geometry_overlay/operation/boolean.rs

geometry_overlay/operation/
boolean.rs
