Merge branch 'main' into mkirk/im-matches

thehappycheese · Aug 7, 2023 · 4c64892 · 4c64892
2 parents 06aa146 + 1df7851
commit 4c64892
Show file tree

Hide file tree

Showing 19 changed files with 1,794 additions and 28 deletions.
diff --git a/geo/CHANGES.md b/geo/CHANGES.md
@@ -2,9 +2,10 @@
 
 ## Unreleased
 
+* Add `HausdorffDistance` algorithm trait to calculate the Hausdorff distance between any two geometries. <https://github.com/georust/geo/pull/1041>
+  * <https://github.com/georust/geo/pull/1041>
 * Add `matches` method to IntersectionMatrix for ergonomic de-9im comparisons.
-  <https://github.com/georust/geo/pull/1043>
-
+  * <https://github.com/georust/geo/pull/1043>
 
 ## 0.26.0
 
@@ -17,6 +18,13 @@
 
 - Add `TriangulateEarcut` algorithm trait to triangulate polygons with the earcut algorithm.
   - <https://github.com/georust/geo/pull/1007>
+- 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
 

diff --git a/geo/Cargo.toml b/geo/Cargo.toml
@@ -105,3 +105,7 @@ harness = false
 [[bench]]
 name = "winding_order"
 harness = false
+
+[[bench]]
+name = "monotone_subdiv"
+harness = false
diff --git a/geo/benches/monotone_subdiv.rs b/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);
diff --git a/geo/benches/utils/random.rs b/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())
+    }
 }
diff --git a/geo/src/algorithm/hausdorff_distance.rs b/geo/src/algorithm/hausdorff_distance.rs
@@ -0,0 +1,136 @@
+use crate::algorithm::EuclideanDistance;
+use crate::CoordsIter;
+use crate::GeoFloat;
+use geo_types::{Coord, Point};
+use num_traits::Bounded;
+
+/// Determine the distance between two geometries using the [Hausdorff distance formula].
+///
+/// Hausdorff distance is used to compare two point sets. It measures the maximum euclidean
+/// distance of a point in one set to the nearest point in another set. Hausdorff distance
+/// is often used to measure the amount of mismatch between two sets.
+///
+/// [Hausdorff distance formula]: https://en.wikipedia.org/wiki/Hausdorff_distance
+
+pub trait HausdorffDistance<T>
+where
+    T: GeoFloat,
+{
+    fn hausdorff_distance<Rhs>(&self, rhs: &Rhs) -> T
+    where
+        Rhs: for<'a> CoordsIter<'a, Scalar = T>;
+}
+
+impl<T, G> HausdorffDistance<T> for G
+where
+    T: GeoFloat,
+    G: for<'a> CoordsIter<'a, Scalar = T>,
+{
+    fn hausdorff_distance<Rhs>(&self, rhs: &Rhs) -> T
+    where
+        Rhs: for<'a> CoordsIter<'a, Scalar = T>,
+    {
+        // calculate from A -> B
+        let hd1 = self
+            .coords_iter()
+            .map(|c| {
+                rhs.coords_iter()
+                    .map(|c2| c.euclidean_distance(&c2))
+                    .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
+            })
+            .fold(<T as Bounded>::min_value(), |accum, val| accum.max(val));
+
+        // Calculate from B -> A
+        let hd2 = rhs
+            .coords_iter()
+            .map(|c| {
+                self.coords_iter()
+                    .map(|c2| c.euclidean_distance(&c2))
+                    .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
+            })
+            .fold(<T as Bounded>::min_value(), |accum, val| accum.max(val));
+
+        // The max of the two
+        hd1.max(hd2)
+    }
+}
+
+// ┌───────────────────────────┐
+// │ Implementations for Coord │
+// └───────────────────────────┘
+
+impl<T> HausdorffDistance<T> for Coord<T>
+where
+    T: GeoFloat,
+{
+    fn hausdorff_distance<Rhs>(&self, rhs: &Rhs) -> T
+    where
+        Rhs: for<'a> CoordsIter<'a, Scalar = T>,
+    {
+        Point::from(*self).hausdorff_distance(rhs)
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use crate::HausdorffDistance;
+    use crate::{line_string, polygon, MultiPoint, MultiPolygon};
+
+    #[test]
+    fn hd_mpnt_mpnt() {
+        let p1: MultiPoint<_> = vec![(0., 0.), (1., 2.)].into();
+        let p2: MultiPoint<_> = vec![(2., 3.), (1., 2.)].into();
+        assert_relative_eq!(p1.hausdorff_distance(&p2), 2.236068, epsilon = 1.0e-6);
+    }
+
+    #[test]
+    fn hd_mpnt_poly() {
+        let p1: MultiPoint<_> = vec![(0., 0.), (1., 2.)].into();
+        let poly = polygon![
+        (x: 1., y: -3.1), (x: 3.7, y: 2.7),
+        (x: 0.9, y: 7.6), (x: -4.8, y: 6.7),
+        (x: -7.5, y: 0.9), (x: -4.7, y: -4.),
+        (x: 1., y: -3.1)
+        ];
+
+        assert_relative_eq!(p1.hausdorff_distance(&poly), 7.553807, epsilon = 1.0e-6)
+    }
+
+    #[test]
+    fn hd_mpnt_lns() {
+        let p1: MultiPoint<_> = vec![(0., 0.), (1., 2.)].into();
+        let lns = line_string![
+        (x: 1., y: -3.1), (x: 3.7, y: 2.7),
+        (x: 0.9, y: 7.6), (x: -4.8, y: 6.7),
+        (x: -7.5, y: 0.9), (x: -4.7, y: -4.),
+        (x: 1., y: -3.1)
+        ];
+
+        assert_relative_eq!(p1.hausdorff_distance(&lns), 7.553807, epsilon = 1.0e-6)
+    }
+
+    #[test]
+    fn hd_mpnt_mply() {
+        let p1: MultiPoint<_> = vec![(0., 0.), (1., 2.)].into();
+        let multi_polygon = MultiPolygon::new(vec![
+            polygon![
+              (x: 0.0f32, y: 0.0),
+              (x: 2.0, y: 0.0),
+              (x: 2.0, y: 1.0),
+              (x: 0.0, y: 1.0),
+            ],
+            polygon![
+              (x: 1.0, y: 1.0),
+              (x: -2.0, y: 1.0),
+              (x: -2.0, y: -1.0),
+              (x: 1.0, y: -1.0),
+            ],
+        ]);
+
+        assert_relative_eq!(
+            p1.hausdorff_distance(&multi_polygon),
+            2.236068,
+            epsilon = 1.0e-6
+        )
+    }
+}