Merge pull request #1057 from georust/triangle

Point in triangle and rect performance improvements

geo-types/src/geometry/line_string.rs

@@ -257,15 +257,15 @@ impl<T: CoordNum> LineString<T> {
     /// );
     /// assert!(lines.next().is_none());
     /// ```
-    pub fn lines(&'_ self) -> impl ExactSizeIterator + Iterator<Item = Line<T>> + '_ {
+    pub fn lines(&'_ self) -> impl ExactSizeIterator<Item = Line<T>> + '_ {
         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)) }
         })
     }
 
     /// An iterator which yields the coordinates of a [`LineString`] as [Triangle]s
-    pub fn triangles(&'_ self) -> impl ExactSizeIterator + Iterator<Item = Triangle<T>> + '_ {
+    pub fn triangles(&'_ self) -> impl ExactSizeIterator<Item = Triangle<T>> + '_ {
         self.0.windows(3).map(|w| {
             // slice::windows(N) is guaranteed to yield a slice with exactly N elements
             unsafe {

geo/CHANGES.md

@@ -41,6 +41,8 @@
   implemented for `CoordsIter` without re-allocating (e.g., creating a `MultiPoint`).
 * Add `compose_many` method to `AffineOps`
   * <https://github.com/georust/geo/pull/1148>
+* Point in `Triangle` and `Rect` performance improvemnets
+  * <https://github.com/georust/geo/pull/1057>
 
 ## 0.27.0
 

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

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);

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<Point> = geo_test_fixtures::nl_plots()
+            .iter()
+            .map(|plot| plot.centroid().unwrap())
+            .collect();
+        let zone_bbox: Vec<Rect> = 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);

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<Point> = geo_test_fixtures::nl_plots()
+        .iter()
+        .map(|plot| plot.centroid().unwrap())
+        .collect();
+    let zone_bbox: Vec<Rect> = 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<Point> = geo_test_fixtures::nl_plots()
+        .iter()
+        .map(|plot| plot.centroid().unwrap())
+        .collect();
+    let zone_triangles: Vec<Triangle> = 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
+);

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<T>) -> 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
     }
 }
 

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;
+        }
     }
 }
 

geo/src/algorithm/intersects/collections.rs

@@ -25,6 +25,7 @@ where
 symmetric_intersects_impl!(Coord<T>, Geometry<T>);
 symmetric_intersects_impl!(Line<T>, Geometry<T>);
 symmetric_intersects_impl!(Rect<T>, Geometry<T>);
+symmetric_intersects_impl!(Triangle<T>, Geometry<T>);
 symmetric_intersects_impl!(Polygon<T>, Geometry<T>);
 
 impl<T, G> Intersects<G> for GeometryCollection<T>
@@ -43,4 +44,5 @@ where
 symmetric_intersects_impl!(Coord<T>, GeometryCollection<T>);
 symmetric_intersects_impl!(Line<T>, GeometryCollection<T>);
 symmetric_intersects_impl!(Rect<T>, GeometryCollection<T>);
+symmetric_intersects_impl!(Triangle<T>, GeometryCollection<T>);
 symmetric_intersects_impl!(Polygon<T>, GeometryCollection<T>);

geo/src/algorithm/intersects/line.rs

@@ -67,3 +67,13 @@ where
         }
     }
 }
+
+impl<T> Intersects<Triangle<T>> for Line<T>
+where
+    T: GeoNum,
+{
+    fn intersects(&self, rhs: &Triangle<T>) -> bool {
+        self.intersects(&rhs.to_polygon())
+    }
+}
+symmetric_intersects_impl!(Triangle<T>, Line<T>);