Merge #639

639: Introduce the geomgraph module for DE-9IM Relate trait r=michaelkirk a=michaelkirk

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

Fixes #513, #515

(I'm sorry it's so large)

~~I'm going to leave it as a draft (edit: 🤦 I failed to actually open the PR as a draft) while I wait to merge #636 and #638 and then do some rebasing, but I don't anticipate doing other large changes before review.~~ *update: ready for review!*

Here's some of the earlier work in pursuit of this:

#514
#516
#523
#524 
#538
#552
#561
#611
#628
#629
#636 

Primarily, this introduces the geomgraph module for a 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

You can see some of the reference code if you omit the "REMOVE JTS COMMENTS" commit. In some places the implementation is quite close to the JTS source. 

The overall "flow" is pretty similar to that of JTS, but in the small, there were some divergences. It's not easy (or desirable) to literally translate a Java codebase making heavy use of inheritance and pointers to rust. Additionally, I chose to take advantage of `Option` and rust's enums with associated data to make some case analysis more explicit.

There is a corresponding PR in our [jts-test-runner](georust/jts-test-runner#6) crate which includes the bulk of the tests for the new Relate trait.

## Algorithm Overview 

This functionality is accessed on geometries, via the `Relate` trait, e.g. `line.relate(point)` which returns a DE-9IM [`IntersectionMatrix`](https://en.wikipedia.org/wiki/DE-9IM#Matrix_model).

The `Relate` trait is driven by the `RelateOperation`. The `RelateOperation` builds a `GeometryGraph` for each of the two geometries being related. 

A `GeometryGraph` is a systematic way to organize the "interesting" parts of a geometry's structure - e.g. where its vertices, lines, and areas lie relative to one another.

Once the `RelateOperation` has built the two `GeometryGraph`s, it uses them to efficiently compare the two Geometries's structures, outputting the `IntesectionMatrix`.





Co-authored-by: Michael Kirk <michael.code@endoftheworl.de>
Co-authored-by: bors[bot] <26634292+bors[bot]@users.noreply.github.com>

geo-types/src/private_utils.rs

@@ -16,15 +16,7 @@ pub fn line_bounding_rect<T>(line: Line<T>) -> Rect<T>
 where
     T: CoordNum,
 {
-    let a = line.start;
-    let b = line.end;
-    let (xmin, xmax) = if a.x <= b.x { (a.x, b.x) } else { (b.x, a.x) };
-    let (ymin, ymax) = if a.y <= b.y { (a.y, b.y) } else { (b.y, a.y) };
-
-    Rect::new(
-        Coordinate { x: xmin, y: ymin },
-        Coordinate { x: xmax, y: ymax },
-    )
+    Rect::new(line.start, line.end)
 }
 
 pub fn get_bounding_rect<I, T>(collection: I) -> Option<Rect<T>>

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", branch = "mkirk/relate" }
 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};
 
@@ -66,14 +66,7 @@ where
     type Output = Rect<T>;
 
     fn bounding_rect(&self) -> Self::Output {
-        let a = self.start;
-        let b = self.end;
-        let (xmin, xmax) = if a.x <= b.x { (a.x, b.x) } else { (b.x, a.x) };
-        let (ymin, ymax) = if a.y <= b.y { (a.y, b.y) } else { (b.y, a.y) };
-        Rect::new(
-            Coordinate { x: xmin, y: ymin },
-            Coordinate { x: xmax, y: ymax },
-        )
+        Rect::new(self.start, self.end)
     }
 }
 
@@ -166,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() {