Introduce the geomgraph module for DE-9IM Relate trait

geomgraph implements a topology graph largely inspired by JTS's module
of the same name:
https://github.com/locationtech/jts/tree/jts-1.18.1/modules/core/src/main/java/org/locationtech/jts/geomgraph

geo/CHANGES.md

@@ -4,6 +4,12 @@
 
 * Add `line_intersection` to compute point or segment intersection of two Lines.
   * <https://github.com/georust/geo/pull/636>
+* Add `Relate` trait to topologically relate two geometries based on [DE-9IM](https://en.wikipedia.org/wiki/DE-9IM) semantics.
+  * <https://github.com/georust/geo/pull/639>
+* Fix `Contains` implementation for Polygons to match the OGC spec using the new `Relate` trait
+  * <https://github.com/georust/geo/pull/639>
+* BREAKING: `Contains` no longer supports integer `Polygon` and `Geometry`
+  * <https://github.com/georust/geo/pull/639>
 
 ## 0.17.1
 

geo/Cargo.toml

@@ -19,6 +19,7 @@ num-traits = "0.2"
 serde = { version = "1.0", optional = true, features = ["derive"] }
 rstar = { version = "0.8" }
 geographiclib-rs = { version = "0.2" }
+log = "0.4.11"
 
 proj = { version = "0.20.3", optional = true }
 
@@ -32,11 +33,12 @@ proj-network = ["use-proj", "proj/network"]
 use-serde = ["serde", "geo-types/serde"]
 
 [dev-dependencies]
+pretty_env_logger = "0.4"
 approx = "0.4.0"
 criterion = { version = "0.3" }
 # jts-test-runner is an internal crate which exists only to be part of the geo test suite.
 # As such it's kept unpublished. It's in a separate repo primarily because it's kind of large.
-jts-test-runner = { git = "https://github.com/georust/jts-test-runner", rev = "3294b9af1d3e64fcc9caf9646a93d0116f7d8321" }
+jts-test-runner = { git = "https://github.com/georust/jts-test-runner", rev = "4a69a55" }
 rand = "0.8.0"
 
 [[bench]]
@@ -75,6 +77,10 @@ harness = false
 name = "rotate"
 harness = false
 
+[[bench]]
+name = "relate"
+harness = false
+
 [[bench]]
 name = "simplify"
 harness = false

geo/benches/relate.rs

@@ -0,0 +1,36 @@
+#[macro_use]
+extern crate criterion;
+extern crate geo;
+
+use crate::geo::relate::Relate;
+use criterion::Criterion;
+use geo::{LineString, Polygon};
+
+fn criterion_benchmark(c: &mut Criterion) {
+    c.bench_function("relate overlapping 50-point polygons", |bencher| {
+        let points = include!("../src/algorithm/test_fixtures/norway_main.rs");
+
+        let sub_polygon = {
+            let points = points[0..50].to_vec();
+            let mut exterior = LineString::<f32>::from(points);
+            exterior.close();
+            Polygon::new(exterior, vec![])
+        };
+
+        let polygon = {
+            let points = points[40..90].to_vec();
+            let mut exterior = LineString::<f32>::from(points);
+            exterior.close();
+            Polygon::new(exterior, vec![])
+        };
+
+        bencher.iter(|| {
+            criterion::black_box(
+                criterion::black_box(&polygon).relate(criterion::black_box(&sub_polygon)),
+            );
+        });
+    });
+}
+
+criterion_group!(benches, criterion_benchmark);
+criterion_main!(benches);

geo/src/algorithm/bounding_rect.rs

@@ -1,7 +1,7 @@
 use crate::utils::{partial_max, partial_min};
 use crate::{
-    CoordNum, Coordinate, Geometry, GeometryCollection, Line, LineString, MultiLineString,
-    MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
+    CoordNum, Coordinate, Geometry, GeometryCollection, GeometryCow, Line, LineString,
+    MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
 };
 use geo_types::private_utils::{get_bounding_rect, line_string_bounding_rect};
 
@@ -159,6 +159,17 @@ where
     }
 }
 
+impl<T> BoundingRect<T> for GeometryCow<'_, T>
+where
+    T: CoordNum,
+{
+    type Output = Option<Rect<T>>;
+
+    crate::geometry_cow_delegate_impl! {
+       fn bounding_rect(&self) -> Self::Output;
+    }
+}
+
 impl<T> BoundingRect<T> for GeometryCollection<T>
 where
     T: CoordNum,

geo/src/algorithm/centroid.rs

@@ -615,8 +615,8 @@ mod test {
             .centroid()
             .unwrap()
             .map_coords(|&(x, y)| (x - shift_x, y - shift_y));
-        eprintln!("centroid {:?}", centroid.0);
-        eprintln!("new_centroid {:?}", new_centroid.0);
+        debug!("centroid {:?}", centroid.0);
+        debug!("new_centroid {:?}", new_centroid.0);
         assert_relative_eq!(centroid.0.x, new_centroid.0.x, max_relative = 0.0001);
         assert_relative_eq!(centroid.0.y, new_centroid.0.y, max_relative = 0.0001);
     }

geo/src/algorithm/contains/geometry.rs

@@ -9,17 +9,17 @@ impl<T> Contains<Coordinate<T>> for Geometry<T>
 where
     T: GeoNum,
 {
-    geometry_delegate_impl! {
-        fn contains(&self, coord: &Coordinate<T>) -> bool;
+    fn contains(&self, coord: &Coordinate<T>) -> bool {
+        self.contains(&Point::from(*coord))
     }
 }
 
 impl<T> Contains<Point<T>> for Geometry<T>
 where
     T: GeoNum,
 {
-    fn contains(&self, point: &Point<T>) -> bool {
-        self.contains(&point.0)
+    geometry_delegate_impl! {
+        fn contains(&self, point: &Point<T>) -> bool;
     }
 }
 

geo/src/algorithm/contains/mod.rs

@@ -370,29 +370,45 @@ mod test {
     #[test]
     fn line_in_polygon_test() {
         let c = |x, y| Coordinate { x, y };
-        let line = Line::new(c(0, 10), c(30, 40));
-        let linestring0 = line_string![c(-10, 0), c(50, 0), c(50, 50), c(0, 50), c(-10, 0)];
+        let line = Line::new(c(0.0, 10.0), c(30.0, 40.0));
+        let linestring0 = line_string![
+            c(-10.0, 0.0),
+            c(50.0, 0.0),
+            c(50.0, 50.0),
+            c(0.0, 50.0),
+            c(-10.0, 0.0)
+        ];
         let poly0 = Polygon::new(linestring0, Vec::new());
-        let linestring1 = line_string![c(0, 0), c(0, 20), c(20, 20), c(20, 0), c(0, 0)];
+        let linestring1 = line_string![
+            c(0.0, 0.0),
+            c(0.0, 20.0),
+            c(20.0, 20.0),
+            c(20.0, 0.0),
+            c(0.0, 0.0)
+        ];
         let poly1 = Polygon::new(linestring1, Vec::new());
         assert!(poly0.contains(&line));
         assert!(!poly1.contains(&line));
     }
     #[test]
-    #[ignore]
     fn line_in_polygon_edgecases_test() {
         // Some DE-9IM edge cases for checking line is
         // inside polygon The end points of the line can be
-        // on the boundary of the polygon but we don't allow
-        // that yet.
+        // on the boundary of the polygon.
         let c = |x, y| Coordinate { x, y };
         // A non-convex polygon
-        let linestring0 = line_string![c(0, 0), c(1, 1), c(1, -1), c(-1, -1), c(-1, 1)];
+        let linestring0 = line_string![
+            c(0.0, 0.0),
+            c(1.0, 1.0),
+            c(1.0, -1.0),
+            c(-1.0, -1.0),
+            c(-1.0, 1.0)
+        ];
         let poly = Polygon::new(linestring0, Vec::new());
 
-        assert!(poly.contains(&Line::new(c(0, 0), c(1, -1))));
-        assert!(poly.contains(&Line::new(c(-1, 1), c(1, -1))));
-        assert!(!poly.contains(&Line::new(c(-1, 1), c(1, 1))));
+        assert!(poly.contains(&Line::new(c(0.0, 0.0), c(1.0, -1.0))));
+        assert!(poly.contains(&Line::new(c(-1.0, 1.0), c(1.0, -1.0))));
+        assert!(!poly.contains(&Line::new(c(-1.0, 1.0), c(1.0, 1.0))));
     }
     #[test]
     fn line_in_linestring_edgecases() {

geo/src/algorithm/contains/polygon.rs

@@ -1,11 +1,13 @@
 use super::Contains;
-use crate::intersects::Intersects;
-use crate::{CoordNum, Coordinate, GeoNum, Line, LineString, MultiPolygon, Point, Polygon};
+use crate::relate::Relate;
+use crate::{
+    Coordinate, GeoFloat, GeoNum, GeometryCollection, Line, LineString, MultiLineString,
+    MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
+};
 
 // ┌─────────────────────────────┐
 // │ Implementations for Polygon │
 // └─────────────────────────────┘
-
 impl<T> Contains<Coordinate<T>> for Polygon<T>
 where
     T: GeoNum,
@@ -26,63 +28,129 @@ where
     }
 }
 
-// TODO: ensure DE-9IM compliance: esp., when
-// line.start and line.end is on the boundaries
 impl<T> Contains<Line<T>> for Polygon<T>
 where
-    T: GeoNum,
+    T: GeoFloat,
 {
     fn contains(&self, line: &Line<T>) -> bool {
-        // both endpoints are contained in the polygon and the line
-        // does NOT intersect the exterior or any of the interior boundaries
-        self.contains(&line.start)
-            && self.contains(&line.end)
-            && !self.exterior().intersects(line)
-            && !self.interiors().iter().any(|inner| inner.intersects(line))
+        self.relate(line).is_contains()
     }
 }
 
-// TODO: also check interiors
 impl<T> Contains<Polygon<T>> for Polygon<T>
 where
-    T: GeoNum,
+    T: GeoFloat,
 {
     fn contains(&self, poly: &Polygon<T>) -> bool {
-        // decompose poly's exterior ring into Lines, and check each for containment
-        poly.exterior().lines().all(|line| self.contains(&line))
+        self.relate(poly).is_contains()
     }
 }
 
-// TODO: ensure DE-9IM compliance
 impl<T> Contains<LineString<T>> for Polygon<T>
 where
-    T: GeoNum,
+    T: GeoFloat,
 {
     fn contains(&self, linestring: &LineString<T>) -> bool {
-        // All LineString points must be inside the Polygon
-        if linestring.points_iter().all(|point| self.contains(&point)) {
-            // The Polygon interior is allowed to intersect with the LineString
-            // but the Polygon's rings are not
-            !self
-                .interiors()
-                .iter()
-                .any(|ring| ring.intersects(linestring))
-        } else {
-            false
-        }
+        self.relate(linestring).is_contains()
     }
 }
 
 // ┌──────────────────────────────────┐
 // │ Implementations for MultiPolygon │
 // └──────────────────────────────────┘
-// TODO: ensure DE-9IM compliance
-impl<G, T> Contains<G> for MultiPolygon<T>
+
+impl<T> Contains<Coordinate<T>> for MultiPolygon<T>
+where
+    T: GeoNum,
+{
+    fn contains(&self, coord: &Coordinate<T>) -> bool {
+        self.iter().any(|poly| poly.contains(coord))
+    }
+}
+
+impl<T> Contains<Point<T>> for MultiPolygon<T>
+where
+    T: GeoNum,
+{
+    fn contains(&self, p: &Point<T>) -> bool {
+        self.contains(&p.0)
+    }
+}
+
+impl<T: GeoNum> Contains<MultiPoint<T>> for MultiPolygon<T> {
+    fn contains(&self, rhs: &MultiPoint<T>) -> bool {
+        rhs.iter().all(|point| self.contains(point))
+    }
+}
+
+impl<F> Contains<Line<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &Line<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<LineString<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &LineString<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<MultiLineString<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &MultiLineString<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<Polygon<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &Polygon<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<MultiPolygon<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &MultiPolygon<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<GeometryCollection<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &GeometryCollection<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<Rect<F>> for MultiPolygon<F>
+where
+    F: GeoFloat,
+{
+    fn contains(&self, rhs: &Rect<F>) -> bool {
+        rhs.relate(self).is_within()
+    }
+}
+
+impl<F> Contains<Triangle<F>> for MultiPolygon<F>
 where
-    T: CoordNum,
-    Polygon<T>: Contains<G>,
+    F: GeoFloat,
 {
-    fn contains(&self, rhs: &G) -> bool {
-        self.iter().any(|p| p.contains(rhs))
+    fn contains(&self, rhs: &Triangle<F>) -> bool {
+        rhs.relate(self).is_within()
     }
 }