diff --git a/geo-types/src/geometry/line_string.rs b/geo-types/src/geometry/line_string.rs index 05713ebfc..4512b4ccc 100644 --- a/geo-types/src/geometry/line_string.rs +++ b/geo-types/src/geometry/line_string.rs @@ -257,7 +257,7 @@ impl LineString { /// ); /// assert!(lines.next().is_none()); /// ``` - pub fn lines(&'_ self) -> impl ExactSizeIterator + Iterator> + '_ { + pub fn lines(&'_ self) -> impl ExactSizeIterator> + '_ { self.0.windows(2).map(|w| { // slice::windows(N) is guaranteed to yield a slice with exactly N elements unsafe { Line::new(*w.get_unchecked(0), *w.get_unchecked(1)) } @@ -265,7 +265,7 @@ impl LineString { } /// An iterator which yields the coordinates of a [`LineString`] as [Triangle]s - pub fn triangles(&'_ self) -> impl ExactSizeIterator + Iterator> + '_ { + pub fn triangles(&'_ self) -> impl ExactSizeIterator> + '_ { self.0.windows(3).map(|w| { // slice::windows(N) is guaranteed to yield a slice with exactly N elements unsafe { diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 80630b3f5..d7f29d2bd 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -41,6 +41,8 @@ implemented for `CoordsIter` without re-allocating (e.g., creating a `MultiPoint`). * Add `compose_many` method to `AffineOps` * +* Point in `Triangle` and `Rect` performance improvemnets + * ## 0.27.0 diff --git a/geo/Cargo.toml b/geo/Cargo.toml index f38fb5cbc..45f5f6f4f 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -47,6 +47,10 @@ wkt = "0.10.1" name = "area" harness = false +[[bench]] +name = "coordinate_position" +harness = false + [[bench]] name = "contains" harness = false diff --git a/geo/benches/contains.rs b/geo/benches/contains.rs index 56080b735..3a09d615a 100644 --- a/geo/benches/contains.rs +++ b/geo/benches/contains.rs @@ -1,6 +1,6 @@ use criterion::{criterion_group, criterion_main, Criterion}; use geo::contains::Contains; -use geo::{point, polygon, Line, Point, Polygon}; +use geo::{point, polygon, Line, Point, Polygon, Triangle}; fn criterion_benchmark(c: &mut Criterion) { c.bench_function("point in simple polygon", |bencher| { @@ -110,6 +110,24 @@ fn criterion_benchmark(c: &mut Criterion) { ); }); }); + + c.bench_function("Triangle contains point", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (5., 10.)]); + let point = Point::new(5., 5.); + + bencher.iter(|| { + assert!(criterion::black_box(&triangle).contains(criterion::black_box(&point))); + }); + }); + + c.bench_function("Triangle contains point on edge", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (6., 10.)]); + let point = Point::new(3., 5.); + + bencher.iter(|| { + assert!(!criterion::black_box(&triangle).contains(criterion::black_box(&point))); + }); + }); } criterion_group!(benches, criterion_benchmark); diff --git a/geo/benches/coordinate_position.rs b/geo/benches/coordinate_position.rs new file mode 100644 index 000000000..2b5e2b468 --- /dev/null +++ b/geo/benches/coordinate_position.rs @@ -0,0 +1,69 @@ +#[macro_use] +extern crate criterion; +extern crate geo; + +use geo::{ + coordinate_position::CoordPos, BoundingRect, Centroid, CoordinatePosition, Point, Rect, + Triangle, +}; + +use criterion::Criterion; + +fn criterion_benchmark(c: &mut Criterion) { + c.bench_function("Point position to rect", |bencher| { + let plot_centroids: Vec = geo_test_fixtures::nl_plots() + .iter() + .map(|plot| plot.centroid().unwrap()) + .collect(); + let zone_bbox: Vec = geo_test_fixtures::nl_zones() + .iter() + .map(|plot| plot.bounding_rect().unwrap()) + .collect(); + bencher.iter(|| { + let mut inside = 0; + let mut outsied = 0; + let mut boundary = 0; + + for a in &plot_centroids { + for b in &zone_bbox { + match criterion::black_box(b).coordinate_position(criterion::black_box(&a.0)) { + CoordPos::OnBoundary => boundary += 1, + CoordPos::Inside => inside += 1, + CoordPos::Outside => outsied += 1, + } + } + } + + assert_eq!(inside, 2246); + assert_eq!(outsied, 26510); + assert_eq!(boundary, 0); + }); + }); + + c.bench_function("Point in triangle", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (5., 10.)]); + let point = Point::new(5., 5.); + + bencher.iter(|| { + assert!( + criterion::black_box(&triangle).coordinate_position(criterion::black_box(&point.0)) + != CoordPos::Outside + ); + }); + }); + + c.bench_function("Point on triangle boundary", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (6., 10.)]); + let point = Point::new(3., 5.); + + bencher.iter(|| { + assert!( + criterion::black_box(&triangle).coordinate_position(criterion::black_box(&point.0)) + == CoordPos::OnBoundary + ); + }); + }); +} + +criterion_group!(benches, criterion_benchmark); +criterion_main!(benches); diff --git a/geo/benches/intersection.rs b/geo/benches/intersection.rs index f0bcaeeba..8dda67d0f 100644 --- a/geo/benches/intersection.rs +++ b/geo/benches/intersection.rs @@ -60,11 +60,110 @@ fn rect_intersection(c: &mut Criterion) { }); } +fn point_rect_intersection(c: &mut Criterion) { + use geo::algorithm::{BoundingRect, Centroid}; + use geo::geometry::{Point, Rect}; + let plot_centroids: Vec = geo_test_fixtures::nl_plots() + .iter() + .map(|plot| plot.centroid().unwrap()) + .collect(); + let zone_bbox: Vec = geo_test_fixtures::nl_zones() + .iter() + .map(|plot| plot.bounding_rect().unwrap()) + .collect(); + + c.bench_function("Point intersects rect", |bencher| { + bencher.iter(|| { + let mut intersects = 0; + let mut non_intersects = 0; + + for a in &plot_centroids { + for b in &zone_bbox { + if criterion::black_box(a.intersects(b)) { + intersects += 1; + } else { + non_intersects += 1; + } + } + } + + assert_eq!(intersects, 2246); + assert_eq!(non_intersects, 26510); + }); + }); +} + +fn point_triangle_intersection(c: &mut Criterion) { + use geo::{Centroid, TriangulateEarcut}; + use geo_types::{Point, Triangle}; + let plot_centroids: Vec = geo_test_fixtures::nl_plots() + .iter() + .map(|plot| plot.centroid().unwrap()) + .collect(); + let zone_triangles: Vec = geo_test_fixtures::nl_zones() + .iter() + .flat_map(|plot| plot.earcut_triangles_iter()) + .collect(); + + c.bench_function("Point intersects triangle", |bencher| { + bencher.iter(|| { + let mut intersects = 0; + let mut non_intersects = 0; + + for a in &plot_centroids { + for b in &zone_triangles { + if criterion::black_box(a.intersects(b)) { + intersects += 1; + } else { + non_intersects += 1; + } + } + } + + assert_eq!(intersects, 533); + assert_eq!(non_intersects, 5450151); + }); + }); + + c.bench_function("Triangle intersects point", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (5., 10.)]); + let point = Point::new(5., 5.); + + bencher.iter(|| { + assert!(criterion::black_box(&triangle).intersects(criterion::black_box(&point))); + }); + }); + + c.bench_function("Triangle intersects point on edge", |bencher| { + let triangle = Triangle::from([(0., 0.), (10., 0.), (6., 10.)]); + let point = Point::new(3., 5.); + + bencher.iter(|| { + assert!(criterion::black_box(&triangle).intersects(criterion::black_box(&point))); + }); + }); +} + criterion_group! { name = bench_multi_polygons; config = Criterion::default().sample_size(10); targets = multi_polygon_intersection } criterion_group!(bench_rects, rect_intersection); +criterion_group! { + name = bench_point_rect; + config = Criterion::default().sample_size(50); + targets = point_rect_intersection +} +criterion_group! { + name = bench_point_triangle; + config = Criterion::default().sample_size(50); + targets = point_triangle_intersection +} -criterion_main!(bench_multi_polygons, bench_rects); +criterion_main!( + bench_multi_polygons, + bench_rects, + bench_point_rect, + bench_point_triangle +); diff --git a/geo/src/algorithm/contains/triangle.rs b/geo/src/algorithm/contains/triangle.rs index b1500bb65..e21b079b0 100644 --- a/geo/src/algorithm/contains/triangle.rs +++ b/geo/src/algorithm/contains/triangle.rs @@ -1,6 +1,6 @@ use super::{impl_contains_from_relate, impl_contains_geometry_for, Contains}; use crate::geometry::*; -use crate::{GeoFloat, GeoNum}; +use crate::{kernels::Kernel, GeoFloat, GeoNum, Orientation}; // ┌──────────────────────────────┐ // │ Implementations for Triangle │ @@ -11,9 +11,31 @@ where T: GeoNum, { fn contains(&self, coord: &Coord) -> bool { - let ls = LineString::new(vec![self.0, self.1, self.2, self.0]); - use crate::utils::{coord_pos_relative_to_ring, CoordPos}; - coord_pos_relative_to_ring(*coord, &ls) == CoordPos::Inside + // leverageing robust predicates + self.to_lines() + .map(|l| T::Ker::orient2d(l.start, l.end, *coord)) + .windows(2) + .all(|win| win[0] == win[1] && win[0] != Orientation::Collinear) + + // // neglecting robust prdicates, hence faster + // let p0x = self.0.x.to_f64().unwrap(); + // let p0y = self.0.y.to_f64().unwrap(); + // let p1x = self.1.x.to_f64().unwrap(); + // let p1y = self.1.y.to_f64().unwrap(); + // let p2x = self.2.x.to_f64().unwrap(); + // let p2y = self.2.y.to_f64().unwrap(); + + // let px = coord.x.to_f64().unwrap(); + // let py = coord.y.to_f64().unwrap(); + + // let a = 0.5 * (-p1y * p2x + p0y * (-p1x + p2x) + p0x * (p1y - p2y) + p1x * p2y); + + // let sign = a.signum(); + + // let s = (p0y * p2x - p0x * p2y + (p2y - p0y) * px + (p0x - p2x) * py) * sign; + // let t = (p0x * p1y - p0y * p1x + (p0y - p1y) * px + (p1x - p0x) * py) * sign; + + // s > 0. && t > 0. && (s + t) < 2. * a * sign } } diff --git a/geo/src/algorithm/coordinate_position.rs b/geo/src/algorithm/coordinate_position.rs index ee0de3ea9..91f60e0c2 100644 --- a/geo/src/algorithm/coordinate_position.rs +++ b/geo/src/algorithm/coordinate_position.rs @@ -1,5 +1,7 @@ +use std::cmp::Ordering; + use crate::geometry::*; -use crate::intersects::value_in_between; +use crate::intersects::{point_in_rect, value_in_between}; use crate::kernels::*; use crate::{BoundingRect, HasDimensions, Intersects}; use crate::{GeoNum, GeometryCow}; @@ -184,9 +186,20 @@ where is_inside: &mut bool, boundary_count: &mut usize, ) { - // PERF TODO: I'm sure there's a better way to calculate than converting to a polygon - self.to_polygon() - .calculate_coordinate_position(coord, is_inside, boundary_count); + *is_inside = self + .to_lines() + .map(|l| { + let orientation = T::Ker::orient2d(l.start, l.end, *coord); + if orientation == Orientation::Collinear + && point_in_rect(*coord, l.start, l.end) + && coord.x != l.end.x + { + *boundary_count += 1; + } + orientation + }) + .windows(2) + .all(|win| win[0] == win[1] && win[0] != Orientation::Collinear); } } @@ -201,9 +214,39 @@ where is_inside: &mut bool, boundary_count: &mut usize, ) { - // PERF TODO: I'm sure there's a better way to calculate than converting to a polygon - self.to_polygon() - .calculate_coordinate_position(coord, is_inside, boundary_count); + let mut boundary = false; + + let min = self.min(); + + match coord.x.partial_cmp(&min.x).unwrap() { + Ordering::Less => return, + Ordering::Equal => boundary = true, + Ordering::Greater => {} + } + match coord.y.partial_cmp(&min.y).unwrap() { + Ordering::Less => return, + Ordering::Equal => boundary = true, + Ordering::Greater => {} + } + + let max = self.max(); + + match max.x.partial_cmp(&coord.x).unwrap() { + Ordering::Less => return, + Ordering::Equal => boundary = true, + Ordering::Greater => {} + } + match max.y.partial_cmp(&coord.y).unwrap() { + Ordering::Less => return, + Ordering::Equal => boundary = true, + Ordering::Greater => {} + } + + if boundary { + *boundary_count += 1; + } else { + *is_inside = true; + } } } diff --git a/geo/src/algorithm/intersects/collections.rs b/geo/src/algorithm/intersects/collections.rs index 877aff641..bdd426502 100644 --- a/geo/src/algorithm/intersects/collections.rs +++ b/geo/src/algorithm/intersects/collections.rs @@ -25,6 +25,7 @@ where symmetric_intersects_impl!(Coord, Geometry); symmetric_intersects_impl!(Line, Geometry); symmetric_intersects_impl!(Rect, Geometry); +symmetric_intersects_impl!(Triangle, Geometry); symmetric_intersects_impl!(Polygon, Geometry); impl Intersects for GeometryCollection @@ -43,4 +44,5 @@ where symmetric_intersects_impl!(Coord, GeometryCollection); symmetric_intersects_impl!(Line, GeometryCollection); symmetric_intersects_impl!(Rect, GeometryCollection); +symmetric_intersects_impl!(Triangle, GeometryCollection); symmetric_intersects_impl!(Polygon, GeometryCollection); diff --git a/geo/src/algorithm/intersects/line.rs b/geo/src/algorithm/intersects/line.rs index 95d39338e..4dc88701f 100644 --- a/geo/src/algorithm/intersects/line.rs +++ b/geo/src/algorithm/intersects/line.rs @@ -67,3 +67,13 @@ where } } } + +impl Intersects> for Line +where + T: GeoNum, +{ + fn intersects(&self, rhs: &Triangle) -> bool { + self.intersects(&rhs.to_polygon()) + } +} +symmetric_intersects_impl!(Triangle, Line); diff --git a/geo/src/algorithm/intersects/line_string.rs b/geo/src/algorithm/intersects/line_string.rs index 1d885bfe9..07f278c2d 100644 --- a/geo/src/algorithm/intersects/line_string.rs +++ b/geo/src/algorithm/intersects/line_string.rs @@ -19,6 +19,7 @@ where symmetric_intersects_impl!(Coord, LineString); symmetric_intersects_impl!(Line, LineString); symmetric_intersects_impl!(Rect, LineString); +symmetric_intersects_impl!(Triangle, LineString); // Blanket implementation from LineString impl Intersects for MultiLineString @@ -38,3 +39,4 @@ where symmetric_intersects_impl!(Point, MultiLineString); symmetric_intersects_impl!(Line, MultiLineString); symmetric_intersects_impl!(Rect, MultiLineString); +symmetric_intersects_impl!(Triangle, MultiLineString); diff --git a/geo/src/algorithm/intersects/point.rs b/geo/src/algorithm/intersects/point.rs index e0ad230fa..6c0cef52c 100644 --- a/geo/src/algorithm/intersects/point.rs +++ b/geo/src/algorithm/intersects/point.rs @@ -25,4 +25,5 @@ where symmetric_intersects_impl!(Coord, MultiPoint); symmetric_intersects_impl!(Line, MultiPoint); +symmetric_intersects_impl!(Triangle, MultiPoint); symmetric_intersects_impl!(Polygon, MultiPoint); diff --git a/geo/src/algorithm/intersects/polygon.rs b/geo/src/algorithm/intersects/polygon.rs index 676a9708e..df68b5bd8 100644 --- a/geo/src/algorithm/intersects/polygon.rs +++ b/geo/src/algorithm/intersects/polygon.rs @@ -1,8 +1,9 @@ use super::{has_disjoint_bboxes, Intersects}; -use crate::utils::{coord_pos_relative_to_ring, CoordPos}; -use crate::BoundingRect; +use crate::coordinate_position::CoordPos; +use crate::{BoundingRect, CoordinatePosition}; use crate::{ Coord, CoordNum, GeoNum, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Rect, + Triangle, }; impl Intersects> for Polygon @@ -10,11 +11,7 @@ where T: GeoNum, { fn intersects(&self, p: &Coord) -> bool { - coord_pos_relative_to_ring(*p, self.exterior()) != CoordPos::Outside - && self - .interiors() - .iter() - .all(|int| coord_pos_relative_to_ring(*p, int) != CoordPos::Inside) + self.coordinate_position(p) != CoordPos::Outside } } symmetric_intersects_impl!(Coord, Polygon); @@ -45,6 +42,16 @@ where } symmetric_intersects_impl!(Rect, Polygon); +impl Intersects> for Polygon +where + T: GeoNum, +{ + fn intersects(&self, rect: &Triangle) -> bool { + self.intersects(&rect.to_polygon()) + } +} +symmetric_intersects_impl!(Triangle, Polygon); + impl Intersects> for Polygon where T: GeoNum, @@ -81,6 +88,7 @@ where symmetric_intersects_impl!(Point, MultiPolygon); symmetric_intersects_impl!(Line, MultiPolygon); symmetric_intersects_impl!(Rect, MultiPolygon); +symmetric_intersects_impl!(Triangle, MultiPolygon); symmetric_intersects_impl!(Polygon, MultiPolygon); #[cfg(test)] diff --git a/geo/src/algorithm/intersects/rect.rs b/geo/src/algorithm/intersects/rect.rs index adb17b55d..860af707a 100644 --- a/geo/src/algorithm/intersects/rect.rs +++ b/geo/src/algorithm/intersects/rect.rs @@ -1,4 +1,4 @@ -use super::{value_in_range, Intersects}; +use super::Intersects; use crate::*; impl Intersects> for Rect @@ -6,11 +6,10 @@ where T: CoordNum, { fn intersects(&self, rhs: &Coord) -> bool { - // Funnily, we don't use point_in_rect, as we know - // self.min <= self.max. - let bound_1 = self.min(); - let bound_2 = self.max(); - value_in_range(rhs.x, bound_1.x, bound_2.x) && value_in_range(rhs.y, bound_1.y, bound_2.y) + rhs.x >= self.min().x + && rhs.y >= self.min().y + && rhs.x <= self.max().x + && rhs.y <= self.max().y } } symmetric_intersects_impl!(Coord, Rect); @@ -63,3 +62,13 @@ where } } symmetric_intersects_impl!(Line, Rect); + +impl Intersects> for Rect +where + T: GeoNum, +{ + fn intersects(&self, rhs: &Triangle) -> bool { + self.intersects(&rhs.to_polygon()) + } +} +symmetric_intersects_impl!(Triangle, Rect); diff --git a/geo/src/algorithm/intersects/triangle.rs b/geo/src/algorithm/intersects/triangle.rs index 0af84aa83..823953320 100644 --- a/geo/src/algorithm/intersects/triangle.rs +++ b/geo/src/algorithm/intersects/triangle.rs @@ -1,16 +1,52 @@ use super::Intersects; use crate::*; -impl Intersects for Triangle +impl Intersects> for Triangle where - T: CoordNum, - Polygon: Intersects, + T: GeoNum, { - fn intersects(&self, rhs: &G) -> bool { - self.to_polygon().intersects(rhs) + fn intersects(&self, rhs: &Coord) -> bool { + let mut orientations = self + .to_lines() + .map(|l| T::Ker::orient2d(l.start, l.end, *rhs)); + + orientations.sort(); + + !orientations + .windows(2) + .any(|win| win[0] != win[1] && win[1] != Orientation::Collinear) + + // // neglecting robust predicates, hence faster + // let p0x = self.0.x.to_f64().unwrap(); + // let p0y = self.0.y.to_f64().unwrap(); + // let p1x = self.1.x.to_f64().unwrap(); + // let p1y = self.1.y.to_f64().unwrap(); + // let p2x = self.2.x.to_f64().unwrap(); + // let p2y = self.2.y.to_f64().unwrap(); + + // let px = rhs.x.to_f64().unwrap(); + // let py = rhs.y.to_f64().unwrap(); + + // let s = (p0x - p2x) * (py - p2y) - (p0y - p2y) * (px - p2x); + // let t = (p1x - p0x) * (py - p0y) - (p1y - p0y) * (px - p0x); + + // if (s < 0.) != (t < 0.) && s != 0. && t != 0. { + // return false; + // } + + // let d = (p2x - p1x) * (py - p1y) - (p2y - p1y) * (px - p1x); + // d == 0. || (d < 0.) == (s + t <= 0.) } } + symmetric_intersects_impl!(Coord, Triangle); -symmetric_intersects_impl!(Line, Triangle); -symmetric_intersects_impl!(Rect, Triangle); -symmetric_intersects_impl!(Polygon, Triangle); +symmetric_intersects_impl!(Triangle, Point); + +impl Intersects> for Triangle +where + T: GeoNum, +{ + fn intersects(&self, rhs: &Triangle) -> bool { + self.to_polygon().intersects(&rhs.to_polygon()) + } +}