From 207226b802d7957f02f47f6bf5d3a947facda134 Mon Sep 17 00:00:00 2001
From: Michael Kirk <michael.code@endoftheworl.de>
Date: Thu, 5 May 2022 11:22:05 -0700
Subject: [PATCH] working rtree intersector impl

---
 geo/CHANGES.md                                |  4 +
 .../index/rstar_edge_set_intersector.rs       | 90 ++++++++++++++-----
 2 files changed, 72 insertions(+), 22 deletions(-)

diff --git a/geo/CHANGES.md b/geo/CHANGES.md
index 86e4673594..741c4d227c 100644
--- a/geo/CHANGES.md
+++ b/geo/CHANGES.md
@@ -3,6 +3,10 @@
 
 ## Unreleased
 
+* POSSIBLY BREAKING: `GeoFloat` types must now implement `num_traits::Signed` and `num_traits::Bounded`. This shouldn't 
+  affect you if you are using a standard `Geometry<f64>` or `Geometry<f32>` or `geo::GeoFloat` generically.
+* Speed up `Relate` and `Contains` traits for large LineStrings and Polygons by using an RTree to more efficiently
+  inspect edges in our topology graph.
 * Speed up intersection checks by using a preliminary bbox check
 * Remove unneeded reference for `*MapCoords*` closure parameter.
   * <https://github.com/georust/geo/pull/810>
diff --git a/geo/src/algorithm/relate/geomgraph/index/rstar_edge_set_intersector.rs b/geo/src/algorithm/relate/geomgraph/index/rstar_edge_set_intersector.rs
index 91619f9197..6bb180714a 100644
--- a/geo/src/algorithm/relate/geomgraph/index/rstar_edge_set_intersector.rs
+++ b/geo/src/algorithm/relate/geomgraph/index/rstar_edge_set_intersector.rs
@@ -5,41 +5,71 @@ use crate::GeoFloat;
 use std::cell::RefCell;
 use std::rc::Rc;
 
+use rstar::RTree;
+
 pub(crate) struct RstarEdgeSetIntersector;
 
 impl RstarEdgeSetIntersector {
     pub fn new() -> Self {
         RstarEdgeSetIntersector
     }
+}
 
-    fn compute_intersects<F: GeoFloat>(
-        &mut self,
-        edge0: &Rc<RefCell<Edge<F>>>,
-        edge1: &Rc<RefCell<Edge<F>>>,
-        segment_intersector: &mut SegmentIntersector<F>,
-    ) {
-        let edge0_coords_len = edge0.borrow().coords().len() - 1;
-        let edge1_coords_len = edge1.borrow().coords().len() - 1;
-        for i0 in 0..edge0_coords_len {
-            for i1 in 0..edge1_coords_len {
-                segment_intersector.add_intersections(edge0, i0, edge1, i1);
-            }
+struct Segment<'a, F: GeoFloat + rstar::RTreeNum> {
+    i: usize,
+    edge: &'a RefCell<Edge<F>>,
+    envelope: rstar::AABB<crate::Coordinate<F>>,
+}
+
+impl<'a, F> Segment<'a, F>
+where
+    F: GeoFloat + rstar::RTreeNum,
+{
+    fn new(i: usize, edge: &'a RefCell<Edge<F>>) -> Self {
+        use crate::rstar::RTreeObject;
+        let p1 = edge.borrow().coords()[i];
+        let p2 = edge.borrow().coords()[i + 1];
+        Self {
+            i,
+            edge,
+            envelope: rstar::AABB::from_corners(p1, p2),
         }
     }
 }
 
-impl<F: GeoFloat> EdgeSetIntersector<F> for RstarEdgeSetIntersector {
+impl<'a, F> rstar::RTreeObject for Segment<'a, F>
+where
+    F: GeoFloat + rstar::RTreeNum,
+{
+    type Envelope = rstar::AABB<crate::Coordinate<F>>;
+
+    fn envelope(&self) -> Self::Envelope {
+        self.envelope
+    }
+}
+
+impl<F> EdgeSetIntersector<F> for RstarEdgeSetIntersector
+where
+    F: GeoFloat + rstar::RTreeNum,
+{
     fn compute_intersections_within_set(
         &mut self,
         edges: &[Rc<RefCell<Edge<F>>>],
         check_for_self_intersecting_edges: bool,
         segment_intersector: &mut SegmentIntersector<F>,
     ) {
-        for edge0 in edges.iter() {
-            for edge1 in edges.iter() {
-                if check_for_self_intersecting_edges || edge0.as_ptr() != edge1.as_ptr() {
-                    self.compute_intersects(edge0, edge1, segment_intersector);
-                }
+        let segments: Vec<Segment<F>> = edges
+            .iter()
+            .flat_map(|edge| {
+                let start_of_final_segment: usize = RefCell::borrow(edge).coords().len() - 1;
+                (0..start_of_final_segment).map(|segment_i| Segment::new(segment_i, edge))
+            })
+            .collect();
+        let tree = RTree::bulk_load(segments);
+
+        for (edge0, edge1) in tree.intersection_candidates_with_other_tree(&tree) {
+            if check_for_self_intersecting_edges || edge0.edge.as_ptr() != edge1.edge.as_ptr() {
+                segment_intersector.add_intersections(edge0.edge, edge0.i, edge1.edge, edge1.i);
             }
         }
     }
@@ -50,10 +80,26 @@ impl<F: GeoFloat> EdgeSetIntersector<F> for RstarEdgeSetIntersector {
         edges1: &[Rc<RefCell<Edge<F>>>],
         segment_intersector: &mut SegmentIntersector<F>,
     ) {
-        for edge0 in edges0 {
-            for edge1 in edges1 {
-                self.compute_intersects(edge0, edge1, segment_intersector);
-            }
+        let segments0: Vec<Segment<F>> = edges0
+            .iter()
+            .flat_map(|edge| {
+                let start_of_final_segment: usize = RefCell::borrow(edge).coords().len() - 1;
+                (0..start_of_final_segment).map(|segment_i| Segment::new(segment_i, edge))
+            })
+            .collect();
+        let tree_0 = RTree::bulk_load(segments0);
+
+        let segments1: Vec<Segment<F>> = edges1
+            .iter()
+            .flat_map(|edge| {
+                let start_of_final_segment: usize = RefCell::borrow(edge).coords().len() - 1;
+                (0..start_of_final_segment).map(|segment_i| Segment::new(segment_i, edge))
+            })
+            .collect();
+        let tree_1 = RTree::bulk_load(segments1);
+
+        for (edge0, edge1) in tree_0.intersection_candidates_with_other_tree(&tree_1) {
+            segment_intersector.add_intersections(edge0.edge, edge0.i, edge1.edge, edge1.i);
         }
     }
 }