From 2cb95a06e66bbd6888bab26639b07d438f7ebc00 Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Wed, 26 Feb 2025 15:04:46 -0800 Subject: [PATCH 1/2] Fix release warnings - these vars only used in debug --- .../algorithm/relate/geomgraph/edge_end_bundle_star.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/geo/src/algorithm/relate/geomgraph/edge_end_bundle_star.rs b/geo/src/algorithm/relate/geomgraph/edge_end_bundle_star.rs index fdf33a9ba1..b53cfa03e0 100644 --- a/geo/src/algorithm/relate/geomgraph/edge_end_bundle_star.rs +++ b/geo/src/algorithm/relate/geomgraph/edge_end_bundle_star.rs @@ -83,7 +83,7 @@ impl LabeledEdgeEndBundleStar { debug!("edge_end_bundle_star: {:?}", self); } - fn propagate_side_labels(&mut self, geom_index: usize, geometry_graph: &GeometryGraph) { + fn propagate_side_labels(&mut self, geom_index: usize, _geometry_graph: &GeometryGraph) { let mut start_position = None; for edge_ends in self.edge_end_bundles_iter() { @@ -108,11 +108,11 @@ impl LabeledEdgeEndBundleStar { let left_position = label.position(geom_index, Direction::Left); let right_position = label.position(geom_index, Direction::Right); - if let Some(right_position) = right_position { + if let Some(_right_position) = right_position { #[cfg(debug_assertions)] - if right_position != current_position { + if _right_position != current_position { use crate::algorithm::Validation; - if geometry_graph.geometry().is_valid() { + if _geometry_graph.geometry().is_valid() { debug_assert!(false, "topology position conflict with coordinate — this can happen with invalid geometries. coordinate: {:?}, right_location: {:?}, current_location: {:?}", edge_ends.coordinate(), right_position, current_position); } else { warn!("topology position conflict with coordinate — this can happen with invalid geometries. coordinate: {:?}, right_location: {:?}, current_location: {:?}", edge_ends.coordinate(), right_position, current_position); From 2efea9263c9d2108ca46b51693297545efd63dff Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Wed, 26 Feb 2025 14:21:57 -0800 Subject: [PATCH 2/2] Speed up Relate for PreparedGeometry Most of the speedup comes from caching bounding_rect == Bench $ cargo bench --bench prepared_geometry -- --baseline=main Compiling geo v0.29.3 (/Users/mkirk/src/georust/geo/geo) Compiling jts-test-runner v0.1.0 (/Users/mkirk/src/georust/geo/jts-test-runner) Finished `bench` profile [optimized] target(s) in 4.65s Running benches/prepared_geometry.rs (target/release/deps/prepared_geometry-b11c776cecdacd91) Gnuplot not found, using plotters backend relate prepared polygons time: [32.210 ms 32.290 ms 32.383 ms] change: [-33.711% -33.548% -33.350%] (p = 0.00 < 0.05) Performance has improved. Found 4 outliers among 100 measurements (4.00%) 3 (3.00%) high mild 1 (1.00%) high severe Benchmarking relate unprepared polygons: Warming up for 3.0000 s Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 80.0s, or reduce sample count to 10. relate unprepared polygons time: [799.79 ms 801.23 ms 802.67 ms] change: [-1.5726% -1.2283% -0.9119%] (p = 0.00 < 0.05) Change within noise threshold. Found 2 outliers among 100 measurements (2.00%) 2 (2.00%) high mild We get a little more speedup by skipping GeomtryGraph construction if the bounding-check fails == Bench NOTE: this is vs. *main*, so the prepared_geometry bench was already -33.548% faster vs main. $ cargo bench --bench prepared_geometry -- --baseline=main Compiling geo v0.29.3 (/Users/mkirk/src/georust/geo/geo) Compiling jts-test-runner v0.1.0 (/Users/mkirk/src/georust/geo/jts-test-runner) Finished `bench` profile [optimized] target(s) in 6.35s Running benches/prepared_geometry.rs (target/release/deps/prepared_geometry-b11c776cecdacd91) Gnuplot not found, using plotters backend relate prepared polygons time: [25.355 ms 25.456 ms 25.562 ms] change: [-47.861% -47.612% -47.365%] (p = 0.00 < 0.05) Performance has improved. Found 3 outliers among 100 measurements (3.00%) 3 (3.00%) high mild Benchmarking relate unprepared polygons: Warming up for 3.0000 s Warning: Unable to complete 100 samples in 5.0s. You may wish to increase target time to 82.0s, or reduce sample count to 10. relate unprepared polygons time: [761.98 ms 764.63 ms 767.56 ms] change: [-6.1618% -5.7402% -5.3188%] (p = 0.00 < 0.05) Performance has improved. Found 4 outliers among 100 measurements (4.00%) 3 (3.00%) high mild 1 (1.00%) high severe --- geo/CHANGES.md | 4 + .../geomgraph/index/prepared_geometry.rs | 53 +++++++- .../relate/geomgraph/intersection_matrix.rs | 11 +- geo/src/algorithm/relate/mod.rs | 21 ++-- geo/src/algorithm/relate/relate_operation.rs | 116 +++++++++--------- 5 files changed, 128 insertions(+), 77 deletions(-) diff --git a/geo/CHANGES.md b/geo/CHANGES.md index f21a61b0cd..0d66d51658 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -79,6 +79,10 @@ - - Bump `geo` MSRV to 1.80 and update CI - +- BREAKING: Speed up `Relate` for `PreparedGeometry` - this did require + changing some trait constraints, but they are unlikely to affect you in + practice unless you have your own Relate implementation. + - ## 0.29.3 - 2024.12.03 diff --git a/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs b/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs index 368ae1a04e..701c70118f 100644 --- a/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs +++ b/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs @@ -1,12 +1,13 @@ use super::Segment; use crate::geometry::*; use crate::relate::geomgraph::{GeometryGraph, RobustLineIntersector}; -use crate::GeometryCow; +use crate::{BoundingRect, GeometryCow, HasDimensions}; use crate::{GeoFloat, Relate}; use std::cell::RefCell; use std::rc::Rc; +use crate::dimensions::Dimensions; use rstar::{RTree, RTreeNum}; /// A `PreparedGeometry` can be more efficient than a plain Geometry when performing @@ -27,17 +28,19 @@ use rstar::{RTree, RTreeNum}; /// /// ``` pub struct PreparedGeometry<'a, F: GeoFloat + RTreeNum = f64> { - geometry_graph: GeometryGraph<'a, F>, + pub(crate) geometry_graph: GeometryGraph<'a, F>, + pub(crate) bounding_rect: Option>, } mod conversions { use crate::geometry_cow::GeometryCow; use crate::relate::geomgraph::{GeometryGraph, RobustLineIntersector}; - use crate::{GeoFloat, PreparedGeometry}; + use crate::{BoundingRect, GeoFloat, PreparedGeometry}; use geo_types::{ Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle, }; + use rstar::Envelope; use std::rc::Rc; impl From> for PreparedGeometry<'_, F> { @@ -156,13 +159,29 @@ mod conversions { impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { fn from(geometry: GeometryCow<'a, F>) -> Self { let mut geometry_graph = GeometryGraph::new(0, geometry); - geometry_graph.set_tree(Rc::new(geometry_graph.build_tree())); + let r_tree = geometry_graph.build_tree(); + + let envelope = r_tree.root().envelope(); + + // geo and rstar have different conventions for how to represet an empty bounding boxes + let bounding_rect: Option> = if envelope == rstar::AABB::new_empty() { + None + } else { + Some(Rect::new(envelope.lower(), envelope.upper())) + }; + // They should be equal - but we can save the computation in the `--release` case + // by using the bounding_rect from the RTree + debug_assert_eq!(bounding_rect, geometry_graph.geometry().bounding_rect()); + + geometry_graph.set_tree(Rc::new(r_tree)); // TODO: don't pass in line intersector here - in theory we'll want pluggable line intersectors // and the type (Robust) shouldn't be hard coded here. geometry_graph.compute_self_nodes(Box::new(RobustLineIntersector::new())); - - Self { geometry_graph } + Self { + geometry_graph, + bounding_rect, + } } } } @@ -176,6 +195,28 @@ where } } +impl BoundingRect for PreparedGeometry<'_, F> { + type Output = Option>; + + fn bounding_rect(&self) -> Option> { + self.bounding_rect + } +} + +impl HasDimensions for PreparedGeometry<'_, F> { + fn is_empty(&self) -> bool { + self.geometry_graph.geometry().is_empty() + } + + fn dimensions(&self) -> Dimensions { + self.geometry_graph.geometry().dimensions() + } + + fn boundary_dimensions(&self) -> Dimensions { + self.geometry_graph.geometry().boundary_dimensions() + } +} + impl Relate for PreparedGeometry<'_, F> { /// Efficiently builds a [`GeometryGraph`] which can then be used for topological /// computations. diff --git a/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs b/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs index 8a515887be..f67b59fb35 100644 --- a/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs +++ b/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs @@ -1,4 +1,6 @@ -use crate::{coordinate_position::CoordPos, dimensions::Dimensions, GeoNum, GeometryCow}; +use crate::{ + coordinate_position::CoordPos, dimensions::Dimensions, GeoNum, GeometryCow, HasDimensions, +}; use crate::geometry_cow::GeometryCow::Point; use crate::relate::geomgraph::intersection_matrix::dimension_matcher::DimensionMatcher; @@ -125,12 +127,11 @@ impl IntersectionMatrix { /// If the Geometries are disjoint, we need to enter their dimension and boundary dimension in /// the `Outside` rows in the IM - pub(crate) fn compute_disjoint( + pub(crate) fn compute_disjoint( &mut self, - geometry_a: &GeometryCow, - geometry_b: &GeometryCow, + geometry_a: &(impl HasDimensions + ?Sized), + geometry_b: &(impl HasDimensions + ?Sized), ) { - use crate::algorithm::dimensions::HasDimensions; { let dimensions = geometry_a.dimensions(); if dimensions != Dimensions::Empty { diff --git a/geo/src/algorithm/relate/mod.rs b/geo/src/algorithm/relate/mod.rs index 10c7dc23c0..bb97bddff3 100644 --- a/geo/src/algorithm/relate/mod.rs +++ b/geo/src/algorithm/relate/mod.rs @@ -5,7 +5,7 @@ use relate_operation::RelateOperation; use crate::geometry::*; pub use crate::relate::geomgraph::index::PreparedGeometry; pub use crate::relate::geomgraph::GeometryGraph; -use crate::{GeoFloat, GeometryCow}; +use crate::{BoundingRect, GeoFloat, GeometryCow, HasDimensions}; mod edge_end_builder; mod geomgraph; @@ -54,13 +54,20 @@ mod relate_operation; /// ``` /// /// Note: `Relate` must not be called on geometries containing `NaN` coordinates. -pub trait Relate { - /// Construct a [`GeometryGraph`] - fn geometry_graph(&self, arg_index: usize) -> GeometryGraph; +pub trait Relate: BoundingRect + HasDimensions { + /// Returns a noded topology graph for the geometry. + /// + /// # Params + /// + /// `idx`: 0 or 1, designating A or B (respectively) in the role this geometry plays + /// in the relation. e.g. in `a.relate(b)` + fn geometry_graph(&self, idx: usize) -> GeometryGraph; - fn relate(&self, other: &impl Relate) -> IntersectionMatrix { - RelateOperation::new(self.geometry_graph(0), other.geometry_graph(1)) - .compute_intersection_matrix() + fn relate(&self, other: &impl Relate) -> IntersectionMatrix + where + Self: Sized, + { + RelateOperation::new(self, other).compute_intersection_matrix() } } diff --git a/geo/src/algorithm/relate/relate_operation.rs b/geo/src/algorithm/relate/relate_operation.rs index 061bb29199..d03ee8c24b 100644 --- a/geo/src/algorithm/relate/relate_operation.rs +++ b/geo/src/algorithm/relate/relate_operation.rs @@ -6,9 +6,10 @@ use crate::relate::geomgraph::{ CoordNode, CoordPos, Edge, EdgeEnd, EdgeEndBundleStar, GeometryGraph, LabeledEdgeEndBundleStar, RobustLineIntersector, }; -use crate::CoordinatePosition; use crate::{Coord, GeoFloat, GeometryCow}; +use crate::{CoordinatePosition, Relate}; +use geo_types::Rect; use std::cell::RefCell; use std::rc::Rc; @@ -21,12 +22,14 @@ use std::rc::Rc; /// This implementation relies heavily on the functionality of [`GeometryGraph`]. /// /// Based on [JTS's `RelateComputer` as of 1.18.1](https://github.com/locationtech/jts/blob/jts-1.18.1/modules/core/src/main/java/org/locationtech/jts/operation/relate/RelateComputer.java) -pub(crate) struct RelateOperation<'a, F> +pub(crate) struct RelateOperation<'a, F, BBOX1, BBOX2> where F: GeoFloat, + BBOX1: Into>>, + BBOX2: Into>>, { - graph_a: GeometryGraph<'a, F>, - graph_b: GeometryGraph<'a, F>, + geometry_a: &'a dyn Relate, + geometry_b: &'a dyn Relate, nodes: NodeMap, line_intersector: RobustLineIntersector, isolated_edges: Vec>>>, @@ -44,14 +47,19 @@ where } } -impl<'a, F> RelateOperation<'a, F> +impl<'a, F, BBOX1, BBOX2> RelateOperation<'a, F, BBOX1, BBOX2> where F: GeoFloat, + BBOX1: Into>>, + BBOX2: Into>>, { - pub(crate) fn new(graph_a: GeometryGraph<'a, F>, graph_b: GeometryGraph<'a, F>) -> Self { + pub(crate) fn new( + geometry_a: &'a impl Relate, + geometry_b: &'a impl Relate, + ) -> Self { Self { - graph_a, - graph_b, + geometry_a, + geometry_b, nodes: NodeMap::new(), isolated_edges: vec![], line_intersector: RobustLineIntersector::new(), @@ -61,58 +69,61 @@ where pub(crate) fn compute_intersection_matrix(&mut self) -> IntersectionMatrix { let mut intersection_matrix = IntersectionMatrix::empty_disjoint(); - use crate::BoundingRect; use crate::Intersects; match ( - self.graph_a.geometry().bounding_rect(), - self.graph_b.geometry().bounding_rect(), + self.geometry_a.bounding_rect().into(), + self.geometry_b.bounding_rect().into(), ) { (Some(bounding_rect_a), Some(bounding_rect_b)) if bounding_rect_a.intersects(&bounding_rect_b) => {} _ => { // since Geometries don't overlap, we can skip most of the work - intersection_matrix - .compute_disjoint(self.graph_a.geometry(), self.graph_b.geometry()); + intersection_matrix.compute_disjoint(self.geometry_a, self.geometry_b); return intersection_matrix; } } + let mut graph_a = self.geometry_a.geometry_graph(0); + let mut graph_b = self.geometry_b.geometry_graph(1); + // Since changes to topology are inspected at nodes, we must crate a node for each // intersection. - self.graph_a - .compute_self_nodes(Box::new(self.line_intersector.clone())); - self.graph_b - .compute_self_nodes(Box::new(self.line_intersector.clone())); + graph_a.compute_self_nodes(Box::new(self.line_intersector.clone())); + graph_b.compute_self_nodes(Box::new(self.line_intersector.clone())); // compute intersections between edges of the two input geometries - let segment_intersector = self - .graph_a - .compute_edge_intersections(&self.graph_b, Box::new(self.line_intersector.clone())); + let segment_intersector = + graph_a.compute_edge_intersections(&graph_b, Box::new(self.line_intersector.clone())); - self.compute_intersection_nodes(0); - self.compute_intersection_nodes(1); + self.compute_intersection_nodes(&graph_a, 0); + self.compute_intersection_nodes(&graph_b, 1); // Copy the labelling for the nodes in the parent Geometries. These override any labels // determined by intersections between the geometries. - self.copy_nodes_and_labels(0); - self.copy_nodes_and_labels(1); + self.copy_nodes_and_labels(&graph_a, 0); + self.copy_nodes_and_labels(&graph_b, 1); // complete the labelling for any nodes which only have a label for a single geometry - self.label_isolated_nodes(); + self.label_isolated_nodes(&graph_a, &graph_b); // If a proper intersection was found, we can set a lower bound on the IM. - self.compute_proper_intersection_im(&segment_intersector, &mut intersection_matrix); + Self::compute_proper_intersection_im( + &graph_a, + &graph_b, + &segment_intersector, + &mut intersection_matrix, + ); // Now process improper intersections // (eg where one or other of the geometries has a vertex at the intersection point) // We need to compute the edge graph at all nodes to determine the IM. let edge_end_builder = EdgeEndBuilder::new(); - let edge_ends_a: Vec<_> = edge_end_builder.compute_ends_for_edges(self.graph_a.edges()); + let edge_ends_a: Vec<_> = edge_end_builder.compute_ends_for_edges(graph_a.edges()); self.insert_edge_ends(edge_ends_a); - let edge_ends_b: Vec<_> = edge_end_builder.compute_ends_for_edges(self.graph_b.edges()); + let edge_ends_b: Vec<_> = edge_end_builder.compute_ends_for_edges(graph_b.edges()); self.insert_edge_ends(edge_ends_b); let mut nodes = NodeMap::new(); std::mem::swap(&mut self.nodes, &mut nodes); let labeled_node_edges = nodes .into_iter() - .map(|(node, edges)| (node, edges.into_labeled(&self.graph_a, &self.graph_b))) + .map(|(node, edges)| (node, edges.into_labeled(&graph_a, &graph_b))) .collect(); // Compute the labeling for "isolated" components @@ -125,8 +136,8 @@ where // We only need to check components contained in the input graphs, since, by definition, // isolated components will not have been replaced by new components formed by // intersections. - self.label_isolated_edges(0, 1); - self.label_isolated_edges(1, 0); + self.label_isolated_edges(&graph_a, &graph_b, 1); + self.label_isolated_edges(&graph_b, &graph_a, 0); debug!( "before update_intersection_matrix: {:?}", @@ -147,13 +158,14 @@ where } fn compute_proper_intersection_im( - &mut self, + graph_a: &GeometryGraph, + graph_b: &GeometryGraph, segment_intersector: &SegmentIntersector, intersection_matrix: &mut IntersectionMatrix, ) { // If a proper intersection is found, we can set a lower bound on the IM. - let dim_a = self.graph_a.geometry().dimensions(); - let dim_b = self.graph_b.geometry().dimensions(); + let dim_a = graph_a.geometry().dimensions(); + let dim_b = graph_b.geometry().dimensions(); let has_proper = segment_intersector.has_proper_intersection(); let has_proper_interior = segment_intersector.has_proper_interior_intersection(); @@ -231,13 +243,7 @@ where /// argIndex. (E.g. a node may be an intersection node with a computed label of BOUNDARY, but /// in the original arg Geometry it is actually in the interior due to the Boundary /// Determination Rule) - fn copy_nodes_and_labels(&mut self, geom_index: usize) { - let graph = if geom_index == 0 { - &self.graph_a - } else { - assert_eq!(geom_index, 1); - &self.graph_b - }; + fn copy_nodes_and_labels(&mut self, graph: &GeometryGraph, geom_index: usize) { for graph_node in graph.nodes_iter() { let new_node = self .nodes @@ -259,14 +265,7 @@ where /// labeled. /// /// Endpoint nodes will already be labeled from when they were inserted. - fn compute_intersection_nodes(&mut self, geom_index: usize) { - let graph = if geom_index == 0 { - &self.graph_a - } else { - assert_eq!(geom_index, 1); - &self.graph_b - }; - + fn compute_intersection_nodes(&mut self, graph: &GeometryGraph, geom_index: usize) { for edge in graph.edges() { let edge = edge.borrow(); @@ -317,13 +316,12 @@ where /// By definition, "isolated" edges are guaranteed not to touch the boundary of the target /// (since if they did, they would have caused an intersection to be computed and hence would /// not be isolated). - fn label_isolated_edges(&mut self, this_index: usize, target_index: usize) { - let (this_graph, target_graph) = if this_index == 0 { - (&self.graph_a, &self.graph_b) - } else { - (&self.graph_b, &self.graph_a) - }; - + fn label_isolated_edges( + &mut self, + this_graph: &GeometryGraph, + target_graph: &GeometryGraph, + target_index: usize, + ) { for edge in this_graph.edges() { let mut mut_edge = edge.borrow_mut(); if mut_edge.is_isolated() { @@ -356,9 +354,9 @@ where /// are not completely labeled by the initial process of adding nodes to the nodeList. To /// complete the labelling we need to check for nodes that lie in the interior of edges, and in /// the interior of areas. - fn label_isolated_nodes(&mut self) { - let geometry_a = self.graph_a.geometry(); - let geometry_b = self.graph_b.geometry(); + fn label_isolated_nodes(&mut self, graph_a: &GeometryGraph, graph_b: &GeometryGraph) { + let geometry_a = graph_a.geometry(); + let geometry_b = graph_b.geometry(); for (node, _edges) in self.nodes.iter_mut() { let label = node.label(); // isolated nodes should always have at least one geometry in their label