Merge #1004

1004: Performance fix for point in polygon r=michaelkirk a=b4l

- [x] I agree to follow the project's [code of conduct](https://github.com/georust/geo/blob/main/CODE_OF_CONDUCT.md).
- [ ] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---


Co-authored-by: Balthasar Teuscher <[email protected]>
Co-authored-by: Michael Kirk <[email protected]>

geo/benches/contains.rs

@@ -43,15 +43,57 @@ fn criterion_benchmark(c: &mut Criterion) {
         });
     });
 
-    c.bench_function("point outside complex polygon", |bencher| {
+    c.bench_function(
+        "point outside, but within bbox, of complex polygon",
+        |bencher| {
+            let polygon = Polygon::<f64>::new(geo_test_fixtures::louisiana(), vec![]);
+            // lake borgne - near and mostly surrounded by, but not inside of, Louisiana
+            let point = point!(x: -89.641854, y: 30.026283);
+            bencher.iter(|| {
+                assert!(!criterion::black_box(&polygon).contains(criterion::black_box(&point)));
+            });
+        },
+    );
+
+    c.bench_function("point outside bbox of complex polygon", |bencher| {
         let polygon = Polygon::<f64>::new(geo_test_fixtures::louisiana(), vec![]);
-        // lake borgne - near and mostly surrounded by, but not inside of, Louisiana
-        let point = point!(x: -89.641854, y: 30.026283);
+        let point = point!(x: 2.3522, y: 48.8566);
         bencher.iter(|| {
             assert!(!criterion::black_box(&polygon).contains(criterion::black_box(&point)));
         });
     });
 
+    c.bench_function(
+        "point horizontal to comb teeth aka bart's haircut",
+        |bencher| {
+            // Testing a pathological case where the point is horizontal to lots of edges.
+            //
+            // comb teeth -> |\/\/\/\/\/\/|      *  <---point
+            //               |____________|
+            let polygon = polygon!(
+                (x: 0 ,y: 0),
+                (x: 0 ,y: 10),
+                (x: 1 ,y: 5),
+                (x: 2 ,y: 10),
+                (x: 3 ,y: 5),
+                (x: 4 ,y: 10),
+                (x: 5 ,y: 5),
+                (x: 6 ,y: 10),
+                (x: 7 ,y: 5),
+                (x: 8 ,y: 10),
+                (x: 9 ,y: 10),
+                (x: 10,y:  10),
+                (x: 10,y:  0),
+                (x: 0 ,y: 0)
+            );
+            let point = point!(x: 20, y: 7);
+
+            bencher.iter(|| {
+                assert!(!criterion::black_box(&polygon).contains(criterion::black_box(&point)));
+            })
+        },
+    );
+
     c.bench_function("line across complex polygon", |bencher| {
         let polygon = Polygon::<f64>::new(geo_test_fixtures::louisiana(), vec![]);
         // crossing part of, but not contained by Louisiana

geo/src/algorithm/coordinate_position.rs

@@ -1,5 +1,5 @@
 use crate::geometry::*;
-use crate::intersects::point_in_rect;
+use crate::intersects::value_in_between;
 use crate::kernels::*;
 use crate::{BoundingRect, HasDimensions, Intersects};
 use crate::{GeoNum, GeometryCow};
@@ -239,11 +239,6 @@ where
             return;
         }
 
-        // Ok to `unwrap` since we know `self` is non-empty, so bounding-rect is non-null
-        if !self.bounding_rect().unwrap().intersects(coord) {
-            return;
-        }
-
         match coord_pos_relative_to_ring(*coord, self.exterior()) {
             CoordPos::Outside => {}
             CoordPos::OnBoundary => {
@@ -371,21 +366,30 @@ where
     // See: https://en.wikipedia.org/wiki/Point_in_polygon#Winding_number_algorithm
     let mut winding_number = 0;
     for line in linestring.lines() {
+        // Edge Crossing Rules:
+        //   1. an upward edge includes its starting endpoint, and excludes its final endpoint;
+        //   2. a downward edge excludes its starting endpoint, and includes its final endpoint;
+        //   3. horizontal edges are excluded
+        //   4. the edge-ray intersection point must be strictly right of the coord.
         if line.start.y <= coord.y {
-            let o = T::Ker::orient2d(line.start, line.end, coord);
-            if o == Orientation::Collinear && point_in_rect(coord, line.start, line.end) {
-                return CoordPos::OnBoundary;
-            }
-            if line.end.y > coord.y && o == Orientation::CounterClockwise {
-                winding_number += 1;
-            }
+            if line.end.y >= coord.y {
+                let o = T::Ker::orient2d(line.start, line.end, coord);
+                if o == Orientation::CounterClockwise && line.end.y != coord.y {
+                    winding_number += 1
+                } else if o == Orientation::Collinear
+                    && value_in_between(coord.x, line.start.x, line.end.x)
+                {
+                    return CoordPos::OnBoundary;
+                }
+            };
         } else if line.end.y <= coord.y {
             let o = T::Ker::orient2d(line.start, line.end, coord);
-            if o == Orientation::Collinear && point_in_rect(coord, line.start, line.end) {
-                return CoordPos::OnBoundary;
-            }
             if o == Orientation::Clockwise {
-                winding_number -= 1;
+                winding_number -= 1
+            } else if o == Orientation::Collinear
+                && value_in_between(coord.x, line.start.x, line.end.x)
+            {
+                return CoordPos::OnBoundary;
             }
         }
     }

geo/src/algorithm/intersects/mod.rs

@@ -86,7 +86,7 @@ where
 // Helper function to check value lies between two bounds,
 // where the ordering of the bounds is not known
 #[inline]
-fn value_in_between<T>(value: T, bound_1: T, bound_2: T) -> bool
+pub(crate) fn value_in_between<T>(value: T, bound_1: T, bound_2: T) -> bool
 where
     T: std::cmp::PartialOrd,
 {