From 08b5d30fe5662d069ef1d34c4400b05da8031d17 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Fri, 18 Feb 2022 17:21:58 +0000 Subject: [PATCH 1/2] Fix fast path The fix replaces the previous, buggy iteration logic and base case with logic that relies upon the fact that this is a linear algorithm: the minimum distance can be calculated by spinning the calipers at most n times, where n is the number of edges in the larger of the two polygons. --- geo/src/algorithm/euclidean_distance.rs | 39 ++++++++++++++ .../algorithm/polygon_distance_fast_path.rs | 54 +++++++++++-------- 2 files changed, 71 insertions(+), 22 deletions(-) diff --git a/geo/src/algorithm/euclidean_distance.rs b/geo/src/algorithm/euclidean_distance.rs index e11572646..1243fad2f 100644 --- a/geo/src/algorithm/euclidean_distance.rs +++ b/geo/src/algorithm/euclidean_distance.rs @@ -576,6 +576,7 @@ mod test { use crate::algorithm::convex_hull::ConvexHull; use crate::algorithm::euclidean_distance::EuclideanDistance; use crate::{Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon}; + use geo_types::Rect; use geo_types::{polygon, private_utils::line_segment_distance, Coordinate}; #[test] @@ -1026,6 +1027,44 @@ mod test { assert_relative_eq!(dist2, 12.0); } #[test] + // test vertex-vertex minimum distance + fn test_minimum_polygon_distance_4() { + let poly1 = Rect::new((0., 0.), (1., 1.)).to_polygon(); + let poly2 = Rect::new((2., 2.), (3., 3.)).to_polygon(); + let dist = min_poly_dist(&poly1, &poly2); + let dist2 = nearest_neighbour_distance(poly1.exterior(), poly2.exterior()); + assert_eq!(dist, dist2); + } + #[test] + fn test_minimum_polygon_distance_5() { + let poly1 = polygon!( + exterior: [ + (x: -1.5350399, y: -0.18596762), + (x: -0.33726108, y: 0.2206358), + (x: -0.13836576, y: 0.42173427), + (x: -0.33946428, y: 0.6206298), + (x: -1.1405537, y: 0.8162205), + (x: -1.3405507, y: 0.81511897), + (x: -1.5372429, y: 0.21402647), + ], + interiors: [] + ); + + let poly2 = polygon!( + exterior: [ + (x: 0.6649223, y: -0.17384565), + (x: 0.8649193, y: -0.1727441), + (x: 1.0616114, y: 0.4283484), + (x: 0.8616146, y: 0.4272468), + (x: 0.66271913, y: 0.22614834), + ], + interiors: [] + ); + let dist = min_poly_dist(&poly1, &poly2); + let dist2 = nearest_neighbour_distance(poly1.exterior(), poly2.exterior()); + assert_eq!(dist, dist2); + } + #[test] fn test_large_polygon_distance() { let ls = geo_test_fixtures::norway_main::(); let poly1 = Polygon::new(ls, vec![]); diff --git a/geo/src/algorithm/polygon_distance_fast_path.rs b/geo/src/algorithm/polygon_distance_fast_path.rs index d04fb9b65..5b68d9464 100644 --- a/geo/src/algorithm/polygon_distance_fast_path.rs +++ b/geo/src/algorithm/polygon_distance_fast_path.rs @@ -23,8 +23,6 @@ where poly1, poly2, dist: T::infinity(), - ymin1, - ymax2, // initial polygon 1 min y idx p1_idx: poly1_extremes.y_min.index, // initial polygon 2 max y idx @@ -40,15 +38,23 @@ where ap1: T::zero(), aq2: T::zero(), start: None, - finished: false, ip1: false, iq2: false, slope: T::zero(), vertical: false, + // minimum distance can be calculated in at most this many iterations + // we only need to spin the calipers equal to the number of edges in the largest polygon + max_iterations: poly1 + .exterior() + .lines() + .len() + .max(poly2.exterior().lines().len()), }; - while !state.finished { + let mut iterations = 0usize; + while iterations <= state.max_iterations { nextpoints(&mut state); computemin(&mut state); + iterations += 1; } state.dist } @@ -94,9 +100,7 @@ where poly1: &'a Polygon, poly2: &'a Polygon, dist: T, - ymin1: Point, p1_idx: usize, - ymax2: Point, q2_idx: usize, p1: Point, q2: Point, @@ -108,11 +112,11 @@ where ap1: T, aq2: T, start: Option, - finished: bool, ip1: bool, iq2: bool, slope: T, vertical: bool, + max_iterations: usize, } // much of the following code is ported from Java, copyright 1999 Hormoz Pirzadeh, available at: @@ -411,7 +415,7 @@ fn nextpoints(state: &mut Polydist) where T: GeoFloat + FloatConst + Signed, { - state.alignment = None; + state.alignment = Some(AlignedEdge::VertexP); state.ip1 = false; state.iq2 = false; state.ap1 = vertex_line_angle( @@ -464,7 +468,7 @@ where state.alignment = match state.alignment { None => Some(AlignedEdge::VertexQ), Some(_) => Some(AlignedEdge::Edge), - } + }; } if state.ip1 { if state.p1.x() == state.p1next.x() { @@ -473,28 +477,29 @@ where state.slope = T::zero(); } else { state.vertical = false; - state.slope = Line::new(state.p1next.0, state.p1.0).slope(); + if state.p1.x() > state.p1next.x() { + state.slope = (state.p1.y() - state.p1next.y()) / (state.p1.x() - state.p1next.x()); + } else { + state.slope = (state.p1next.y() - state.p1.y()) / (state.p1next.x() - state.p1.x()); + } } - } - if state.iq2 { + } else if state.iq2 { if state.q2.x() == state.q2next.x() { // The Q line of support is vertical state.vertical = true; state.slope = T::zero(); } else { state.vertical = false; - state.slope = Line::new(state.q2next.0, state.q2.0).slope(); + if state.q2.x() > state.q2next.x() { + state.slope = (state.q2.y() - state.q2next.y()) / (state.q2.x() - state.q2next.x()); + } else { + state.slope = (state.q2next.y() - state.q2.y()) / (state.q2next.x() - state.q2.x()); + } } } - // A start value's been set, and both polygon indices are in their initial - // positions -- we're finished, so return the minimum distance - if state.p1 == state.ymin1 && state.q2 == state.ymax2 && state.start.is_some() { - state.finished = true; - } else { - state.start = Some(false); - state.p1 = state.p1next; - state.q2 = state.q2next; - } + state.start = Some(false); + state.p1 = state.p1next; + state.q2 = state.q2next; } /// compute the minimum distance between entities (edges or vertices) @@ -567,6 +572,11 @@ where } Some(AlignedEdge::Edge) => { // both lines of support coincide with edges (i.e. they're parallel) + newdist = state.p1.euclidean_distance(&state.q2); + if newdist <= state.dist { + // New minimum distance is between p1 and q2 + state.dist = newdist; + } newdist = state.p1.euclidean_distance(&state.q2prev); if newdist <= state.dist { // New minimum distance is between p1 and q2prev From 0072cfe17ac53229e8d87bd15526d4a087c06e3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stephan=20H=C3=BCgel?= Date: Sat, 19 Feb 2022 20:52:21 +0000 Subject: [PATCH 2/2] Fix max iterations --- geo/src/algorithm/polygon_distance_fast_path.rs | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/geo/src/algorithm/polygon_distance_fast_path.rs b/geo/src/algorithm/polygon_distance_fast_path.rs index 5b68d9464..9a0aa8f0d 100644 --- a/geo/src/algorithm/polygon_distance_fast_path.rs +++ b/geo/src/algorithm/polygon_distance_fast_path.rs @@ -43,12 +43,9 @@ where slope: T::zero(), vertical: false, // minimum distance can be calculated in at most this many iterations - // we only need to spin the calipers equal to the number of edges in the largest polygon - max_iterations: poly1 - .exterior() - .lines() - .len() - .max(poly2.exterior().lines().len()), + // we only need to spin the calipers equal to the total number of vertices in both polygons + // alternatively, we could accumulate the total rotation angle and stop when it exceeds 2pi + max_iterations: poly1.exterior().0.len() + poly2.exterior().0.len(), }; let mut iterations = 0usize; while iterations <= state.max_iterations {