diff --git a/geo/benches/contains.rs b/geo/benches/contains.rs index 327615848..564b22441 100644 --- a/geo/benches/contains.rs +++ b/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::::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::::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::::new(geo_test_fixtures::louisiana(), vec![]); // crossing part of, but not contained by Louisiana diff --git a/geo/src/algorithm/coordinate_position.rs b/geo/src/algorithm/coordinate_position.rs index 4c89291c0..ee0de3ea9 100644 --- a/geo/src/algorithm/coordinate_position.rs +++ b/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; } } } diff --git a/geo/src/algorithm/intersects/mod.rs b/geo/src/algorithm/intersects/mod.rs index 09931a4be..d015835f4 100644 --- a/geo/src/algorithm/intersects/mod.rs +++ b/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(value: T, bound_1: T, bound_2: T) -> bool +pub(crate) fn value_in_between(value: T, bound_1: T, bound_2: T) -> bool where T: std::cmp::PartialOrd, {