Add PreparedGeometry to speed up repeated Relate operations.

Much of the cost of the `Relate` operation comes from finding all
intersections between all segments of the two geometries.

To speed this up, the edges of each geometry are put into an R-Tree.

We also have to "self node" each geometry, meaning we need to split any
segments at the point that they'd intersect with an edge from the other
geometry.

None of that work is new to this commit, but what is new is that we now
cache the self-noding work and the geometry's r-tree.

    relate prepared polygons
                            time:   [49.036 ms 49.199 ms 49.373 ms]
    relate unprepared polygons
                            time:   [842.91 ms 844.32 ms 845.82 ms]

geo/CHANGES.md

@@ -10,6 +10,8 @@
   * <https://github.com/georust/geo/pull/1192>
 * Fix `AffineTransform::compose` ordering to be conventional - such that the argument is applied *after* self.
   * <https://github.com/georust/geo/pull/1196>
+* Add `PreparedGeometry` to speed up repeated `Relate` operations.
+  * <https://github.com/georust/geo/pull/1197>
 
 ## 0.28.0
 

geo/benches/prepared_geometry.rs

@@ -1,11 +1,45 @@
 use criterion::{criterion_group, criterion_main, Criterion};
 use geo::algorithm::Relate;
-use geo::geometry::MultiPolygon;
+use geo::PreparedGeometry;
+use geo_types::MultiPolygon;
 
 fn criterion_benchmark(c: &mut Criterion) {
+    c.bench_function("relate prepared polygons", |bencher| {
+        let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots();
+        let zone_polygons = geo_test_fixtures::nl_zones();
+
+        bencher.iter(|| {
+            let mut intersects = 0;
+            let mut non_intersects = 0;
+
+            let plot_polygons = plot_polygons
+                .iter()
+                .map(PreparedGeometry::from)
+                .collect::<Vec<_>>();
+
+            let zone_polygons = zone_polygons
+                .iter()
+                .map(PreparedGeometry::from)
+                .collect::<Vec<_>>();
+
+            for a in &plot_polygons {
+                for b in &zone_polygons {
+                    if criterion::black_box(a.relate(b).is_intersects()) {
+                        intersects += 1;
+                    } else {
+                        non_intersects += 1;
+                    }
+                }
+            }
+
+            assert_eq!(intersects, 974);
+            assert_eq!(non_intersects, 27782);
+        });
+    });
+
     c.bench_function("relate unprepared polygons", |bencher| {
         let plot_polygons: MultiPolygon = geo_test_fixtures::nl_plots();
-        let zone_polygons: MultiPolygon = geo_test_fixtures::nl_zones();
+        let zone_polygons = geo_test_fixtures::nl_zones();
 
         bencher.iter(|| {
             let mut intersects = 0;

geo/src/algorithm/relate/geomgraph/edge.rs

@@ -7,7 +7,7 @@ use std::collections::BTreeSet;
 /// An `Edge` represents a one dimensional line in a geometry.
 ///
 /// This is based on [JTS's `Edge` as of 1.18.1](https://github.com/locationtech/jts/blob/jts-1.18.1/modules/core/src/main/java/org/locationtech/jts/geomgraph/Edge.java)
-#[derive(Debug)]
+#[derive(Debug, PartialEq, Clone)]
 pub(crate) struct Edge<F: GeoFloat> {
     /// `coordinates` of the line geometry
     coords: Vec<Coord<F>>,
@@ -48,13 +48,21 @@ impl<F: GeoFloat> Edge<F> {
         &mut self.label
     }
 
+    /// When comparing two prepared geometries, we cache each geometry's topology graph. Depending
+    /// on the order of the operation - `a.relate(b)` vs `b.relate(a)` - we may need to swap the
+    /// label.
+    pub fn swap_label_args(&mut self) {
+        self.label.swap_args()
+    }
+
     pub fn coords(&self) -> &[Coord<F>] {
         &self.coords
     }
 
     pub fn is_isolated(&self) -> bool {
         self.is_isolated
     }
+
     pub fn mark_as_unisolated(&mut self) {
         self.is_isolated = false;
     }

geo/src/algorithm/relate/geomgraph/edge_intersection.rs

@@ -6,7 +6,7 @@ use crate::{Coord, GeoFloat};
 /// the start of the line segment) The intersection point must be precise.
 ///
 /// This is based on [JTS's EdgeIntersection as of 1.18.1](https://github.com/locationtech/jts/blob/jts-1.18.1/modules/core/src/main/java/org/locationtech/jts/geomgraph/EdgeIntersection.java)
-#[derive(Debug)]
+#[derive(Debug, Clone)]
 pub(crate) struct EdgeIntersection<F: GeoFloat> {
     coord: Coord<F>,
     segment_index: usize,

geo/src/algorithm/relate/geomgraph/geometry_graph.rs

@@ -1,38 +1,43 @@
 use super::{
     index::{
-        EdgeSetIntersector, RstarEdgeSetIntersector, SegmentIntersector, SimpleEdgeSetIntersector,
+        EdgeSetIntersector, RStarEdgeSetIntersector, Segment, SegmentIntersector,
+        SimpleEdgeSetIntersector,
     },
     CoordNode, CoordPos, Direction, Edge, Label, LineIntersector, PlanarGraph, TopologyPosition,
 };
 
 use crate::HasDimensions;
 use crate::{Coord, GeoFloat, GeometryCow, Line, LineString, Point, Polygon};
 
+use rstar::{RTree, RTreeNum};
 use std::cell::RefCell;
 use std::rc::Rc;
 
-/// The computation of the [`IntersectionMatrix`] relies on the use of a
-/// structure called a "topology graph". The topology graph contains [nodes](CoordNode) and
-/// [edges](Edge) corresponding to the nodes and line segments of a [`Geometry`]. Each
+/// The computation of the [`IntersectionMatrix`](crate::algorithm::relate::IntersectionMatrix) relies on the use of a
+/// structure called a "topology graph". The topology graph contains nodes (CoordNode) and
+/// edges (Edge) corresponding to the nodes and line segments of a [`Geometry`](crate::Geometry). Each
 /// node and edge in the graph is labeled with its topological location
 /// relative to the source geometry.
 ///
 /// Note that there is no requirement that points of self-intersection be a
-/// vertex. Thus to obtain a correct topology graph, [`Geometry`] must be
+/// vertex. Thus, to obtain a correct topology graph, [`Geometry`](crate::Geometry) must be
 /// self-noded before constructing their graphs.
 ///
 /// Two fundamental operations are supported by topology graphs:
 ///   - Computing the intersections between all the edges and nodes of a single graph
 ///   - Computing the intersections between the edges and nodes of two different graphs
 ///
 /// GeometryGraph is based on [JTS's `GeomGraph` as of 1.18.1](https://github.com/locationtech/jts/blob/jts-1.18.1/modules/core/src/main/java/org/locationtech/jts/geomgraph/GeometryGraph.java)
-pub(crate) struct GeometryGraph<'a, F>
+#[derive(Clone)]
+pub struct GeometryGraph<'a, F>
 where
     F: GeoFloat,
 {
     arg_index: usize,
-    parent_geometry: &'a GeometryCow<'a, F>,
+    parent_geometry: GeometryCow<'a, F>,
+    tree: Option<Rc<RTree<Segment<F>>>>,
     use_boundary_determination_rule: bool,
+    has_computed_self_nodes: bool,
     planar_graph: PlanarGraph<F>,
 }
 
@@ -44,44 +49,102 @@ impl<F> GeometryGraph<'_, F>
 where
     F: GeoFloat,
 {
-    pub fn edges(&self) -> &[Rc<RefCell<Edge<F>>>] {
+    pub(crate) fn set_tree(&mut self, tree: Rc<RTree<Segment<F>>>) {
+        self.tree = Some(tree);
+    }
+
+    pub(crate) fn get_or_build_tree(&self) -> Rc<RTree<Segment<F>>> {
+        self.tree
+            .clone()
+            .unwrap_or_else(|| Rc::new(self.build_tree()))
+    }
+
+    pub(crate) fn build_tree(&self) -> RTree<Segment<F>> {
+        let segments: Vec<Segment<F>> = self
+            .edges()
+            .iter()
+            .enumerate()
+            .flat_map(|(edge_idx, edge)| {
+                let edge = RefCell::borrow(edge);
+                let start_of_final_segment: usize = edge.coords().len() - 1;
+                (0..start_of_final_segment).map(move |segment_idx| {
+                    let p1 = edge.coords()[segment_idx];
+                    let p2 = edge.coords()[segment_idx + 1];
+                    Segment::new(edge_idx, segment_idx, p1, p2)
+                })
+            })
+            .collect();
+        RTree::bulk_load(segments)
+    }
+
+    pub(crate) fn assert_eq_graph(&self, other: &Self) {
+        assert_eq!(self.arg_index, other.arg_index);
+        assert_eq!(
+            self.use_boundary_determination_rule,
+            other.use_boundary_determination_rule
+        );
+        assert_eq!(self.parent_geometry, other.parent_geometry);
+        self.planar_graph.assert_eq_graph(&other.planar_graph);
+    }
+
+    pub(crate) fn clone_for_arg_index(&self, arg_index: usize) -> Self {
+        debug_assert!(
+            self.has_computed_self_nodes,
+            "should only be called after computing self nodes"
+        );
+        let planar_graph = self
+            .planar_graph
+            .clone_for_arg_index(self.arg_index, arg_index);
+        Self {
+            arg_index,
+            parent_geometry: self.parent_geometry.clone(),
+            tree: self.tree.clone(),
+            use_boundary_determination_rule: self.use_boundary_determination_rule,
+            has_computed_self_nodes: true,
+            planar_graph,
+        }
+    }
+
+    pub(crate) fn edges(&self) -> &[Rc<RefCell<Edge<F>>>] {
         self.planar_graph.edges()
     }
 
-    pub fn insert_edge(&mut self, edge: Edge<F>) {
+    pub(crate) fn insert_edge(&mut self, edge: Edge<F>) {
         self.planar_graph.insert_edge(edge)
     }
 
-    pub fn is_boundary_node(&self, coord: Coord<F>) -> bool {
+    pub(crate) fn is_boundary_node(&self, coord: Coord<F>) -> bool {
         self.planar_graph.is_boundary_node(self.arg_index, coord)
     }
 
-    pub fn add_node_with_coordinate(&mut self, coord: Coord<F>) -> &mut CoordNode<F> {
+    pub(crate) fn add_node_with_coordinate(&mut self, coord: Coord<F>) -> &mut CoordNode<F> {
         self.planar_graph.add_node_with_coordinate(coord)
     }
 
-    pub fn nodes_iter(&self) -> impl Iterator<Item = &CoordNode<F>> {
+    pub(crate) fn nodes_iter(&self) -> impl Iterator<Item = &CoordNode<F>> {
         self.planar_graph.nodes.iter()
     }
 }
 
 impl<'a, F> GeometryGraph<'a, F>
 where
-    F: GeoFloat,
+    F: GeoFloat + RTreeNum,
 {
-    pub fn new(arg_index: usize, parent_geometry: &'a GeometryCow<F>) -> Self {
+    pub(crate) fn new(arg_index: usize, parent_geometry: GeometryCow<'a, F>) -> Self {
         let mut graph = GeometryGraph {
             arg_index,
             parent_geometry,
             use_boundary_determination_rule: true,
+            tree: None,
+            has_computed_self_nodes: false,
             planar_graph: PlanarGraph::new(),
         };
-        graph.add_geometry(parent_geometry);
+        graph.add_geometry(&graph.parent_geometry.clone());
         graph
     }
 
-    pub fn geometry(&self) -> &GeometryCow<F> {
-        self.parent_geometry
+    pub(crate) fn geometry(&self) -> &GeometryCow<F> {
+        &self.parent_geometry
     }
 
     /// Determine whether a component (node or edge) that appears multiple times in elements
@@ -96,22 +159,11 @@ where
         }
     }
 
-    fn create_edge_set_intersector() -> Box<dyn EdgeSetIntersector<F>> {
-        // PERF: faster algorithms exist. This one was chosen for simplicity of implementation and
-        //       debugging
-        // Slow, but simple and good for debugging
-        // Box::new(SimpleEdgeSetIntersector::new())
-
-        // Should be much faster for sparse intersections, while not much slower than
-        // SimpleEdgeSetIntersector in the dense case
-        Box::new(RstarEdgeSetIntersector::new())
-    }
-
     fn boundary_nodes(&self) -> impl Iterator<Item = &CoordNode<F>> {
         self.planar_graph.boundary_nodes(self.arg_index)
     }
 
-    pub fn add_geometry(&mut self, geometry: &GeometryCow<F>) {
+    pub(crate) fn add_geometry(&mut self, geometry: &GeometryCow<F>) {
         if geometry.is_empty() {
             return;
         }
@@ -273,13 +325,13 @@ where
     /// assumed to be valid).
     ///
     /// `line_intersector` the [`LineIntersector`] to use to determine intersection
-    pub fn compute_self_nodes(
-        &mut self,
-        line_intersector: Box<dyn LineIntersector<F>>,
-    ) -> SegmentIntersector<F> {
-        let mut segment_intersector = SegmentIntersector::new(line_intersector, true);
+    pub(crate) fn compute_self_nodes(&mut self, line_intersector: Box<dyn LineIntersector<F>>) {
+        if self.has_computed_self_nodes {
+            return;
+        }
+        self.has_computed_self_nodes = true;
 
-        let mut edge_set_intersector = Self::create_edge_set_intersector();
+        let mut segment_intersector = SegmentIntersector::new(line_intersector, true);
 
         // optimize intersection search for valid Polygons and LinearRings
         let is_rings = match self.geometry() {
@@ -290,18 +342,16 @@ where
         };
         let check_for_self_intersecting_edges = !is_rings;
 
+        let edge_set_intersector = RStarEdgeSetIntersector;
         edge_set_intersector.compute_intersections_within_set(
-            self.edges(),
+            self,
             check_for_self_intersecting_edges,
             &mut segment_intersector,
         );
-
         self.add_self_intersection_nodes();
-
-        segment_intersector
     }
 
-    pub fn compute_edge_intersections(
+    pub(crate) fn compute_edge_intersections(
         &self,
         other: &GeometryGraph<F>,
         line_intersector: Box<dyn LineIntersector<F>>,
@@ -312,10 +362,10 @@ where
             other.boundary_nodes().cloned().collect(),
         );
 
-        let mut edge_set_intersector = Self::create_edge_set_intersector();
+        let edge_set_intersector = RStarEdgeSetIntersector;
         edge_set_intersector.compute_intersections_between_sets(
-            self.edges(),
-            other.edges(),
+            self,
+            other,
             &mut segment_intersector,
         );
 

geo/src/algorithm/relate/geomgraph/index/edge_set_intersector.rs

@@ -1,4 +1,4 @@
-use super::super::Edge;
+use super::super::{Edge, GeometryGraph};
 use super::SegmentIntersector;
 use crate::{Coord, GeoFloat};
 
@@ -13,18 +13,18 @@ pub(crate) trait EdgeSetIntersector<F: GeoFloat> {
     /// `check_for_self_intersecting_edges`: if false, an edge is not checked for intersections with itself.
     /// `segment_intersector`: the SegmentIntersector to use
     fn compute_intersections_within_set(
-        &mut self,
-        edges: &[Rc<RefCell<Edge<F>>>],
+        &self,
+        graph: &GeometryGraph<F>,
         check_for_self_intersecting_edges: bool,
         segment_intersector: &mut SegmentIntersector<F>,
     );
 
     /// Compute all intersections between two sets of edges, recording those intersections on
     /// the intersecting edges.
-    fn compute_intersections_between_sets(
-        &mut self,
-        edges0: &[Rc<RefCell<Edge<F>>>],
-        edges1: &[Rc<RefCell<Edge<F>>>],
+    fn compute_intersections_between_sets<'a>(
+        &self,
+        graph_0: &GeometryGraph<'a, F>,
+        graph_1: &GeometryGraph<'a, F>,
         segment_intersector: &mut SegmentIntersector<F>,
     );
 }

geo/src/algorithm/relate/geomgraph/index/mod.rs

@@ -1,9 +1,13 @@
 mod edge_set_intersector;
+mod prepared_geometry;
 mod rstar_edge_set_intersector;
+mod segment;
 mod segment_intersector;
 mod simple_edge_set_intersector;
 
 pub(crate) use edge_set_intersector::EdgeSetIntersector;
-pub(crate) use rstar_edge_set_intersector::RstarEdgeSetIntersector;
+pub use prepared_geometry::PreparedGeometry;
+pub(crate) use rstar_edge_set_intersector::RStarEdgeSetIntersector;
+pub(crate) use segment::Segment;
 pub(crate) use segment_intersector::SegmentIntersector;
 pub(crate) use simple_edge_set_intersector::SimpleEdgeSetIntersector;