Merge #1018

1018: Monotone Polygon Subdivision r=rmanoka a=rmanoka

- [x] I agree to follow the project's [code of conduct](https://github.com/georust/geo/blob/main/CODE_OF_CONDUCT.md).
- [x] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---


This PR provides a pre-processed representation for a polygon (port from my `geo-crossings` unpublished crate), that allows fast point-in-polygon queries.  This implementation is based on these awesome [lecture notes] by David Mount.  The results are promising:  the pre-processed queries are faster across the board, about 1.5-2x in the worst-case, to upto 100x in the best cases.

In this PR, we provide:

- A `MonoPoly<T>` structure that represents a monotone polygon, supporting fast point-in-polygon queries.

- A builder algorithm: `monotone_subdivision(Polygon<T>)` that sub-divides a polygon into a collection of `MonoPoly<T>`s.

- Benchmarks that show that the pre-processed point-in-polygon queries are fast.

# Benchmarks

We use three types of polygons:
  
  - A steppy polygon that zig-zags along the X-axis as it moves up along Y-axis, then again zig-zags down the Y-axis.  This is the worst case for our algorithm as it gives many monotone polygons along the X-axis.
  
  - Same as above but along Y-axis.  This is a good case as it is representible by a single monotone polygon.
  
  - Circular polygons.  These also give many monotone-polygons on sub-division, and are a good test case.

With each polygon class, we generate polygons with different number of vertices, and measure query times for checking if a random point in the bounding box lies in the polygon.

Try out with `cargo criterion --bench monotone_subdiv`  .

# Remarks / Caveats

- This uses a simpler version of sweep used for the boolean ops.  I've re-implemented the sweep to try to move away from the complicated float issues we've been facing.  For instance, this whole algorithm works without requiring `T: Float`.

- The algorithm currently only supports a single polygon, and no two segments in it can intersect anywhere except at their endpoints.  This is almost the SFS spec, but I think SFS also allows intersections at interior that are "tangential".  I believe this can be handled without involving float ops as the intersection is "exact", but needs more work.

- This data-structure is currently not a stand-in replacement for a `Polygon` for `CoordinatePosition` algorithm as we don't currently remember the new edges added during subdivision.  These edges actually lie inside the polygon, and points on these edges would be on the border of two of the sub-divided `MonoPoly`.  This could be handled by some clever `mod 2` counting but I wasn't sure of the logic.  As it stands, we can get fast equivalent `Intersects` implementation with this impl via a wrapper around a `Vec<MonoPoly<T>>`.

[lecture notes]: //www.cs.umd.edu/class/spring2020/cmsc754/Lects/lect05-triangulate.pdf

Co-authored-by: Rajsekar Manokaran <rajsekar@gmail.com>
Co-authored-by: rmanoka <rajsekar@gmail.com>

geo/CHANGES.md

@@ -18,6 +18,11 @@
 - Add `Vector2DOps` trait to algorithims module and implemented it for `Coord<T::CoordFloat>`
   - <https://github.com/georust/geo/pull/1025>
 
+- Add a fast point-in-polygon query datastructure that pre-processes a `Polygon` as a set of monotone polygons. Ref. `crate::algorithm::MonotonicPolygons`.
+  - <https://github.com/georust/geo/pull/1018>
+
+
+
 ## 0.25.0
 
 - Added `CrossTrackDistance` trait to calculate the distance from a point

geo/Cargo.toml

@@ -105,3 +105,7 @@ harness = false
 [[bench]]
 name = "winding_order"
 harness = false
+
+[[bench]]
+name = "monotone_subdiv"
+harness = false

geo/benches/monotone_subdiv.rs

@@ -0,0 +1,140 @@
+#[macro_use]
+extern crate criterion;
+
+extern crate geo;
+
+use std::fmt::Display;
+use std::panic::catch_unwind;
+
+use criterion::measurement::Measurement;
+use geo::coordinate_position::CoordPos;
+use geo::monotone::monotone_subdivision;
+use geo::{CoordinatePosition, MapCoords, Polygon};
+
+use criterion::{BatchSize, BenchmarkGroup, BenchmarkId, Criterion};
+use geo_types::{Coord, Rect};
+use wkt::ToWkt;
+
+#[path = "utils/random.rs"]
+mod random;
+use rand::thread_rng;
+use random::*;
+
+fn criterion_benchmark_pt_in_poly(c: &mut Criterion) {
+    let pt_samples = Samples::from_fn(512, || {
+        uniform_point(&mut thread_rng(), Rect::new((-1., -1.), (1., 1.)))
+    });
+
+    for size in [16, 64, 512, 1024, 2048] {
+        let mut grp = c.benchmark_group("rand pt-in-poly steppy-polygon (worst case)".to_string());
+        let poly = steppy_polygon(&mut thread_rng(), size);
+        bench_pt_in_poly(&mut grp, poly, size, &pt_samples)
+    }
+    for size in [16, 64, 512, 1024, 2048] {
+        let mut grp = c.benchmark_group("rand pt-in-poly steppy-polygon (best case)".to_string());
+        let poly = steppy_polygon(&mut thread_rng(), size).map_coords(|c| (c.y, c.x).into());
+        bench_pt_in_poly(&mut grp, poly, size, &pt_samples)
+    }
+    for size in [16, 64, 512, 1024, 2048] {
+        let mut grp = c.benchmark_group("rand pt-in-poly circular-polygon".to_string());
+        let poly = circular_polygon(&mut thread_rng(), size);
+        bench_pt_in_poly(&mut grp, poly, size, &pt_samples)
+    }
+}
+
+fn bench_pt_in_poly<T, I>(
+    g: &mut BenchmarkGroup<T>,
+    polygon: Polygon<f64>,
+    param: I,
+    samples: &Samples<Coord<f64>>,
+) where
+    T: Measurement,
+    I: Display + Copy,
+{
+    let mon = match catch_unwind(|| monotone_subdivision([polygon.clone()])) {
+        Ok(m) => m,
+        Err(_) => {
+            panic!(
+                "Monotone subdivision failed for polygon: {}",
+                polygon.to_wkt()
+            );
+        }
+    };
+
+    g.bench_with_input(
+        BenchmarkId::new("Simple point-in-poly", param),
+        &(),
+        |b, _| {
+            b.iter_batched(
+                samples.sampler(),
+                |pt| {
+                    polygon.coordinate_position(pt);
+                },
+                BatchSize::SmallInput,
+            );
+        },
+    );
+    g.bench_with_input(
+        BenchmarkId::new("Pre-processed point-in-poly", param),
+        &(),
+        |b, _| {
+            b.iter_batched(
+                samples.sampler(),
+                |pt| {
+                    mon.iter()
+                        .filter(|mp| mp.coordinate_position(pt) == CoordPos::Inside)
+                        .count();
+                },
+                BatchSize::SmallInput,
+            );
+        },
+    );
+}
+
+fn criterion_benchmark_monotone_subdiv(c: &mut Criterion) {
+    for size in [16, 64, 2048, 32000] {
+        let mut grp = c.benchmark_group("monotone_subdiv steppy-polygon (worst case)".to_string());
+        let poly_fn = |size| steppy_polygon(&mut thread_rng(), size);
+        bench_monotone_subdiv(&mut grp, poly_fn, size)
+    }
+    for size in [16, 64, 2048, 32000] {
+        let mut grp = c.benchmark_group("monotone_subdiv steppy-polygon (best case)".to_string());
+        let poly_fn =
+            |size| steppy_polygon(&mut thread_rng(), size).map_coords(|c| (c.y, c.x).into());
+        bench_monotone_subdiv(&mut grp, poly_fn, size)
+    }
+    for size in [16, 64, 2048, 32000] {
+        let mut grp = c.benchmark_group("monotone_subdiv circular-polygon".to_string());
+        let poly_fn = |size| circular_polygon(&mut thread_rng(), size);
+        bench_monotone_subdiv(&mut grp, poly_fn, size)
+    }
+}
+
+fn bench_monotone_subdiv<T, F>(g: &mut BenchmarkGroup<T>, mut f: F, param: usize)
+where
+    T: Measurement,
+    F: FnMut(usize) -> Polygon<f64>,
+{
+    let samples = Samples::from_fn(16, || f(param));
+    g.bench_with_input(
+        BenchmarkId::new("Montone subdivision", param),
+        &(),
+        |b, _| {
+            b.iter_batched(
+                samples.sampler(),
+                |pt| {
+                    let mon = monotone_subdivision([pt.clone()]);
+                    mon.len();
+                },
+                BatchSize::SmallInput,
+            );
+        },
+    );
+}
+
+criterion_group!(
+    benches,
+    criterion_benchmark_pt_in_poly,
+    criterion_benchmark_monotone_subdiv
+);
+criterion_main!(benches);

geo/benches/utils/random.rs

@@ -1,7 +1,7 @@
 #![allow(unused)]
 use std::f64::consts::PI;
 
-use geo::algorithm::{ConcaveHull, ConvexHull, MapCoords, Rotate};
+use geo::algorithm::{BoundingRect, ConcaveHull, ConvexHull, MapCoords, Rotate};
 use geo::geometry::*;
 
 use rand::{thread_rng, Rng};
@@ -76,32 +76,50 @@ pub fn circular_polygon<R: Rng>(mut rng: R, steps: usize) -> Polygon<f64> {
 pub fn steppy_polygon<R: Rng>(mut rng: R, steps: usize) -> Polygon<f64> {
     let mut ring = Vec::with_capacity(2 * steps);
 
-    let ystep = 1.0;
-    let nudge_std = ystep / 1000.0;
+    let y_step = 10.0;
+    let nudge_std = y_step / 1000.0;
     let mut y = 0.0;
     let normal = Normal::new(0.0, nudge_std * nudge_std).unwrap();
-    let shift = 50.0;
+    let x_shift = 100.0;
 
     ring.push((0.0, 0.0).into());
     (0..steps).for_each(|_| {
-        let x: f64 = rng.sample::<f64, _>(Standard) * shift / 2.;
-        let x = (x * 10.) as i64 as f64 / 10.;
-        y += ystep;
-        // y += normal.sample(&mut rng);
+        let x: f64 = rng.sample::<f64, _>(Standard);
+        y += y_step;
         ring.push((x, y).into());
     });
-    ring.push((shift, y).into());
+    ring.push((x_shift, y).into());
     (0..steps).for_each(|_| {
-        let x: f64 = rng.sample::<f64, _>(Standard) * shift;
-        let x = (x * 10.) as i64 as f64 / 10.;
-        y -= ystep;
+        let x: f64 = rng.sample::<f64, _>(Standard);
+        y -= y_step;
         // y += normal.sample(&mut rng);
-        ring.push((shift + x, y).into());
+        ring.push((x_shift + x, y).into());
     });
 
-    Polygon::new(LineString(ring), vec![])
+    normalize_polygon(Polygon::new(LineString(ring), vec![]))
 }
 
+/// Normalizes polygon to fit and fill `[-1, 1] X [-1, 1]` square.
+///
+/// Uses `MapCoord` and `BoundingRect`
+pub fn normalize_polygon(poly: Polygon<f64>) -> Polygon<f64> {
+    let bounds = poly.bounding_rect().unwrap();
+    let dims = bounds.max() - bounds.min();
+    let x_scale = 2. / dims.x;
+    let y_scale = 2. / dims.y;
+
+    let x_shift = -bounds.min().x * x_scale - 1.;
+    let y_shift = -bounds.min().y * y_scale - 1.;
+    poly.map_coords(|mut c| {
+        c.x *= x_scale;
+        c.x += x_shift;
+        c.y *= y_scale;
+        c.y += y_shift;
+        c
+    })
+}
+
+#[derive(Debug, Clone)]
 pub struct Samples<T>(Vec<T>);
 impl<T> Samples<T> {
     pub fn sampler<'a>(&'a self) -> impl FnMut() -> &'a T {
@@ -116,4 +134,8 @@ impl<T> Samples<T> {
     pub fn from_fn<F: FnMut() -> T>(size: usize, mut proc: F) -> Self {
         Self((0..size).map(|_| proc()).collect())
     }
+
+    pub fn map<U, F: FnMut(T) -> U>(self, mut proc: F) -> Samples<U> {
+        Samples(self.0.into_iter().map(proc).collect())
+    }
 }

geo/src/algorithm/area.rs

@@ -526,6 +526,10 @@ mod test {
             ],
         ];
         // Value from shapely
-        assert_relative_eq!(poly.unsigned_area(), 0.006547948219252177, max_relative = 0.0001);
+        assert_relative_eq!(
+            poly.unsigned_area(),
+            0.006547948219252177,
+            max_relative = 0.0001
+        );
     }
 }

geo/src/algorithm/mod.rs

@@ -267,3 +267,7 @@ pub mod sweep;
 pub mod outlier_detection;
 
 pub use outlier_detection::OutlierDetection;
+
+/// Monotonic polygon subdivision
+pub mod monotone;
+pub use monotone::{monotone_subdivision, MonoPoly, MonotonicPolygons};