Merge pull request #1246 from urschrei/cascaded_union

Cascaded / Unary union

Author: michaelkirk

geo-test-fixtures/fixtures/nl_plots_epsg_28992.wkt

@@ -126,13 +126,21 @@ where
 }
 
 // From https://afnemers.ruimtelijkeplannen.nl/afnemers/services?request=GetFeature&service=WFS&srsName=EPSG:4326&typeName=Enkelbestemming&version=2.0.0&bbox=165618,480983,166149,481542";
-pub fn nl_plots<T>() -> MultiPolygon<T>
+pub fn nl_plots_wgs84<T>() -> MultiPolygon<T>
 where
     T: WktFloat + Default + FromStr,
 {
     multi_polygon("nl_plots.wkt")
 }
 
+pub fn nl_plots_epsg_28992<T>() -> MultiPolygon<T>
+where
+    T: WktFloat + Default + FromStr,
+{
+    // https://epsg.io/28992
+    multi_polygon("nl_plots_epsg_28992.wkt")
+}
+
 fn line_string<T>(name: &str) -> LineString<T>
 where
     T: WktFloat + Default + FromStr,

geo-test-fixtures/src/lib.rs

@@ -2,6 +2,8 @@
 
 ## Unreleased
 
+- Add Unary Union algorithm for fast union ops on adjacent / overlapping geometries
+  - <https://github.com/georust/geo/pull/1246>
 - Loosen bounds on `RemoveRepeatedPoints` trait (`num_traits::FromPrimitive` isn't required)
   - <https://github.com/georust/geo/pull/1278>
 

geo/CHANGES.md

@@ -11,7 +11,7 @@ use criterion::Criterion;
 
 fn criterion_benchmark(c: &mut Criterion) {
     c.bench_function("Point position to rect", |bencher| {
-        let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots()
+        let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots_wgs84()
             .iter()
             .map(|plot| plot.centroid().unwrap())
             .collect();

geo/benches/coordinate_position.rs

@@ -3,7 +3,7 @@ use geo::intersects::Intersects;
 use geo::MultiPolygon;
 
 fn multi_polygon_intersection(c: &mut Criterion) {
-    let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots();
+    let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots_wgs84();
     let zone_polygons: MultiPolygon = geo_test_fixtures::nl_zones();
 
     c.bench_function("MultiPolygon intersects", |bencher| {
@@ -30,7 +30,7 @@ fn multi_polygon_intersection(c: &mut Criterion) {
 fn rect_intersection(c: &mut Criterion) {
     use geo::algorithm::BoundingRect;
     use geo::geometry::Rect;
-    let plot_bbox: Vec<Rect> = geo_test_fixtures::nl_plots()
+    let plot_bbox: Vec<Rect> = geo_test_fixtures::nl_plots_wgs84()
         .iter()
         .map(|plot| plot.bounding_rect().unwrap())
         .collect();
@@ -63,7 +63,7 @@ fn rect_intersection(c: &mut Criterion) {
 fn point_rect_intersection(c: &mut Criterion) {
     use geo::algorithm::{BoundingRect, Centroid};
     use geo::geometry::{Point, Rect};
-    let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots()
+    let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots_wgs84()
         .iter()
         .map(|plot| plot.centroid().unwrap())
         .collect();
@@ -96,7 +96,7 @@ fn point_rect_intersection(c: &mut Criterion) {
 fn point_triangle_intersection(c: &mut Criterion) {
     use geo::{Centroid, TriangulateEarcut};
     use geo_types::{Point, Triangle};
-    let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots()
+    let plot_centroids: Vec<Point> = geo_test_fixtures::nl_plots_wgs84()
         .iter()
         .map(|plot| plot.centroid().unwrap())
         .collect();

geo/benches/intersection.rs

@@ -5,7 +5,7 @@ use geo_types::MultiPolygon;
 
 fn criterion_benchmark(c: &mut Criterion) {
     c.bench_function("relate prepared polygons", |bencher| {
-        let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots();
+        let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots_wgs84();
         let zone_polygons = geo_test_fixtures::nl_zones();
 
         bencher.iter(|| {
@@ -38,7 +38,7 @@ fn criterion_benchmark(c: &mut Criterion) {
     });
 
     c.bench_function("relate unprepared polygons", |bencher| {
-        let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots();
+        let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots_wgs84();
         let zone_polygons = geo_test_fixtures::nl_zones();
 
         bencher.iter(|| {

geo/benches/prepared_geometry.rs

@@ -2,23 +2,25 @@ mod i_overlay_integration;
 #[cfg(test)]
 mod tests;
 
-use crate::bool_ops::i_overlay_integration::convert::{
-    multi_polygon_from_shapes, ring_to_shape_path,
-};
-use crate::bool_ops::i_overlay_integration::BoolOpsCoord;
+use i_overlay_integration::convert::{multi_polygon_from_shapes, ring_to_shape_path};
+use i_overlay_integration::BoolOpsCoord;
+pub use i_overlay_integration::BoolOpsNum;
+
+use crate::geometry::{LineString, MultiLineString, MultiPolygon, Polygon};
+use crate::winding_order::{Winding, WindingOrder};
+
 use i_overlay::core::fill_rule::FillRule;
+use i_overlay::core::overlay_rule::OverlayRule;
 use i_overlay::float::clip::FloatClip;
+use i_overlay::float::overlay::FloatOverlay;
 use i_overlay::float::single::SingleFloatOverlay;
 use i_overlay::string::clip::ClipRule;
-pub use i_overlay_integration::BoolOpsNum;
-
-use crate::geometry::{LineString, MultiLineString, MultiPolygon, Polygon};
 
 /// Boolean Operations on geometry.
 ///
 /// Boolean operations are set operations on geometries considered as a subset
-/// of the 2-D plane. The operations supported are: intersection, union, xor or
-/// symmetric difference, and set-difference on pairs of 2-D geometries and
+/// of the 2-D plane. The operations supported are: intersection, union,
+/// symmetric difference (xor), and set-difference on pairs of 2-D geometries and
 /// clipping a 1-D geometry with self.
 ///
 /// These operations are implemented on [`Polygon`] and the [`MultiPolygon`]
@@ -34,6 +36,11 @@ use crate::geometry::{LineString, MultiLineString, MultiPolygon, Polygon};
 /// In particular, taking `union` with an empty geom should remove degeneracies
 /// and fix invalid polygons as long the interior-exterior requirement above is
 /// satisfied.
+///
+/// # Performance
+///
+/// For union operations on a large number of [`Polygon`]s or [`MultiPolygons`],
+/// using [`unary_union`] will yield far better performance.
 pub trait BooleanOps {
     type Scalar: BoolOpsNum;
 
@@ -57,18 +64,26 @@ pub trait BooleanOps {
         multi_polygon_from_shapes(shapes)
     }
 
+    /// Returns the overlapping regions shared by both `self` and `other`.
     fn intersection(
         &self,
         other: &impl BooleanOps<Scalar = Self::Scalar>,
     ) -> MultiPolygon<Self::Scalar> {
         self.boolean_op(other, OpType::Intersection)
     }
+
+    /// Combines the regions of both `self` and `other` into a single geometry, removing
+    /// overlaps and merging boundaries.
     fn union(&self, other: &impl BooleanOps<Scalar = Self::Scalar>) -> MultiPolygon<Self::Scalar> {
         self.boolean_op(other, OpType::Union)
     }
+
+    /// The regions that are in either `self` or `other`, but not in both.
     fn xor(&self, other: &impl BooleanOps<Scalar = Self::Scalar>) -> MultiPolygon<Self::Scalar> {
         self.boolean_op(other, OpType::Xor)
     }
+
+    /// The regions of `self` which are not in `other`.
     fn difference(
         &self,
         other: &impl BooleanOps<Scalar = Self::Scalar>,
@@ -109,6 +124,75 @@ pub enum OpType {
     Xor,
 }
 
+/// Efficient [union](BooleanOps::union) of many adjacent / overlapping geometries
+///
+/// This is typically much faster than `union`ing a bunch of geometries together one at a time.
+///
+/// Note: Geometries can be wound in either direction, but the winding order must be consistent,
+/// and the polygon's interiors must be wound opposite to its exterior.
+///
+/// See [Orient] for more information.
+///
+/// [Orient]: crate::algorithm::orient::Orient
+///
+/// # Arguments
+///
+/// `boppables`: A collection of `Polygon` or `MultiPolygons` to union together.
+///
+/// returns the union of all the inputs.
+///
+/// # Examples
+///
+/// ```
+/// use geo::algorithm::unary_union;
+/// use geo::wkt;
+///
+/// let right_piece = wkt!(POLYGON((4. 0.,4. 4.,8. 4.,8. 0.,4. 0.)));
+/// let left_piece = wkt!(POLYGON((0. 0.,0. 4.,4. 4.,4. 0.,0. 0.)));
+///
+/// // touches neither right nor left piece
+/// let separate_piece = wkt!(POLYGON((14. 10.,14. 14.,18. 14.,18. 10.,14. 10.)));
+///
+/// let polygons = vec![left_piece, separate_piece, right_piece];
+/// let actual_output = unary_union(&polygons);
+///
+/// let expected_output = wkt!(MULTIPOLYGON(
+///     // left and right piece have been combined
+///     ((0. 0., 0. 4., 8. 4., 8. 0.,  0. 0.)),
+///     // separate piece remains separate
+///     ((14. 10., 14. 14., 18. 14.,18. 10., 14. 10.))
+/// ));
+/// assert_eq!(actual_output, expected_output);
+/// ```
+pub fn unary_union<'a, B: BooleanOps + 'a>(
+    boppables: impl IntoIterator<Item = &'a B>,
+) -> MultiPolygon<B::Scalar> {
+    let mut winding_order: Option<WindingOrder> = None;
+    let subject = boppables
+        .into_iter()
+        .flat_map(|boppable| {
+            let rings = boppable.rings();
+            rings
+                .map(|ring| {
+                    if winding_order.is_none() {
+                        winding_order = ring.winding_order();
+                    }
+                    ring_to_shape_path(ring)
+                })
+                .collect::<Vec<_>>()
+        })
+        .collect::<Vec<_>>();
+
+    let fill_rule = if winding_order == Some(WindingOrder::Clockwise) {
+        FillRule::Positive
+    } else {
+        FillRule::Negative
+    };
+
+    let shapes = FloatOverlay::with_subj(&subject).overlay(OverlayRule::Subject, fill_rule);
+    multi_polygon_from_shapes(shapes)
+}
+
 impl<T: BoolOpsNum> BooleanOps for Polygon<T> {
     type Scalar = T;
 

geo/src/algorithm/bool_ops/mod.rs

@@ -1,7 +1,94 @@
-use super::BooleanOps;
-use crate::{wkt, Convert, MultiPolygon, Relate};
+use super::{unary_union, BooleanOps};
+use crate::{wkt, Convert, MultiPolygon, Polygon, Relate};
+use std::time::Instant;
 use wkt::ToWkt;
 
+#[test]
+fn test_unary_union() {
+    let poly1: Polygon = wkt!(POLYGON((204.0 287.0,203.69670020700084 288.2213844497616,200.38308697914755 288.338793163584,204.0 287.0)));
+    let poly2: Polygon = wkt!(POLYGON((210.0 290.0,204.07584923592933 288.2701221108328,212.24082541367974 285.47846008552216,210.0 290.0)));
+    let poly3: Polygon = wkt!(POLYGON((211.0 292.0,202.07584923592933 288.2701221108328,212.24082541367974 285.47846008552216,210.0 290.0)));
+
+    let polys = vec![poly1.clone(), poly2.clone(), poly3.clone()];
+    let poly_union = unary_union(&polys);
+    assert_eq!(poly_union.0.len(), 1);
+
+    let multi_poly_12 = MultiPolygon::new(vec![poly1.clone(), poly2.clone()]);
+    let multi_poly_3 = MultiPolygon::new(vec![poly3]);
+    let multi_polys = vec![multi_poly_12.clone(), multi_poly_3.clone()];
+    let multi_poly_union = unary_union(&multi_polys);
+    assert_eq!(multi_poly_union.0.len(), 1);
+}
+
+#[test]
+fn test_unary_union_errors() {
+    let input: MultiPolygon = geo_test_fixtures::nl_plots_epsg_28992();
+
+    assert_eq!(input.0.len(), 316);
+
+    let input_area = input.signed_area();
+    assert_relative_eq!(input_area, 763889.4732974821);
+
+    let naive_union = {
+        let start = Instant::now();
+        let mut output = MultiPolygon::new(Vec::new());
+        for poly in input.iter() {
+            output = output.union(poly);
+        }
+        let union = output;
+        let duration = start.elapsed();
+        println!("Time elapsed (naive): {:.2?}", duration);
+        union
+    };
+
+    let simplified_union = {
+        let start = Instant::now();
+        let union = unary_union(input.iter());
+        let duration = start.elapsed();
+        println!("Time elapsed (simplification): {:.2?}", duration);
+        union
+    };
+
+    use crate::algorithm::Area;
+    let naive_area = naive_union.unsigned_area();
+    let simplified_area = simplified_union.unsigned_area();
+    assert_relative_eq!(naive_area, simplified_area, max_relative = 1e-5);
+
+    // Serial vs. parallel are expected to have slightly different results.
+    //
+    // Each boolean operation scales the floating point to a discrete
+    // integer grid, which introduces some error, and this error factor depends on the magnitude
+    // of the input.
+    //
+    // Because the serial vs. parallel approaches group inputs differently, error is accumulated
+    // differently - hence the slightly different outputs.
+    //
+    // xor'ing the two shapes represents the magnitude of the difference between the two outputs.
+    //
+    // We want to verify that this error is small - it should be near 0, but the
+    // magnitude of the error is relative to the magnitude of the input geometries, so we offset
+    // both the error and 0 by `input_area` to make a scale relative comparison.
+    let naive_vs_simplified_discrepancy = simplified_union.xor(&naive_union);
+    assert_relative_eq!(
+        input_area + naive_vs_simplified_discrepancy.unsigned_area(),
+        0.0 + input_area,
+        max_relative = 1e-5
+    );
+
+    assert_eq!(simplified_union.0.len(), 1);
+    assert_relative_eq!(simplified_area, input_area, max_relative = 1e-5);
+}
+
+#[test]
+fn test_unary_union_winding() {
+    let input: MultiPolygon = geo_test_fixtures::nl_plots_epsg_28992();
+
+    use crate::orient::{Direction, Orient};
+    let default_winding_union = unary_union(input.orient(Direction::Default).iter());
+    let reversed_winding_union = unary_union(input.orient(Direction::Reversed).iter());
+    assert_eq!(default_winding_union, reversed_winding_union);
+}
+
 #[test]
 fn jts_test_overlay_la_1() {
     // From TestOverlayLA.xml test case with description "mLmA - A and B complex, overlapping and touching #1"

geo/src/algorithm/bool_ops/tests.rs

@@ -6,9 +6,9 @@ pub use kernels::{Kernel, Orientation};
 pub mod area;
 pub use area::Area;
 
-/// Boolean Ops such as union, xor, difference.
+/// Boolean Operations such as the union, xor, or difference of two geometries.
 pub mod bool_ops;
-pub use bool_ops::{BooleanOps, OpType};
+pub use bool_ops::{unary_union, BooleanOps, OpType};
 
 /// Calculate the bounding rectangle of a `Geometry`.
 pub mod bounding_rect;