From f64ae17dc8bb826f3e9a77465a0fbf0d6d3b4ccf Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Thu, 14 Jul 2022 21:26:15 +0100 Subject: [PATCH 1/2] existing bench --- geo/Cargo.toml | 4 ++++ geo/benches/prepared_geometry.rs | 31 +++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+) create mode 100644 geo/benches/prepared_geometry.rs diff --git a/geo/Cargo.toml b/geo/Cargo.toml index 0d2f58a7f..082db4af6 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -83,6 +83,10 @@ harness = false name = "euclidean_distance" harness = false +[[bench]] +name = "prepared_geometry" +harness = false + [[bench]] name = "rotate" harness = false diff --git a/geo/benches/prepared_geometry.rs b/geo/benches/prepared_geometry.rs new file mode 100644 index 000000000..c1c22b565 --- /dev/null +++ b/geo/benches/prepared_geometry.rs @@ -0,0 +1,31 @@ +use criterion::{criterion_group, criterion_main, Criterion}; +use geo::algorithm::Relate; +use geo::geometry::MultiPolygon; + +fn criterion_benchmark(c: &mut Criterion) { + 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(); + + bencher.iter(|| { + let mut intersects = 0; + let mut non_intersects = 0; + + 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); + }); + }); +} + +criterion_group!(benches, criterion_benchmark); +criterion_main!(benches); From 9452543f3764167d641ccd603cb4a29b214f81fe Mon Sep 17 00:00:00 2001 From: Michael Kirk Date: Thu, 14 Jul 2022 22:13:33 +0100 Subject: [PATCH 2/2] 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 | 2 + geo/benches/prepared_geometry.rs | 38 ++- geo/src/algorithm/relate/geomgraph/edge.rs | 10 +- .../relate/geomgraph/edge_intersection.rs | 2 +- .../relate/geomgraph/geometry_graph.rs | 136 +++++++---- .../geomgraph/index/edge_set_intersector.rs | 14 +- .../algorithm/relate/geomgraph/index/mod.rs | 6 +- .../geomgraph/index/prepared_geometry.rs | 227 ++++++++++++++++++ .../index/rstar_edge_set_intersector.rs | 123 +++------- .../relate/geomgraph/index/segment.rs | 33 +++ .../index/simple_edge_set_intersector.rs | 24 +- .../relate/geomgraph/intersection_matrix.rs | 52 +++- geo/src/algorithm/relate/geomgraph/label.rs | 6 +- geo/src/algorithm/relate/geomgraph/mod.rs | 2 +- geo/src/algorithm/relate/geomgraph/node.rs | 6 +- .../algorithm/relate/geomgraph/node_map.rs | 3 +- .../relate/geomgraph/planar_graph.rs | 33 +++ .../relate/geomgraph/topology_position.rs | 2 +- geo/src/algorithm/relate/mod.rs | 77 ++---- geo/src/algorithm/relate/relate_operation.rs | 78 +----- geo/src/geometry_cow.rs | 81 ++++++- geo/src/lib.rs | 1 + 22 files changed, 681 insertions(+), 275 deletions(-) create mode 100644 geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs create mode 100644 geo/src/algorithm/relate/geomgraph/index/segment.rs diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 47f60803b..ce1d36e19 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -10,6 +10,8 @@ * * Fix `AffineTransform::compose` ordering to be conventional - such that the argument is applied *after* self. * +* Add `PreparedGeometry` to speed up repeated `Relate` operations. + * ## 0.28.0 diff --git a/geo/benches/prepared_geometry.rs b/geo/benches/prepared_geometry.rs index c1c22b565..96253aa6b 100644 --- a/geo/benches/prepared_geometry.rs +++ b/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::>(); + + let zone_polygons = zone_polygons + .iter() + .map(PreparedGeometry::from) + .collect::>(); + + 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; diff --git a/geo/src/algorithm/relate/geomgraph/edge.rs b/geo/src/algorithm/relate/geomgraph/edge.rs index f21561d4c..c4d372ac8 100644 --- a/geo/src/algorithm/relate/geomgraph/edge.rs +++ b/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 { /// `coordinates` of the line geometry coords: Vec>, @@ -48,6 +48,13 @@ impl Edge { &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] { &self.coords } @@ -55,6 +62,7 @@ impl Edge { pub fn is_isolated(&self) -> bool { self.is_isolated } + pub fn mark_as_unisolated(&mut self) { self.is_isolated = false; } diff --git a/geo/src/algorithm/relate/geomgraph/edge_intersection.rs b/geo/src/algorithm/relate/geomgraph/edge_intersection.rs index b88908711..b8e4c8c04 100644 --- a/geo/src/algorithm/relate/geomgraph/edge_intersection.rs +++ b/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 { coord: Coord, segment_index: usize, diff --git a/geo/src/algorithm/relate/geomgraph/geometry_graph.rs b/geo/src/algorithm/relate/geomgraph/geometry_graph.rs index 714b38e66..6746dc8ac 100644 --- a/geo/src/algorithm/relate/geomgraph/geometry_graph.rs +++ b/geo/src/algorithm/relate/geomgraph/geometry_graph.rs @@ -1,6 +1,7 @@ use super::{ index::{ - EdgeSetIntersector, RstarEdgeSetIntersector, SegmentIntersector, SimpleEdgeSetIntersector, + EdgeSetIntersector, RStarEdgeSetIntersector, Segment, SegmentIntersector, + SimpleEdgeSetIntersector, }, CoordNode, CoordPos, Direction, Edge, Label, LineIntersector, PlanarGraph, TopologyPosition, }; @@ -8,17 +9,18 @@ use super::{ 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: @@ -26,13 +28,16 @@ use std::rc::Rc; /// - 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>>>, use_boundary_determination_rule: bool, + has_computed_self_nodes: bool, planar_graph: PlanarGraph, } @@ -44,44 +49,102 @@ impl GeometryGraph<'_, F> where F: GeoFloat, { - pub fn edges(&self) -> &[Rc>>] { + pub(crate) fn set_tree(&mut self, tree: Rc>>) { + self.tree = Some(tree); + } + + pub(crate) fn get_or_build_tree(&self) -> Rc>> { + self.tree + .clone() + .unwrap_or_else(|| Rc::new(self.build_tree())) + } + + pub(crate) fn build_tree(&self) -> RTree> { + let segments: Vec> = 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>>] { self.planar_graph.edges() } - pub fn insert_edge(&mut self, edge: Edge) { + pub(crate) fn insert_edge(&mut self, edge: Edge) { self.planar_graph.insert_edge(edge) } - pub fn is_boundary_node(&self, coord: Coord) -> bool { + pub(crate) fn is_boundary_node(&self, coord: Coord) -> bool { self.planar_graph.is_boundary_node(self.arg_index, coord) } - pub fn add_node_with_coordinate(&mut self, coord: Coord) -> &mut CoordNode { + pub(crate) fn add_node_with_coordinate(&mut self, coord: Coord) -> &mut CoordNode { self.planar_graph.add_node_with_coordinate(coord) } - pub fn nodes_iter(&self) -> impl Iterator> { + pub(crate) fn nodes_iter(&self) -> impl Iterator> { 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) -> 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 { - self.parent_geometry + pub(crate) fn geometry(&self) -> &GeometryCow { + &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> { - // 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> { self.planar_graph.boundary_nodes(self.arg_index) } - pub fn add_geometry(&mut self, geometry: &GeometryCow) { + pub(crate) fn add_geometry(&mut self, geometry: &GeometryCow) { 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>, - ) -> SegmentIntersector { - let mut segment_intersector = SegmentIntersector::new(line_intersector, true); + pub(crate) fn compute_self_nodes(&mut self, line_intersector: Box>) { + 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, line_intersector: Box>, @@ -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, ); diff --git a/geo/src/algorithm/relate/geomgraph/index/edge_set_intersector.rs b/geo/src/algorithm/relate/geomgraph/index/edge_set_intersector.rs index 709eb2bf1..a1f81af66 100644 --- a/geo/src/algorithm/relate/geomgraph/index/edge_set_intersector.rs +++ b/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 { /// `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>>], + &self, + graph: &GeometryGraph, check_for_self_intersecting_edges: bool, segment_intersector: &mut SegmentIntersector, ); /// Compute all intersections between two sets of edges, recording those intersections on /// the intersecting edges. - fn compute_intersections_between_sets( - &mut self, - edges0: &[Rc>>], - edges1: &[Rc>>], + fn compute_intersections_between_sets<'a>( + &self, + graph_0: &GeometryGraph<'a, F>, + graph_1: &GeometryGraph<'a, F>, segment_intersector: &mut SegmentIntersector, ); } diff --git a/geo/src/algorithm/relate/geomgraph/index/mod.rs b/geo/src/algorithm/relate/geomgraph/index/mod.rs index ae48b3db9..5695878e2 100644 --- a/geo/src/algorithm/relate/geomgraph/index/mod.rs +++ b/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; diff --git a/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs b/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs new file mode 100644 index 000000000..2fbc199e2 --- /dev/null +++ b/geo/src/algorithm/relate/geomgraph/index/prepared_geometry.rs @@ -0,0 +1,227 @@ +use super::Segment; +use crate::geometry::*; +use crate::relate::geomgraph::{GeometryGraph, RobustLineIntersector}; +use crate::GeometryCow; +use crate::{GeoFloat, Relate}; + +use std::cell::RefCell; +use std::rc::Rc; + +use rstar::{RTree, RTreeNum}; + +/// A `PreparedGeometry` can be more efficient than a plain Geometry when performing +/// multiple topological comparisons against the `PreparedGeometry`. +/// +/// ``` +/// use geo::{Relate, PreparedGeometry, wkt}; +/// +/// let polygon = wkt! { POLYGON((2.0 2.0,2.0 6.0,4.0 6.0)) }; +/// let touching_line = wkt! { LINESTRING(0.0 0.0,2.0 2.0) }; +/// let intersecting_line = wkt! { LINESTRING(0.0 0.0,3.0 3.0) }; +/// let contained_line = wkt! { LINESTRING(2.0 2.0,3.0 5.0) }; +/// +/// let prepared_polygon = PreparedGeometry::from(polygon); +/// assert!(prepared_polygon.relate(&touching_line).is_touches()); +/// assert!(prepared_polygon.relate(&intersecting_line).is_intersects()); +/// assert!(prepared_polygon.relate(&contained_line).is_contains()); +/// +/// ``` +pub struct PreparedGeometry<'a, F: GeoFloat + RTreeNum = f64> { + geometry_graph: GeometryGraph<'a, F>, +} + +mod conversions { + use crate::geometry_cow::GeometryCow; + use crate::relate::geomgraph::{GeometryGraph, RobustLineIntersector}; + use crate::{GeoFloat, PreparedGeometry}; + use geo_types::{ + Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, + Point, Polygon, Rect, Triangle, + }; + use std::rc::Rc; + + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(point: Point) -> Self { + PreparedGeometry::from(GeometryCow::from(point)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(line: Line) -> Self { + PreparedGeometry::from(GeometryCow::from(line)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(line_string: LineString) -> Self { + PreparedGeometry::from(GeometryCow::from(line_string)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(polygon: Polygon) -> Self { + PreparedGeometry::from(GeometryCow::from(polygon)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(multi_point: MultiPoint) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_point)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(multi_line_string: MultiLineString) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_line_string)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(multi_polygon: MultiPolygon) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_polygon)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(rect: Rect) -> Self { + PreparedGeometry::from(GeometryCow::from(rect)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(triangle: Triangle) -> Self { + PreparedGeometry::from(GeometryCow::from(triangle)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(geometry_collection: GeometryCollection) -> Self { + PreparedGeometry::from(GeometryCow::from(geometry_collection)) + } + } + impl<'a, F: GeoFloat> From> for PreparedGeometry<'a, F> { + fn from(geometry: Geometry) -> Self { + PreparedGeometry::from(GeometryCow::from(geometry)) + } + } + + impl<'a, F: GeoFloat> From<&'a Point> for PreparedGeometry<'a, F> { + fn from(point: &'a Point) -> Self { + PreparedGeometry::from(GeometryCow::from(point)) + } + } + impl<'a, F: GeoFloat> From<&'a Line> for PreparedGeometry<'a, F> { + fn from(line: &'a Line) -> Self { + PreparedGeometry::from(GeometryCow::from(line)) + } + } + impl<'a, F: GeoFloat> From<&'a LineString> for PreparedGeometry<'a, F> { + fn from(line_string: &'a LineString) -> Self { + PreparedGeometry::from(GeometryCow::from(line_string)) + } + } + impl<'a, F: GeoFloat> From<&'a Polygon> for PreparedGeometry<'a, F> { + fn from(polygon: &'a Polygon) -> Self { + PreparedGeometry::from(GeometryCow::from(polygon)) + } + } + impl<'a, F: GeoFloat> From<&'a MultiPoint> for PreparedGeometry<'a, F> { + fn from(multi_point: &'a MultiPoint) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_point)) + } + } + impl<'a, F: GeoFloat> From<&'a MultiLineString> for PreparedGeometry<'a, F> { + fn from(multi_line_string: &'a MultiLineString) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_line_string)) + } + } + impl<'a, F: GeoFloat> From<&'a MultiPolygon> for PreparedGeometry<'a, F> { + fn from(multi_polygon: &'a MultiPolygon) -> Self { + PreparedGeometry::from(GeometryCow::from(multi_polygon)) + } + } + impl<'a, F: GeoFloat> From<&'a GeometryCollection> for PreparedGeometry<'a, F> { + fn from(geometry_collection: &'a GeometryCollection) -> Self { + PreparedGeometry::from(GeometryCow::from(geometry_collection)) + } + } + impl<'a, F: GeoFloat> From<&'a Rect> for PreparedGeometry<'a, F> { + fn from(rect: &'a Rect) -> Self { + PreparedGeometry::from(GeometryCow::from(rect)) + } + } + impl<'a, F: GeoFloat> From<&'a Triangle> for PreparedGeometry<'a, F> { + fn from(triangle: &'a Triangle) -> Self { + PreparedGeometry::from(GeometryCow::from(triangle)) + } + } + + impl<'a, F: GeoFloat> From<&'a Geometry> for PreparedGeometry<'a, F> { + fn from(geometry: &'a Geometry) -> Self { + PreparedGeometry::from(GeometryCow::from(geometry)) + } + } + + 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())); + + // 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 } + } + } +} + +impl<'a, F> PreparedGeometry<'a, F> +where + F: GeoFloat + RTreeNum, +{ + pub(crate) fn geometry(&self) -> &GeometryCow { + self.geometry_graph.geometry() + } +} + +impl Relate for PreparedGeometry<'_, F> { + /// Efficiently builds a [`GeometryGraph`] which can then be used for topological + /// computations. + fn geometry_graph(&self, arg_index: usize) -> GeometryGraph { + self.geometry_graph.clone_for_arg_index(arg_index) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::algorithm::Relate; + use crate::polygon; + + #[test] + fn relate() { + let p1 = polygon![(x: 0.0, y: 0.0), (x: 2.0, y: 0.0), (x: 1.0, y: 1.0)]; + let p2 = polygon![(x: 0.5, y: 0.0), (x: 2.0, y: 0.0), (x: 1.0, y: 1.0)]; + let prepared_1 = PreparedGeometry::from(&p1); + let prepared_2 = PreparedGeometry::from(&p2); + assert!(prepared_1.relate(&prepared_2).is_contains()); + assert!(prepared_2.relate(&prepared_1).is_within()); + } + + #[test] + fn prepared_with_unprepared() { + let p1 = polygon![(x: 0.0, y: 0.0), (x: 2.0, y: 0.0), (x: 1.0, y: 1.0)]; + let p2 = polygon![(x: 0.5, y: 0.0), (x: 2.0, y: 0.0), (x: 1.0, y: 1.0)]; + let prepared_1 = PreparedGeometry::from(&p1); + assert!(prepared_1.relate(&p2).is_contains()); + assert!(p2.relate(&prepared_1).is_within()); + } + + #[test] + fn swap_arg_index() { + let poly = polygon![(x: 0.0, y: 0.0), (x: 2.0, y: 0.0), (x: 1.0, y: 1.0)]; + let prepared_geom = PreparedGeometry::from(&poly); + + let poly_cow = GeometryCow::from(&poly); + + let cached_graph = prepared_geom.geometry_graph(0); + let fresh_graph = GeometryGraph::new(0, poly_cow.clone()); + cached_graph.assert_eq_graph(&fresh_graph); + + let cached_graph = prepared_geom.geometry_graph(1); + let fresh_graph = GeometryGraph::new(1, poly_cow); + cached_graph.assert_eq_graph(&fresh_graph); + } +} 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 63c7b2b69..5c0190d33 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 @@ -1,105 +1,62 @@ -use super::super::Edge; -use super::{EdgeSetIntersector, SegmentIntersector}; +use super::super::{Edge, GeometryGraph}; +use super::{EdgeSetIntersector, Segment, SegmentIntersector}; use crate::GeoFloat; use std::cell::RefCell; use std::rc::Rc; -use rstar::RTree; +use rstar::{RTree, RTreeNum}; -pub(crate) struct RstarEdgeSetIntersector; +pub(crate) struct RStarEdgeSetIntersector; -impl RstarEdgeSetIntersector { - pub fn new() -> Self { - RstarEdgeSetIntersector - } -} - -struct Segment<'a, F: GeoFloat + rstar::RTreeNum> { - i: usize, - edge: &'a RefCell>, - envelope: rstar::AABB>, -} - -impl<'a, F> Segment<'a, F> -where - F: GeoFloat + rstar::RTreeNum, -{ - fn new(i: usize, edge: &'a RefCell>) -> Self { - use 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<'a, F> rstar::RTreeObject for Segment<'a, F> -where - F: GeoFloat + rstar::RTreeNum, -{ - type Envelope = rstar::AABB>; - - fn envelope(&self) -> Self::Envelope { - self.envelope - } -} - -impl EdgeSetIntersector for RstarEdgeSetIntersector +impl EdgeSetIntersector for RStarEdgeSetIntersector where - F: GeoFloat + rstar::RTreeNum, + F: GeoFloat + RTreeNum, { fn compute_intersections_within_set( - &mut self, - edges: &[Rc>>], + &self, + graph: &GeometryGraph, check_for_self_intersecting_edges: bool, segment_intersector: &mut SegmentIntersector, ) { - let segments: Vec> = 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); + let edges = graph.edges(); + + let tree = graph.get_or_build_tree(); + for (segment_0, segment_1) in tree.intersection_candidates_with_other_tree(&tree) { + if check_for_self_intersecting_edges || segment_0.edge_idx != segment_1.edge_idx { + let edge_0 = &edges[segment_0.edge_idx]; + let edge_1 = &edges[segment_1.edge_idx]; + segment_intersector.add_intersections( + edge_0, + segment_0.segment_idx, + edge_1, + segment_1.segment_idx, + ); } } } - fn compute_intersections_between_sets( - &mut self, - edges0: &[Rc>>], - edges1: &[Rc>>], + fn compute_intersections_between_sets<'a>( + &self, + graph_0: &GeometryGraph<'a, F>, + graph_1: &GeometryGraph<'a, F>, segment_intersector: &mut SegmentIntersector, ) { - let segments0: Vec> = 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> = 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); + let edges_0 = graph_0.edges(); + let edges_1 = graph_1.edges(); + + let tree_0 = graph_0.get_or_build_tree(); + let tree_1 = graph_1.get_or_build_tree(); + + for (segment_0, segment_1) in tree_0.intersection_candidates_with_other_tree(&tree_1) { + let edge_0 = &edges_0[segment_0.edge_idx]; + let edge_1 = &edges_1[segment_1.edge_idx]; + segment_intersector.add_intersections( + edge_0, + segment_0.segment_idx, + edge_1, + segment_1.segment_idx, + ); } } } diff --git a/geo/src/algorithm/relate/geomgraph/index/segment.rs b/geo/src/algorithm/relate/geomgraph/index/segment.rs new file mode 100644 index 000000000..11c5c2543 --- /dev/null +++ b/geo/src/algorithm/relate/geomgraph/index/segment.rs @@ -0,0 +1,33 @@ +use crate::Coord; +use crate::GeoFloat; + +#[derive(Debug, Clone)] +pub(crate) struct Segment { + pub edge_idx: usize, + pub segment_idx: usize, + pub envelope: rstar::AABB>, +} + +impl Segment +where + F: GeoFloat + rstar::RTreeNum, +{ + pub fn new(edge_idx: usize, segment_idx: usize, p1: Coord, p2: Coord) -> Self { + Self { + edge_idx, + segment_idx, + envelope: rstar::AABB::from_corners(p1, p2), + } + } +} + +impl rstar::RTreeObject for Segment +where + F: GeoFloat + rstar::RTreeNum, +{ + type Envelope = rstar::AABB>; + + fn envelope(&self) -> Self::Envelope { + self.envelope + } +} diff --git a/geo/src/algorithm/relate/geomgraph/index/simple_edge_set_intersector.rs b/geo/src/algorithm/relate/geomgraph/index/simple_edge_set_intersector.rs index db9cc9845..a73593446 100644 --- a/geo/src/algorithm/relate/geomgraph/index/simple_edge_set_intersector.rs +++ b/geo/src/algorithm/relate/geomgraph/index/simple_edge_set_intersector.rs @@ -1,4 +1,4 @@ -use super::super::Edge; +use super::super::{Edge, GeometryGraph}; use super::{EdgeSetIntersector, SegmentIntersector}; use crate::GeoFloat; @@ -13,7 +13,7 @@ impl SimpleEdgeSetIntersector { } fn compute_intersects( - &mut self, + &self, edge0: &Rc>>, edge1: &Rc>>, segment_intersector: &mut SegmentIntersector, @@ -30,11 +30,12 @@ impl SimpleEdgeSetIntersector { impl EdgeSetIntersector for SimpleEdgeSetIntersector { fn compute_intersections_within_set( - &mut self, - edges: &[Rc>>], + &self, + graph: &GeometryGraph, check_for_self_intersecting_edges: bool, segment_intersector: &mut SegmentIntersector, ) { + let edges = graph.edges(); for edge0 in edges.iter() { for edge1 in edges.iter() { if check_for_self_intersecting_edges || edge0.as_ptr() != edge1.as_ptr() { @@ -44,14 +45,17 @@ impl EdgeSetIntersector for SimpleEdgeSetIntersector { } } - fn compute_intersections_between_sets( - &mut self, - edges0: &[Rc>>], - edges1: &[Rc>>], + fn compute_intersections_between_sets<'a>( + &self, + graph_0: &GeometryGraph<'a, F>, + graph_1: &GeometryGraph<'a, F>, segment_intersector: &mut SegmentIntersector, ) { - for edge0 in edges0 { - for edge1 in edges1 { + let edges_0 = graph_0.edges(); + let edges_1 = graph_1.edges(); + + for edge0 in edges_0 { + for edge1 in edges_1 { self.compute_intersects(edge0, edge1, segment_intersector); } } diff --git a/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs b/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs index 1cd9f30ac..c365edf80 100644 --- a/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs +++ b/geo/src/algorithm/relate/geomgraph/intersection_matrix.rs @@ -1,6 +1,7 @@ -use crate::{coordinate_position::CoordPos, dimensions::Dimensions}; +use crate::{coordinate_position::CoordPos, dimensions::Dimensions, GeoNum, GeometryCow}; use crate::geometry_cow::GeometryCow::Point; +use crate::relate::geomgraph::intersection_matrix::dimension_matcher::DimensionMatcher; use std::str::FromStr; /// Models a *Dimensionally Extended Nine-Intersection Model (DE-9IM)* matrix. @@ -104,10 +105,57 @@ impl std::fmt::Debug for IntersectionMatrix { } impl IntersectionMatrix { - pub fn empty() -> Self { + pub const fn empty() -> Self { IntersectionMatrix(LocationArray([LocationArray([Dimensions::Empty; 3]); 3])) } + pub(crate) const fn empty_disjoint() -> Self { + IntersectionMatrix(LocationArray([ + LocationArray([Dimensions::Empty, Dimensions::Empty, Dimensions::Empty]), + LocationArray([Dimensions::Empty, Dimensions::Empty, Dimensions::Empty]), + // since Geometries are finite and embedded in a 2-D space, + // the `(Outside, Outside)` element must always be 2-D + LocationArray([ + Dimensions::Empty, + Dimensions::Empty, + Dimensions::TwoDimensional, + ]), + ])) + } + + /// 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( + &mut self, + geometry_a: &GeometryCow, + geometry_b: &GeometryCow, + ) { + use crate::algorithm::dimensions::HasDimensions; + { + let dimensions = geometry_a.dimensions(); + if dimensions != Dimensions::Empty { + self.set(CoordPos::Inside, CoordPos::Outside, dimensions); + + let boundary_dimensions = geometry_a.boundary_dimensions(); + if boundary_dimensions != Dimensions::Empty { + self.set(CoordPos::OnBoundary, CoordPos::Outside, boundary_dimensions); + } + } + } + + { + let dimensions = geometry_b.dimensions(); + if dimensions != Dimensions::Empty { + self.set(CoordPos::Outside, CoordPos::Inside, dimensions); + + let boundary_dimensions = geometry_b.boundary_dimensions(); + if boundary_dimensions != Dimensions::Empty { + self.set(CoordPos::Outside, CoordPos::OnBoundary, boundary_dimensions); + } + } + } + } + /// Set `dimensions` of the cell specified by the positions. /// /// `position_a`: which position `dimensions` applies to within the first geometry diff --git a/geo/src/algorithm/relate/geomgraph/label.rs b/geo/src/algorithm/relate/geomgraph/label.rs index ee7ff64b7..85da82fed 100644 --- a/geo/src/algorithm/relate/geomgraph/label.rs +++ b/geo/src/algorithm/relate/geomgraph/label.rs @@ -13,7 +13,7 @@ use std::fmt; /// /// If the component has *no* incidence with one of the geometries, than the `Label`'s /// `TopologyPosition` for that geometry is called `empty`. -#[derive(Clone)] +#[derive(Clone, PartialEq)] pub(crate) struct Label { geometry_topologies: [TopologyPosition; 2], } @@ -29,6 +29,10 @@ impl fmt::Debug for Label { } impl Label { + pub fn swap_args(&mut self) { + self.geometry_topologies.swap(0, 1) + } + /// Construct an empty `Label` for relating a 1-D line or 0-D point to both geometries. pub fn empty_line_or_point() -> Label { Label { diff --git a/geo/src/algorithm/relate/geomgraph/mod.rs b/geo/src/algorithm/relate/geomgraph/mod.rs index 6405d6dc0..62af0bee0 100644 --- a/geo/src/algorithm/relate/geomgraph/mod.rs +++ b/geo/src/algorithm/relate/geomgraph/mod.rs @@ -8,7 +8,7 @@ pub(crate) use edge_end::{EdgeEnd, EdgeEndKey}; pub(crate) use edge_end_bundle::{EdgeEndBundle, LabeledEdgeEndBundle}; pub(crate) use edge_end_bundle_star::{EdgeEndBundleStar, LabeledEdgeEndBundleStar}; pub(crate) use edge_intersection::EdgeIntersection; -pub(crate) use geometry_graph::GeometryGraph; +pub use geometry_graph::GeometryGraph; pub(crate) use intersection_matrix::IntersectionMatrix; pub(crate) use label::Label; pub(crate) use line_intersector::{LineIntersection, LineIntersector}; diff --git a/geo/src/algorithm/relate/geomgraph/node.rs b/geo/src/algorithm/relate/geomgraph/node.rs index 550fca210..ed77d524a 100644 --- a/geo/src/algorithm/relate/geomgraph/node.rs +++ b/geo/src/algorithm/relate/geomgraph/node.rs @@ -1,7 +1,7 @@ use super::{CoordPos, Dimensions, EdgeEnd, EdgeEndBundleStar, IntersectionMatrix, Label}; use crate::{Coord, GeoFloat}; -#[derive(Debug, Clone)] +#[derive(Debug, Clone, PartialEq)] pub(crate) struct CoordNode where F: GeoFloat, @@ -11,6 +11,10 @@ where } impl CoordNode { + pub fn swap_label_args(&mut self) { + self.label.swap_args() + } + pub(crate) fn label(&self) -> &Label { &self.label } diff --git a/geo/src/algorithm/relate/geomgraph/node_map.rs b/geo/src/algorithm/relate/geomgraph/node_map.rs index 3ae78cbf9..b77e67565 100644 --- a/geo/src/algorithm/relate/geomgraph/node_map.rs +++ b/geo/src/algorithm/relate/geomgraph/node_map.rs @@ -6,6 +6,7 @@ use std::fmt; use std::marker::PhantomData; /// A map of nodes, indexed by the coordinate of the node +#[derive(Clone, PartialEq)] pub(crate) struct NodeMap where F: GeoFloat, @@ -16,7 +17,7 @@ where } /// Creates the node stored in `NodeMap` -pub(crate) trait NodeFactory { +pub(crate) trait NodeFactory: PartialEq { type Node; fn create_node(coordinate: Coord) -> Self::Node; } diff --git a/geo/src/algorithm/relate/geomgraph/planar_graph.rs b/geo/src/algorithm/relate/geomgraph/planar_graph.rs index 0acf8e405..454e410f1 100644 --- a/geo/src/algorithm/relate/geomgraph/planar_graph.rs +++ b/geo/src/algorithm/relate/geomgraph/planar_graph.rs @@ -7,6 +7,7 @@ use crate::{Coord, GeoFloat}; use std::cell::RefCell; use std::rc::Rc; +#[derive(Clone, PartialEq)] pub(crate) struct PlanarGraphNode; /// The basic node constructor does not allow for incident edges @@ -20,12 +21,44 @@ where } } +#[derive(Clone, PartialEq)] pub(crate) struct PlanarGraph { pub(crate) nodes: NodeMap, edges: Vec>>>, } impl PlanarGraph { + pub fn clone_for_arg_index(&self, from_arg_index: usize, to_arg_index: usize) -> Self { + let mut graph = Self { + nodes: self.nodes.clone(), + // deep copy edges + edges: self + .edges + .iter() + .map(|e| Rc::new(RefCell::new(e.borrow().clone()))) + .collect(), + }; + assert_eq!(from_arg_index, 0); + if from_arg_index != to_arg_index { + graph.swap_labels(); + } + graph + } + + fn swap_labels(&mut self) { + for node in self.nodes.iter_mut() { + node.swap_label_args(); + } + for edge in &mut self.edges { + edge.borrow_mut().swap_label_args(); + } + } + + pub fn assert_eq_graph(&self, other: &Self) { + assert_eq!(self.nodes, other.nodes); + assert_eq!(self.edges, other.edges); + } + pub fn edges(&self) -> &[Rc>>] { &self.edges } diff --git a/geo/src/algorithm/relate/geomgraph/topology_position.rs b/geo/src/algorithm/relate/geomgraph/topology_position.rs index 11cda7f57..17435ff56 100644 --- a/geo/src/algorithm/relate/geomgraph/topology_position.rs +++ b/geo/src/algorithm/relate/geomgraph/topology_position.rs @@ -14,7 +14,7 @@ use std::fmt; /// topological relationship attribute for the [`On`](Direction::On) position. /// /// See [`CoordPos`] for the possible values. -#[derive(Copy, Clone)] +#[derive(Copy, Clone, PartialEq)] pub(crate) enum TopologyPosition { Area { on: Option, diff --git a/geo/src/algorithm/relate/mod.rs b/geo/src/algorithm/relate/mod.rs index f6ba3f851..1959d254d 100644 --- a/geo/src/algorithm/relate/mod.rs +++ b/geo/src/algorithm/relate/mod.rs @@ -1,7 +1,10 @@ pub(crate) use edge_end_builder::EdgeEndBuilder; pub use geomgraph::intersection_matrix::IntersectionMatrix; +use relate_operation::RelateOperation; use crate::geometry::*; +pub use crate::relate::geomgraph::index::PreparedGeometry; +pub use crate::relate::geomgraph::GeometryGraph; use crate::{GeoFloat, GeometryCow}; mod edge_end_builder; @@ -51,66 +54,38 @@ mod relate_operation; /// ``` /// /// Note: `Relate` must not be called on geometries containing `NaN` coordinates. -pub trait Relate { - fn relate(&self, other: &T) -> IntersectionMatrix; -} +pub trait Relate { + /// Construct a [`GeometryGraph`] + fn geometry_graph(&self, arg_index: usize) -> GeometryGraph; -impl Relate> for GeometryCow<'_, F> { - fn relate(&self, other: &GeometryCow) -> IntersectionMatrix { - let mut relate_computer = relate_operation::RelateOperation::new(self, other); - relate_computer.compute_intersection_matrix() + fn relate(&self, other: &impl Relate) -> IntersectionMatrix { + RelateOperation::new(self.geometry_graph(0), other.geometry_graph(1)) + .compute_intersection_matrix() } } macro_rules! relate_impl { - ($k:ty, $t:ty) => { - relate_impl![($k, $t),]; - }; - ($(($k:ty, $t:ty),)*) => { + ($($t:ty ,)*) => { $( - impl Relate for $k { - fn relate(&self, other: &$t) -> IntersectionMatrix { - GeometryCow::from(self).relate(&GeometryCow::from(other)) + impl Relate for $t { + fn geometry_graph(&self, arg_index: usize) -> GeometryGraph { + GeometryGraph::new(arg_index, GeometryCow::from(self)) } } )* }; } -/// Call the given macro with every pair of inputs -/// -/// # Examples -/// -/// ```ignore -/// cartesian_pairs!(foo, [Bar, Baz, Qux]); -/// ``` -/// Is akin to calling: -/// ```ignore -/// foo![(Bar, Bar), (Bar, Baz), (Bar, Qux), (Baz, Bar), (Baz, Baz), (Baz, Qux), (Qux, Bar), (Qux, Baz), (Qux, Qux)]; -/// ``` -macro_rules! cartesian_pairs { - ($macro_name:ident, [$($a:ty),*]) => { - cartesian_pairs_helper! { [] [$($a,)*] [$($a,)*] [$($a,)*] $macro_name} - }; -} - -macro_rules! cartesian_pairs_helper { - // popped all a's - we're done. Use the accumulated output as the input to relate macro. - ([$($out_pairs:tt)*] [] [$($b:ty,)*] $init_b:tt $macro_name:ident) => { - $macro_name!{$($out_pairs)*} - }; - // finished one loop of b, pop next a and reset b - ($out_pairs:tt [$a_car:ty, $($a_cdr:ty,)*] [] $init_b:tt $macro_name:ident) => { - cartesian_pairs_helper!{$out_pairs [$($a_cdr,)*] $init_b $init_b $macro_name} - }; - // pop b through all of b with head of a - ([$($out_pairs:tt)*] [$a_car:ty, $($a_cdr:ty,)*] [$b_car:ty, $($b_cdr:ty,)*] $init_b:tt $macro_name:ident) => { - cartesian_pairs_helper!{[$($out_pairs)* ($a_car, $b_car),] [$a_car, $($a_cdr,)*] [$($b_cdr,)*] $init_b $macro_name} - }; -} - -// Implement Relate for every combination of Geometry. Alternatively we could do something like -// `impl Relate> for Into { }` -// but I don't know that we want to make GeometryCow public (yet?). -cartesian_pairs!(relate_impl, [Point, Line, LineString, Polygon, MultiPoint, MultiLineString, MultiPolygon, Rect, Triangle, GeometryCollection]); -relate_impl!(Geometry, Geometry); +relate_impl![ + Point, + Line, + LineString, + Polygon, + MultiPoint, + MultiLineString, + MultiPolygon, + Rect, + Triangle, + GeometryCollection, + Geometry, +]; diff --git a/geo/src/algorithm/relate/relate_operation.rs b/geo/src/algorithm/relate/relate_operation.rs index 4d8676e41..061bb2919 100644 --- a/geo/src/algorithm/relate/relate_operation.rs +++ b/geo/src/algorithm/relate/relate_operation.rs @@ -32,6 +32,7 @@ where isolated_edges: Vec>>>, } +#[derive(PartialEq)] pub(crate) struct RelateNodeFactory; impl NodeFactory for RelateNodeFactory where @@ -47,13 +48,10 @@ impl<'a, F> RelateOperation<'a, F> where F: GeoFloat, { - pub(crate) fn new( - geom_a: &'a GeometryCow<'a, F>, - geom_b: &'a GeometryCow<'a, F>, - ) -> RelateOperation<'a, F> { + pub(crate) fn new(graph_a: GeometryGraph<'a, F>, graph_b: GeometryGraph<'a, F>) -> Self { Self { - graph_a: GeometryGraph::new(0, geom_a), - graph_b: GeometryGraph::new(1, geom_b), + graph_a, + graph_b, nodes: NodeMap::new(), isolated_edges: vec![], line_intersector: RobustLineIntersector::new(), @@ -61,14 +59,7 @@ where } pub(crate) fn compute_intersection_matrix(&mut self) -> IntersectionMatrix { - let mut intersection_matrix = IntersectionMatrix::empty(); - // since Geometries are finite and embedded in a 2-D space, - // the `(Outside, Outside)` element must always be 2-D - intersection_matrix.set( - CoordPos::Outside, - CoordPos::Outside, - Dimensions::TwoDimensional, - ); + let mut intersection_matrix = IntersectionMatrix::empty_disjoint(); use crate::BoundingRect; use crate::Intersects; @@ -80,7 +71,8 @@ where if bounding_rect_a.intersects(&bounding_rect_b) => {} _ => { // since Geometries don't overlap, we can skip most of the work - self.compute_disjoint_intersection_matrix(&mut intersection_matrix); + intersection_matrix + .compute_disjoint(self.graph_a.geometry(), self.graph_b.geometry()); return intersection_matrix; } } @@ -293,44 +285,6 @@ where } } - /// If the Geometries are disjoint, we need to enter their dimension and boundary dimension in - /// the `Outside` rows in the IM - fn compute_disjoint_intersection_matrix(&self, intersection_matrix: &mut IntersectionMatrix) { - { - let geometry_a = self.graph_a.geometry(); - let dimensions = geometry_a.dimensions(); - if dimensions != Dimensions::Empty { - intersection_matrix.set(CoordPos::Inside, CoordPos::Outside, dimensions); - - let boundary_dimensions = geometry_a.boundary_dimensions(); - if boundary_dimensions != Dimensions::Empty { - intersection_matrix.set( - CoordPos::OnBoundary, - CoordPos::Outside, - boundary_dimensions, - ); - } - } - } - - { - let geometry_b = self.graph_b.geometry(); - let dimensions = geometry_b.dimensions(); - if dimensions != Dimensions::Empty { - intersection_matrix.set(CoordPos::Outside, CoordPos::Inside, dimensions); - - let boundary_dimensions = geometry_b.boundary_dimensions(); - if boundary_dimensions != Dimensions::Empty { - intersection_matrix.set( - CoordPos::Outside, - CoordPos::OnBoundary, - boundary_dimensions, - ); - } - } - } - } - fn update_intersection_matrix( &self, labeled_node_edges: Vec<(CoordNode, LabeledEdgeEndBundleStar)>, @@ -457,12 +411,8 @@ mod test { ] .into(); - let gc1 = GeometryCow::from(&square_a); - let gc2 = GeometryCow::from(&square_b); - let mut relate_computer = RelateOperation::new(&gc1, &gc2); - let intersection_matrix = relate_computer.compute_intersection_matrix(); assert_eq!( - intersection_matrix, + square_a.relate(&square_b), IntersectionMatrix::from_str("FF2FF1212").unwrap() ); } @@ -487,12 +437,8 @@ mod test { ] .into(); - let gca = GeometryCow::from(&square_a); - let gcb = GeometryCow::from(&square_b); - let mut relate_computer = RelateOperation::new(&gca, &gcb); - let intersection_matrix = relate_computer.compute_intersection_matrix(); assert_eq!( - intersection_matrix, + square_a.relate(&square_b), IntersectionMatrix::from_str("212FF1FF2").unwrap() ); } @@ -517,12 +463,8 @@ mod test { ] .into(); - let gca = &GeometryCow::from(&square_a); - let gcb = &GeometryCow::from(&square_b); - let mut relate_computer = RelateOperation::new(gca, gcb); - let intersection_matrix = relate_computer.compute_intersection_matrix(); assert_eq!( - intersection_matrix, + square_a.relate(&square_b), IntersectionMatrix::from_str("212101212").unwrap() ); } diff --git a/geo/src/geometry_cow.rs b/geo/src/geometry_cow.rs index 5cd1b2b51..108c24ab4 100644 --- a/geo/src/geometry_cow.rs +++ b/geo/src/geometry_cow.rs @@ -10,7 +10,7 @@ use std::borrow::Cow; /// This is a way to "upgrade" an inner type to something like a `Geometry` without `moving` it. /// /// As an example, see the [`Relate`] trait which uses `GeometryCow`. -#[derive(PartialEq, Debug, Hash)] +#[derive(PartialEq, Debug, Hash, Clone)] pub(crate) enum GeometryCow<'a, T> where T: CoordNum, @@ -103,3 +103,82 @@ impl<'a, T: CoordNum> From<&'a Triangle> for GeometryCow<'a, T> { GeometryCow::Triangle(Cow::Borrowed(triangle)) } } + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(point: Point) -> Self { + GeometryCow::Point(Cow::Owned(point)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(line_string: LineString) -> Self { + GeometryCow::LineString(Cow::Owned(line_string)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(line: Line) -> Self { + GeometryCow::Line(Cow::Owned(line)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(polygon: Polygon) -> Self { + GeometryCow::Polygon(Cow::Owned(polygon)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(multi_point: MultiPoint) -> GeometryCow<'a, T> { + GeometryCow::MultiPoint(Cow::Owned(multi_point)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(multi_line_string: MultiLineString) -> Self { + GeometryCow::MultiLineString(Cow::Owned(multi_line_string)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(multi_polygon: MultiPolygon) -> Self { + GeometryCow::MultiPolygon(Cow::Owned(multi_polygon)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(geometry_collection: GeometryCollection) -> Self { + GeometryCow::GeometryCollection(Cow::Owned(geometry_collection)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(rect: Rect) -> Self { + GeometryCow::Rect(Cow::Owned(rect)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(triangle: Triangle) -> Self { + GeometryCow::Triangle(Cow::Owned(triangle)) + } +} + +impl<'a, T: CoordNum> From> for GeometryCow<'a, T> { + fn from(geometry: Geometry) -> Self { + match geometry { + Geometry::Point(point) => GeometryCow::from(point), + Geometry::Line(line) => GeometryCow::from(line), + Geometry::LineString(line_string) => GeometryCow::from(line_string), + Geometry::Polygon(polygon) => GeometryCow::from(polygon), + Geometry::MultiPoint(multi_point) => GeometryCow::from(multi_point), + Geometry::MultiLineString(multi_line_string) => GeometryCow::from(multi_line_string), + Geometry::MultiPolygon(multi_polygon) => GeometryCow::from(multi_polygon), + Geometry::GeometryCollection(geometry_collection) => { + GeometryCow::from(geometry_collection) + } + Geometry::Rect(rect) => GeometryCow::from(rect), + Geometry::Triangle(triangle) => GeometryCow::from(triangle), + } + } +} diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 66968ad15..f6ecfae66 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -222,6 +222,7 @@ pub use crate::algorithm::*; pub use crate::types::Closest; use std::cmp::Ordering; +pub use crate::relate::PreparedGeometry; pub use geo_types::{coord, line_string, point, polygon, wkt, CoordFloat, CoordNum}; pub mod geometry;