Skip to content

Commit

Permalink
Merge georust#731
Browse files Browse the repository at this point in the history
731: Fix polygon distance fast path r=frewsxcv a=urschrei

This replaces the previous, buggy iteration logic: we now rely upon the fact that this is a linear algorithm (thank you for the inspiration for the fix `@msalib)` so the minimum distance can be calculated by spinning the calipers at most n times, where n is the <s>number of edges in the larger of the two polygons</s> total number of vertices in both polygons' exteriors.

Other improvements include re-calculating the slope according to the same logic in the original algorithm, fixing a branching difference in our port, and new tests that cover both the original and minimised errors that were found.

# THAT HAVING BEEN SAID:

This algorithm's implementation is **extremely** magical and if it were proposed for addition today I would be dubious about its inclusion: why does the `100.00` value work in various functions, but `float::max_value()` doesn't? Why is the epsilon (`0.002`) comparatively large, and why does substituting `float::min_value()` not work? Nobody knows.

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



Co-authored-by: Stephan Hügel <shugel@tcd.ie>
  • Loading branch information
bors[bot] and urschrei authored Feb 20, 2022
2 parents 497d67d + 0072cfe commit 1d3b329
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 22 deletions.
39 changes: 39 additions & 0 deletions geo/src/algorithm/euclidean_distance.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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::<f64>();
let poly1 = Polygon::new(ls, vec![]);
Expand Down
51 changes: 29 additions & 22 deletions geo/src/algorithm/polygon_distance_fast_path.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -40,15 +38,20 @@ 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 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(),
};
while !state.finished {
let mut iterations = 0usize;
while iterations <= state.max_iterations {
nextpoints(&mut state);
computemin(&mut state);
iterations += 1;
}
state.dist
}
Expand Down Expand Up @@ -94,9 +97,7 @@ where
poly1: &'a Polygon<T>,
poly2: &'a Polygon<T>,
dist: T,
ymin1: Point<T>,
p1_idx: usize,
ymax2: Point<T>,
q2_idx: usize,
p1: Point<T>,
q2: Point<T>,
Expand All @@ -108,11 +109,11 @@ where
ap1: T,
aq2: T,
start: Option<bool>,
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:
Expand Down Expand Up @@ -411,7 +412,7 @@ fn nextpoints<T>(state: &mut Polydist<T>)
where
T: GeoFloat + FloatConst + Signed,
{
state.alignment = None;
state.alignment = Some(AlignedEdge::VertexP);
state.ip1 = false;
state.iq2 = false;
state.ap1 = vertex_line_angle(
Expand Down Expand Up @@ -464,7 +465,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() {
Expand All @@ -473,28 +474,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)
Expand Down Expand Up @@ -567,6 +569,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
Expand Down

0 comments on commit 1d3b329

Please sign in to comment.