Merge #226

226: Polygon-Line euclidean distance fix r=urschrei a=urschrei See here: #219 (review) WIP until I add tests for a Line contained in a Polygon interior ring and intersection Co-authored-by: Stephan Hügel <stephan.hugel.12@ucl.ac.uk>
georust · May 16, 2018 · af8c95b · af8c95b
2 parents 29d7cdd + 8b71df6
commit af8c95b
Show file tree

Hide file tree

Showing 4 changed files with 112 additions and 35 deletions.
diff --git a/CHANGES.md b/CHANGES.md
@@ -1,5 +1,10 @@
 # Changes
 
+## geo 0.9.1 (unreleased)
+
+* Fix Line-Polygon euclidean distance
+  * <https://github.com/georust/rust-geo/pull/226>
+
 ## geo 0.9.0
 
 * Make serde an optional dependency for `geo`, rename feature to `use-serde`

diff --git a/geo/Cargo.toml b/geo/Cargo.toml
@@ -1,7 +1,7 @@
 [package]
 name = "geo"
 description = "Geospatial primitives and algorithms"
-version = "0.9.0"
+version = "0.9.1"
 authors = ["Corey Farwell <coreyf@rwell.org>"]
 license = "MIT/Apache-2.0"
 repository = "https://github.com/georust/rust-geo"

diff --git a/geo/src/algorithm/euclidean_distance.rs b/geo/src/algorithm/euclidean_distance.rs
@@ -298,18 +298,15 @@ where
     fn euclidean_distance(&self, other: &Polygon<T>) -> T {
         if self.intersects(other) || other.contains(self) {
             T::zero()
-        } else {
-            // still a possibility that the LineString's inside an interior ring
-            if !other.interiors.is_empty() && ring_contains_point(other, &self.0[0]) {
-                // check each ring distance, returning the minimum
-                let mut mindist: T = Float::max_value();
-                for ring in &other.interiors {
-                    mindist = mindist.min(nearest_neighbour_distance(self, ring))
-                }
-                mindist
-            } else {
-                nearest_neighbour_distance(self, &other.exterior)
+        } else if !other.interiors.is_empty() && ring_contains_point(other, &self.0[0]) {
+            // check each ring distance, returning the minimum
+            let mut mindist: T = Float::max_value();
+            for ring in &other.interiors {
+                mindist = mindist.min(nearest_neighbour_distance(self, ring))
             }
+            mindist
+        } else {
+            nearest_neighbour_distance(self, &other.exterior)
         }
     }
 }
@@ -324,18 +321,50 @@ where
     }
 }
 
-// Line to Polygon distance
-impl<T> EuclideanDistance<T, Polygon<T>> for Line<T>
+impl<T> EuclideanDistance<T, Line<T>> for Line<T>
 where
     T: Float + FloatConst + Signed + SpadeFloat,
 {
-    fn euclidean_distance(&self, other: &Polygon<T>) -> T {
+    fn euclidean_distance(&self, other: &Line<T>) -> T {
+        if self.intersects(other) || self.contains(other) {
+            return T::zero();
+        }
+        // minimum of two Point-Line distances
         self.start
             .euclidean_distance(other)
             .min(self.end.euclidean_distance(other))
     }
 }
 
+// Line to Polygon distance
+impl<T> EuclideanDistance<T, Polygon<T>> for Line<T>
+where
+    T: Float + Signed + SpadeFloat + FloatConst,
+{
+    fn euclidean_distance(&self, other: &Polygon<T>) -> T {
+        if other.contains(self) || self.intersects(other) {
+            return T::zero();
+        }
+        // point-line distance between each exterior polygon point and the line
+        let exterior_min = other.exterior.points().fold(T::max_value(), |acc, point| {
+            acc.min(self.euclidean_distance(point))
+        });
+        // point-line distance between each interior ring point and the line
+        // if there are no rings this just evaluates to max_float
+        let interior_min = other
+            .interiors
+            .iter()
+            .map(|ring| {
+                ring.lines().fold(T::max_value(), |acc, line| {
+                    acc.min(self.euclidean_distance(&line))
+                })
+            })
+            .fold(T::max_value(), |acc, ring_min| acc.min(ring_min));
+        // return smaller of the two values
+        exterior_min.min(interior_min)
+    }
+}
+
 // Polygon to Line distance
 impl<T> EuclideanDistance<T, Line<T>> for Polygon<T>
 where
@@ -373,7 +402,7 @@ where
             }
             return mindist;
         }
-        if !(self.is_convex() && !poly2.is_convex()) {
+        if poly2.is_convex() || !self.is_convex() {
             // fall back to R* nearest neighbour method
             nearest_neighbour_distance(&self.exterior, &poly2.exterior)
         } else {
@@ -384,7 +413,7 @@ where
 
 /// Uses an R* tree and nearest-neighbour lookups to calculate minimum distances
 // This is somewhat slow and memory-inefficient, but certainly better than quadratic time
-fn nearest_neighbour_distance<T>(geom1: &LineString<T>, geom2: &LineString<T>) -> T
+pub fn nearest_neighbour_distance<T>(geom1: &LineString<T>, geom2: &LineString<T>) -> T
 where
     T: Float + SpadeFloat,
 {
@@ -831,4 +860,49 @@ mod test {
         let in_ring_ls: LineString<f64> = in_ring.into();
         assert_eq!(ring_ls.euclidean_distance(&in_ring_ls), 5.992772737231033);
     }
+    #[test]
+    // Line-Polygon test: closest point on Polygon is NOT nearest to a Line end-point
+    fn test_line_polygon_simple() {
+        let line = Line::new(Point::new(0.0, 0.0), Point::new(0.0, 3.0));
+        let v = vec![
+            (5.0, 1.0),
+            (5.0, 2.0),
+            (0.25, 1.5),
+            (5.0, 1.0)
+        ];
+        let poly = Polygon::new(v.into(), vec![]);
+        assert_eq!(line.euclidean_distance(&poly), 0.25);
+    }
+    #[test]
+    // Line-Polygon test: Line intersects Polygon
+    fn test_line_polygon_intersects() {
+        let line = Line::new(Point::new(0.5, 0.0), Point::new(0.0, 3.0));
+        let v = vec![
+            (5.0, 1.0),
+            (5.0, 2.0),
+            (0.25, 1.5),
+            (5.0, 1.0)
+        ];
+        let poly = Polygon::new(v.into(), vec![]);
+        assert_eq!(line.euclidean_distance(&poly), 0.0);
+    }
+    #[test]
+    // Line-Polygon test: Line contained by interior ring
+    fn test_line_polygon_inside_ring() {
+        let line = Line::new(Point::new(4.4, 1.5), Point::new(4.45, 1.5));
+        let v = vec![
+            (5.0, 1.0),
+            (5.0, 2.0),
+            (0.25, 1.0),
+            (5.0, 1.0)
+        ];
+        let v2 = vec![
+            (4.5, 1.2),
+            (4.5, 1.8),
+            (3.5, 1.2),
+            (4.5, 1.2)
+        ];
+        let poly = Polygon::new(v.into(), vec![v2.into()]);
+        assert_eq!(line.euclidean_distance(&poly), 0.04999999999999982);
+    }
 }
diff --git a/geo/src/algorithm/polygon_distance_fast_path.rs b/geo/src/algorithm/polygon_distance_fast_path.rs
@@ -20,11 +20,11 @@ where
     let ymax2 = poly2.exterior.0[poly2_extremes.ymax];
 
     let mut state = Polydist {
-        poly1: poly1,
-        poly2: poly2,
+        poly1,
+        poly2,
         dist: T::infinity(),
-        ymin1: ymin1,
-        ymax2: ymax2,
+        ymin1,
+        ymax2,
         // initial polygon 1 min y idx
         p1_idx: poly1_extremes.ymin,
         // initial polygon 2 max y idx
@@ -345,23 +345,21 @@ where
             // implies p.x() < pprev.x()
             punit = Point::new(p.x(), p.y() + hundred);
         }
-    } else {
-        if p.x() > pprev.x() {
+    } else if p.x() > pprev.x() {
+        punit = Point::new(p.x(), p.y() + hundred);
+    } else if p.x() == pprev.x() {
+        if p.y() > pprev.y() {
             punit = Point::new(p.x(), p.y() + hundred);
-        } else if p.x() == pprev.x() {
-            if p.y() > pprev.y() {
-                punit = Point::new(p.x(), p.y() + hundred);
-            } else {
-                // implies p.y() < pprev.y()
-                // it's safe not to explicitly cover p.y() == pprev.y() because that
-                // implies that the x values are equal, and the y values are equal,
-                // and this is impossible
-                punit = Point::new(p.x(), p.y() - hundred);
-            }
         } else {
-            // implies p.x() < pprev.x()
+            // implies p.y() < pprev.y()
+            // it's safe not to explicitly cover p.y() == pprev.y() because that
+            // implies that the x values are equal, and the y values are equal,
+            // and this is impossible
             punit = Point::new(p.x(), p.y() - hundred);
         }
+    } else {
+        // implies p.x() < pprev.x()
+        punit = Point::new(p.x(), p.y() - hundred);
     }
     let triarea = triangle_area(p, &punit, &pnext);
     let edgelen = p.euclidean_distance(&pnext);
@@ -524,7 +522,7 @@ where
     }
     // 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_none() {
+    if state.p1 == state.ymin1 && state.q2 == state.ymax2 && state.start.is_some() {
         state.finished = true;
     } else {
         state.start = Some(false);