geo/algorithm/
euclidean_distance.rs

1use crate::utils::{coord_pos_relative_to_ring, CoordPos};
2use crate::EuclideanLength;
3use crate::Intersects;
4use crate::{
5    Coord, GeoFloat, GeoNum, Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point,
6    Polygon, Triangle,
7};
8use num_traits::{float::FloatConst, Bounded, Float, Signed};
9
10use rstar::RTree;
11use rstar::RTreeNum;
12
13/// Returns the distance between two geometries.
14
15pub trait EuclideanDistance<T, Rhs = Self> {
16    /// Returns the distance between two geometries
17    ///
18    /// If a `Point` is contained by a `Polygon`, the distance is `0.0`
19    ///
20    /// If a `Point` lies on a `Polygon`'s exterior or interior rings, the distance is `0.0`
21    ///
22    /// If a `Point` lies on a `LineString`, the distance is `0.0`
23    ///
24    /// The distance between a `Point` and an empty `LineString` is `0.0`
25    ///
26    /// # Examples
27    ///
28    /// `Point` to `Point`:
29    ///
30    /// ```
31    /// use approx::assert_relative_eq;
32    /// use geo::EuclideanDistance;
33    /// use geo::point;
34    ///
35    /// let p1 = point!(x: -72.1235, y: 42.3521);
36    /// let p2 = point!(x: -72.1260, y: 42.45);
37    ///
38    /// let distance = p1.euclidean_distance(&p2);
39    ///
40    /// assert_relative_eq!(distance, 0.09793191512474639);
41    /// ```
42    ///
43    /// `Point` to `Polygon`:
44    ///
45    /// ```
46    /// use approx::assert_relative_eq;
47    /// use geo::EuclideanDistance;
48    /// use geo::{point, polygon};
49    ///
50    /// let polygon = polygon![
51    ///     (x: 5., y: 1.),
52    ///     (x: 4., y: 2.),
53    ///     (x: 4., y: 3.),
54    ///     (x: 5., y: 4.),
55    ///     (x: 6., y: 4.),
56    ///     (x: 7., y: 3.),
57    ///     (x: 7., y: 2.),
58    ///     (x: 6., y: 1.),
59    ///     (x: 5., y: 1.),
60    /// ];
61    ///
62    /// let point = point!(x: 2.5, y: 0.5);
63    ///
64    /// let distance = point.euclidean_distance(&polygon);
65    ///
66    /// assert_relative_eq!(distance, 2.1213203435596424);
67    /// ```
68    ///
69    /// `Point` to `LineString`:
70    ///
71    /// ```
72    /// use approx::assert_relative_eq;
73    /// use geo::EuclideanDistance;
74    /// use geo::{point, line_string};
75    ///
76    /// let line_string = line_string![
77    ///     (x: 5., y: 1.),
78    ///     (x: 4., y: 2.),
79    ///     (x: 4., y: 3.),
80    ///     (x: 5., y: 4.),
81    ///     (x: 6., y: 4.),
82    ///     (x: 7., y: 3.),
83    ///     (x: 7., y: 2.),
84    ///     (x: 6., y: 1.),
85    /// ];
86    ///
87    /// let point = point!(x: 5.5, y: 2.1);
88    ///
89    /// let distance = point.euclidean_distance(&line_string);
90    ///
91    /// assert_relative_eq!(distance, 1.1313708498984762);
92    /// ```
93    fn euclidean_distance(&self, rhs: &Rhs) -> T;
94}
95
96// ┌───────────────────────────┐
97// │ Implementations for Coord │
98// └───────────────────────────┘
99
100impl<T> EuclideanDistance<T, Coord<T>> for Coord<T>
101where
102    T: GeoFloat,
103{
104    /// Minimum distance between two `Coord`s
105    fn euclidean_distance(&self, c: &Coord<T>) -> T {
106        Line::new(*self, *c).euclidean_length()
107    }
108}
109
110impl<T> EuclideanDistance<T, Line<T>> for Coord<T>
111where
112    T: GeoFloat,
113{
114    /// Minimum distance from a `Coord` to a `Line`
115    fn euclidean_distance(&self, line: &Line<T>) -> T {
116        line.euclidean_distance(self)
117    }
118}
119
120// ┌───────────────────────────┐
121// │ Implementations for Point │
122// └───────────────────────────┘
123
124impl<T> EuclideanDistance<T, Point<T>> for Point<T>
125where
126    T: GeoFloat,
127{
128    /// Minimum distance between two Points
129    fn euclidean_distance(&self, p: &Point<T>) -> T {
130        self.0.euclidean_distance(&p.0)
131    }
132}
133
134impl<T> EuclideanDistance<T, MultiPoint<T>> for Point<T>
135where
136    T: GeoFloat,
137{
138    /// Minimum distance from a Point to a MultiPoint
139    fn euclidean_distance(&self, points: &MultiPoint<T>) -> T {
140        points
141            .0
142            .iter()
143            .map(|p| self.euclidean_distance(p))
144            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
145    }
146}
147
148impl<T> EuclideanDistance<T, Line<T>> for Point<T>
149where
150    T: GeoFloat,
151{
152    /// Minimum distance from a Line to a Point
153    fn euclidean_distance(&self, line: &Line<T>) -> T {
154        self.0.euclidean_distance(line)
155    }
156}
157
158impl<T> EuclideanDistance<T, LineString<T>> for Point<T>
159where
160    T: GeoFloat,
161{
162    /// Minimum distance from a Point to a LineString
163    fn euclidean_distance(&self, linestring: &LineString<T>) -> T {
164        ::geo_types::private_utils::point_line_string_euclidean_distance(*self, linestring)
165    }
166}
167
168impl<T> EuclideanDistance<T, MultiLineString<T>> for Point<T>
169where
170    T: GeoFloat,
171{
172    /// Minimum distance from a Point to a MultiLineString
173    fn euclidean_distance(&self, mls: &MultiLineString<T>) -> T {
174        mls.0
175            .iter()
176            .map(|ls| self.euclidean_distance(ls))
177            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
178    }
179}
180
181impl<T> EuclideanDistance<T, Polygon<T>> for Point<T>
182where
183    T: GeoFloat,
184{
185    /// Minimum distance from a Point to a Polygon
186    fn euclidean_distance(&self, polygon: &Polygon<T>) -> T {
187        // No need to continue if the polygon intersects the point, or is zero-length
188        if polygon.exterior().0.is_empty() || polygon.intersects(self) {
189            return T::zero();
190        }
191        // fold the minimum interior ring distance if any, followed by the exterior
192        // shell distance, returning the minimum of the two distances
193        polygon
194            .interiors()
195            .iter()
196            .map(|ring| self.euclidean_distance(ring))
197            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
198            .min(
199                polygon
200                    .exterior()
201                    .lines()
202                    .map(|line| {
203                        ::geo_types::private_utils::line_segment_distance(
204                            self.0, line.start, line.end,
205                        )
206                    })
207                    .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val)),
208            )
209    }
210}
211
212impl<T> EuclideanDistance<T, MultiPolygon<T>> for Point<T>
213where
214    T: GeoFloat,
215{
216    /// Minimum distance from a Point to a MultiPolygon
217    fn euclidean_distance(&self, mpolygon: &MultiPolygon<T>) -> T {
218        mpolygon
219            .0
220            .iter()
221            .map(|p| self.euclidean_distance(p))
222            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
223    }
224}
225
226// ┌────────────────────────────────┐
227// │ Implementations for MultiPoint │
228// └────────────────────────────────┘
229
230impl<T> EuclideanDistance<T, Point<T>> for MultiPoint<T>
231where
232    T: GeoFloat,
233{
234    /// Minimum distance from a MultiPoint to a Point
235    fn euclidean_distance(&self, point: &Point<T>) -> T {
236        point.euclidean_distance(self)
237    }
238}
239
240// ┌──────────────────────────┐
241// │ Implementations for Line │
242// └──────────────────────────┘
243
244impl<T> EuclideanDistance<T, Coord<T>> for Line<T>
245where
246    T: GeoFloat,
247{
248    /// Minimum distance from a `Line` to a `Coord`
249    fn euclidean_distance(&self, coord: &Coord<T>) -> T {
250        ::geo_types::private_utils::point_line_euclidean_distance(Point::from(*coord), *self)
251    }
252}
253
254impl<T> EuclideanDistance<T, Point<T>> for Line<T>
255where
256    T: GeoFloat,
257{
258    /// Minimum distance from a Line to a Point
259    fn euclidean_distance(&self, point: &Point<T>) -> T {
260        self.euclidean_distance(&point.0)
261    }
262}
263
264/// Line to Line distance
265impl<T> EuclideanDistance<T, Line<T>> for Line<T>
266where
267    T: GeoFloat + FloatConst + Signed + RTreeNum,
268{
269    fn euclidean_distance(&self, other: &Line<T>) -> T {
270        if self.intersects(other) {
271            return T::zero();
272        }
273        // minimum of all Point-Line distances
274        self.start_point()
275            .euclidean_distance(other)
276            .min(self.end_point().euclidean_distance(other))
277            .min(other.start_point().euclidean_distance(self))
278            .min(other.end_point().euclidean_distance(self))
279    }
280}
281
282/// Line to LineString
283impl<T> EuclideanDistance<T, LineString<T>> for Line<T>
284where
285    T: GeoFloat + FloatConst + Signed + RTreeNum,
286{
287    fn euclidean_distance(&self, other: &LineString<T>) -> T {
288        other.euclidean_distance(self)
289    }
290}
291
292// Line to Polygon distance
293impl<T> EuclideanDistance<T, Polygon<T>> for Line<T>
294where
295    T: GeoFloat + Signed + RTreeNum + FloatConst,
296{
297    fn euclidean_distance(&self, other: &Polygon<T>) -> T {
298        if self.intersects(other) {
299            return T::zero();
300        }
301        // line-line distance between each exterior polygon segment and the line
302        let exterior_min = other
303            .exterior()
304            .lines()
305            .fold(<T as Bounded>::max_value(), |acc, point| {
306                acc.min(self.euclidean_distance(&point))
307            });
308        // line-line distance between each interior ring segment and the line
309        // if there are no rings this just evaluates to max_float
310        let interior_min = other
311            .interiors()
312            .iter()
313            .map(|ring| {
314                ring.lines().fold(<T as Bounded>::max_value(), |acc, line| {
315                    acc.min(self.euclidean_distance(&line))
316                })
317            })
318            .fold(<T as Bounded>::max_value(), |acc, ring_min| {
319                acc.min(ring_min)
320            });
321        // return smaller of the two values
322        exterior_min.min(interior_min)
323    }
324}
325
326/// Line to MultiPolygon distance
327impl<T> EuclideanDistance<T, MultiPolygon<T>> for Line<T>
328where
329    T: GeoFloat + FloatConst + Signed + RTreeNum,
330{
331    fn euclidean_distance(&self, mpolygon: &MultiPolygon<T>) -> T {
332        mpolygon
333            .0
334            .iter()
335            .map(|p| self.euclidean_distance(p))
336            .fold(Bounded::max_value(), |accum, val| accum.min(val))
337    }
338}
339
340// ┌────────────────────────────────┐
341// │ Implementations for LineString │
342// └────────────────────────────────┘
343
344impl<T> EuclideanDistance<T, Point<T>> for LineString<T>
345where
346    T: GeoFloat,
347{
348    /// Minimum distance from a LineString to a Point
349    fn euclidean_distance(&self, point: &Point<T>) -> T {
350        point.euclidean_distance(self)
351    }
352}
353
354/// LineString to Line
355impl<T> EuclideanDistance<T, Line<T>> for LineString<T>
356where
357    T: GeoFloat + FloatConst + Signed + RTreeNum,
358{
359    fn euclidean_distance(&self, other: &Line<T>) -> T {
360        self.lines().fold(Bounded::max_value(), |acc, line| {
361            acc.min(line.euclidean_distance(other))
362        })
363    }
364}
365
366/// LineString-LineString distance
367impl<T> EuclideanDistance<T, LineString<T>> for LineString<T>
368where
369    T: GeoFloat + Signed + RTreeNum,
370{
371    fn euclidean_distance(&self, other: &LineString<T>) -> T {
372        if self.intersects(other) {
373            T::zero()
374        } else {
375            nearest_neighbour_distance(self, other)
376        }
377    }
378}
379
380/// LineString to Polygon
381impl<T> EuclideanDistance<T, Polygon<T>> for LineString<T>
382where
383    T: GeoFloat + FloatConst + Signed + RTreeNum,
384{
385    fn euclidean_distance(&self, other: &Polygon<T>) -> T {
386        if self.intersects(other) {
387            T::zero()
388        } else if !other.interiors().is_empty()
389            && ring_contains_point(other, Point::from(self.0[0]))
390        {
391            // check each ring distance, returning the minimum
392            let mut mindist: T = Float::max_value();
393            for ring in other.interiors() {
394                mindist = mindist.min(nearest_neighbour_distance(self, ring))
395            }
396            mindist
397        } else {
398            nearest_neighbour_distance(self, other.exterior())
399        }
400    }
401}
402
403// ┌─────────────────────────────────────┐
404// │ Implementations for MultiLineString │
405// └─────────────────────────────────────┘
406
407impl<T> EuclideanDistance<T, Point<T>> for MultiLineString<T>
408where
409    T: GeoFloat,
410{
411    /// Minimum distance from a MultiLineString to a Point
412    fn euclidean_distance(&self, point: &Point<T>) -> T {
413        point.euclidean_distance(self)
414    }
415}
416
417// ┌─────────────────────────────┐
418// │ Implementations for Polygon │
419// └─────────────────────────────┘
420
421impl<T> EuclideanDistance<T, Point<T>> for Polygon<T>
422where
423    T: GeoFloat,
424{
425    /// Minimum distance from a Polygon to a Point
426    fn euclidean_distance(&self, point: &Point<T>) -> T {
427        point.euclidean_distance(self)
428    }
429}
430
431// Polygon to Line distance
432impl<T> EuclideanDistance<T, Line<T>> for Polygon<T>
433where
434    T: GeoFloat + FloatConst + Signed + RTreeNum,
435{
436    fn euclidean_distance(&self, other: &Line<T>) -> T {
437        other.euclidean_distance(self)
438    }
439}
440
441/// Polygon to LineString distance
442impl<T> EuclideanDistance<T, LineString<T>> for Polygon<T>
443where
444    T: GeoFloat + FloatConst + Signed + RTreeNum,
445{
446    fn euclidean_distance(&self, other: &LineString<T>) -> T {
447        other.euclidean_distance(self)
448    }
449}
450
451// Polygon to Polygon distance
452impl<T> EuclideanDistance<T, Polygon<T>> for Polygon<T>
453where
454    T: GeoFloat + FloatConst + RTreeNum,
455{
456    fn euclidean_distance(&self, poly2: &Polygon<T>) -> T {
457        if self.intersects(poly2) {
458            return T::zero();
459        }
460        // Containment check
461        if !self.interiors().is_empty()
462            && ring_contains_point(self, Point::from(poly2.exterior().0[0]))
463        {
464            // check each ring distance, returning the minimum
465            let mut mindist: T = Float::max_value();
466            for ring in self.interiors() {
467                mindist = mindist.min(nearest_neighbour_distance(poly2.exterior(), ring))
468            }
469            return mindist;
470        } else if !poly2.interiors().is_empty()
471            && ring_contains_point(poly2, Point::from(self.exterior().0[0]))
472        {
473            let mut mindist: T = Float::max_value();
474            for ring in poly2.interiors() {
475                mindist = mindist.min(nearest_neighbour_distance(self.exterior(), ring))
476            }
477            return mindist;
478        }
479        nearest_neighbour_distance(self.exterior(), poly2.exterior())
480    }
481}
482
483// ┌──────────────────────────────────┐
484// │ Implementations for MultiPolygon │
485// └──────────────────────────────────┘
486
487impl<T> EuclideanDistance<T, Point<T>> for MultiPolygon<T>
488where
489    T: GeoFloat,
490{
491    /// Minimum distance from a MultiPolygon to a Point
492    fn euclidean_distance(&self, point: &Point<T>) -> T {
493        point.euclidean_distance(self)
494    }
495}
496
497/// MultiPolygon to Line distance
498impl<T> EuclideanDistance<T, Line<T>> for MultiPolygon<T>
499where
500    T: GeoFloat + FloatConst + Signed + RTreeNum,
501{
502    fn euclidean_distance(&self, other: &Line<T>) -> T {
503        other.euclidean_distance(self)
504    }
505}
506
507// ┌──────────────────────────────┐
508// │ Implementations for Triangle │
509// └──────────────────────────────┘
510
511impl<T> EuclideanDistance<T, Point<T>> for Triangle<T>
512where
513    T: GeoFloat,
514{
515    fn euclidean_distance(&self, point: &Point<T>) -> T {
516        if self.intersects(point) {
517            return T::zero();
518        }
519
520        [(self.0, self.1), (self.1, self.2), (self.2, self.0)]
521            .iter()
522            .map(|edge| ::geo_types::private_utils::line_segment_distance(point.0, edge.0, edge.1))
523            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
524    }
525}
526
527// ┌───────────┐
528// │ Utilities │
529// └───────────┘
530
531/// This method handles a corner case in which a candidate polygon
532/// is disjoint because it's contained in the inner ring
533/// we work around this by checking that Polygons with inner rings don't
534/// contain a point from the candidate Polygon's outer shell in their simple representations
535fn ring_contains_point<T>(poly: &Polygon<T>, p: Point<T>) -> bool
536where
537    T: GeoNum,
538{
539    match coord_pos_relative_to_ring(p.0, poly.exterior()) {
540        CoordPos::Inside => true,
541        CoordPos::OnBoundary | CoordPos::Outside => false,
542    }
543}
544
545/// Uses an R* tree and nearest-neighbour lookups to calculate minimum distances
546// This is somewhat slow and memory-inefficient, but certainly better than quadratic time
547pub fn nearest_neighbour_distance<T>(geom1: &LineString<T>, geom2: &LineString<T>) -> T
548where
549    T: GeoFloat + RTreeNum,
550{
551    let tree_a: RTree<Line<_>> = RTree::bulk_load(geom1.lines().collect::<Vec<_>>());
552    let tree_b: RTree<Line<_>> = RTree::bulk_load(geom2.lines().collect::<Vec<_>>());
553    // Return minimum distance between all geom a points and geom b lines, and all geom b points and geom a lines
554    geom2
555        .points()
556        .fold(<T as Bounded>::max_value(), |acc, point| {
557            let nearest = tree_a.nearest_neighbor(&point).unwrap();
558            acc.min(nearest.euclidean_distance(&point))
559        })
560        .min(geom1.points().fold(Bounded::max_value(), |acc, point| {
561            let nearest = tree_b.nearest_neighbor(&point).unwrap();
562            acc.min(nearest.euclidean_distance(&point))
563        }))
564}
565
566#[cfg(test)]
567mod test {
568    use super::*;
569    use crate::orient::Direction;
570    use crate::{EuclideanDistance, Orient};
571    use crate::{Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon};
572    use geo_types::{coord, polygon, private_utils::line_segment_distance};
573
574    #[test]
575    fn line_segment_distance_test() {
576        let o1 = Point::new(8.0, 0.0);
577        let o2 = Point::new(5.5, 0.0);
578        let o3 = Point::new(5.0, 0.0);
579        let o4 = Point::new(4.5, 1.5);
580
581        let p1 = Point::new(7.2, 2.0);
582        let p2 = Point::new(6.0, 1.0);
583
584        let dist = line_segment_distance(o1, p1, p2);
585        let dist2 = line_segment_distance(o2, p1, p2);
586        let dist3 = line_segment_distance(o3, p1, p2);
587        let dist4 = line_segment_distance(o4, p1, p2);
588        // Results agree with Shapely
589        assert_relative_eq!(dist, 2.0485900789263356);
590        assert_relative_eq!(dist2, 1.118033988749895);
591        assert_relative_eq!(dist3, std::f64::consts::SQRT_2); // workaround clippy::correctness error approx_constant (1.4142135623730951)
592        assert_relative_eq!(dist4, 1.5811388300841898);
593        // Point is on the line
594        let zero_dist = line_segment_distance(p1, p1, p2);
595        assert_relative_eq!(zero_dist, 0.0);
596    }
597    #[test]
598    // Point to Polygon, outside point
599    fn point_polygon_distance_outside_test() {
600        // an octagon
601        let points = vec![
602            (5., 1.),
603            (4., 2.),
604            (4., 3.),
605            (5., 4.),
606            (6., 4.),
607            (7., 3.),
608            (7., 2.),
609            (6., 1.),
610            (5., 1.),
611        ];
612        let ls = LineString::from(points);
613        let poly = Polygon::new(ls, vec![]);
614        // A Random point outside the octagon
615        let p = Point::new(2.5, 0.5);
616        let dist = p.euclidean_distance(&poly);
617        assert_relative_eq!(dist, 2.1213203435596424);
618    }
619    #[test]
620    // Point to Polygon, inside point
621    fn point_polygon_distance_inside_test() {
622        // an octagon
623        let points = vec![
624            (5., 1.),
625            (4., 2.),
626            (4., 3.),
627            (5., 4.),
628            (6., 4.),
629            (7., 3.),
630            (7., 2.),
631            (6., 1.),
632            (5., 1.),
633        ];
634        let ls = LineString::from(points);
635        let poly = Polygon::new(ls, vec![]);
636        // A Random point inside the octagon
637        let p = Point::new(5.5, 2.1);
638        let dist = p.euclidean_distance(&poly);
639        assert_relative_eq!(dist, 0.0);
640    }
641    #[test]
642    // Point to Polygon, on boundary
643    fn point_polygon_distance_boundary_test() {
644        // an octagon
645        let points = vec![
646            (5., 1.),
647            (4., 2.),
648            (4., 3.),
649            (5., 4.),
650            (6., 4.),
651            (7., 3.),
652            (7., 2.),
653            (6., 1.),
654            (5., 1.),
655        ];
656        let ls = LineString::from(points);
657        let poly = Polygon::new(ls, vec![]);
658        // A point on the octagon
659        let p = Point::new(5.0, 1.0);
660        let dist = p.euclidean_distance(&poly);
661        assert_relative_eq!(dist, 0.0);
662    }
663    #[test]
664    // Point to Polygon, on boundary
665    fn point_polygon_boundary_test2() {
666        let exterior = LineString::from(vec![
667            (0., 0.),
668            (0., 0.0004),
669            (0.0004, 0.0004),
670            (0.0004, 0.),
671            (0., 0.),
672        ]);
673
674        let poly = Polygon::new(exterior, vec![]);
675        let bugged_point = Point::new(0.0001, 0.);
676        assert_relative_eq!(poly.euclidean_distance(&bugged_point), 0.);
677    }
678    #[test]
679    // Point to Polygon, empty Polygon
680    fn point_polygon_empty_test() {
681        // an empty Polygon
682        let points = vec![];
683        let ls = LineString::new(points);
684        let poly = Polygon::new(ls, vec![]);
685        // A point on the octagon
686        let p = Point::new(2.5, 0.5);
687        let dist = p.euclidean_distance(&poly);
688        assert_relative_eq!(dist, 0.0);
689    }
690    #[test]
691    // Point to Polygon with an interior ring
692    fn point_polygon_interior_cutout_test() {
693        // an octagon
694        let ext_points = vec![
695            (4., 1.),
696            (5., 2.),
697            (5., 3.),
698            (4., 4.),
699            (3., 4.),
700            (2., 3.),
701            (2., 2.),
702            (3., 1.),
703            (4., 1.),
704        ];
705        // cut out a triangle inside octagon
706        let int_points = vec![(3.5, 3.5), (4.4, 1.5), (2.6, 1.5), (3.5, 3.5)];
707        let ls_ext = LineString::from(ext_points);
708        let ls_int = LineString::from(int_points);
709        let poly = Polygon::new(ls_ext, vec![ls_int]);
710        // A point inside the cutout triangle
711        let p = Point::new(3.5, 2.5);
712        let dist = p.euclidean_distance(&poly);
713
714        // 0.41036467732879783 <-- Shapely
715        assert_relative_eq!(dist, 0.41036467732879767);
716    }
717
718    #[test]
719    fn line_distance_multipolygon_do_not_intersect_test() {
720        // checks that the distance from the multipolygon
721        // is equal to the distance from the closest polygon
722        // taken in isolation, whatever that distance is
723        let ls1 = LineString::from(vec![
724            (0.0, 0.0),
725            (10.0, 0.0),
726            (10.0, 10.0),
727            (5.0, 15.0),
728            (0.0, 10.0),
729            (0.0, 0.0),
730        ]);
731        let ls2 = LineString::from(vec![
732            (0.0, 30.0),
733            (0.0, 25.0),
734            (10.0, 25.0),
735            (10.0, 30.0),
736            (0.0, 30.0),
737        ]);
738        let ls3 = LineString::from(vec![
739            (15.0, 30.0),
740            (15.0, 25.0),
741            (20.0, 25.0),
742            (20.0, 30.0),
743            (15.0, 30.0),
744        ]);
745        let pol1 = Polygon::new(ls1, vec![]);
746        let pol2 = Polygon::new(ls2, vec![]);
747        let pol3 = Polygon::new(ls3, vec![]);
748        let mp = MultiPolygon::new(vec![pol1.clone(), pol2, pol3]);
749        let pnt1 = Point::new(0.0, 15.0);
750        let pnt2 = Point::new(10.0, 20.0);
751        let ln = Line::new(pnt1.0, pnt2.0);
752        let dist_mp_ln = ln.euclidean_distance(&mp);
753        let dist_pol1_ln = ln.euclidean_distance(&pol1);
754        assert_relative_eq!(dist_mp_ln, dist_pol1_ln);
755    }
756
757    #[test]
758    fn point_distance_multipolygon_test() {
759        let ls1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0), (2.0, 0.0), (0.0, 0.0)]);
760        let ls2 = LineString::from(vec![(3.0, 0.0), (4.0, 10.0), (5.0, 0.0), (3.0, 0.0)]);
761        let p1 = Polygon::new(ls1, vec![]);
762        let p2 = Polygon::new(ls2, vec![]);
763        let mp = MultiPolygon::new(vec![p1, p2]);
764        let p = Point::new(50.0, 50.0);
765        assert_relative_eq!(p.euclidean_distance(&mp), 60.959002616512684);
766    }
767    #[test]
768    // Point to LineString
769    fn point_linestring_distance_test() {
770        // like an octagon, but missing the lowest horizontal segment
771        let points = vec![
772            (5., 1.),
773            (4., 2.),
774            (4., 3.),
775            (5., 4.),
776            (6., 4.),
777            (7., 3.),
778            (7., 2.),
779            (6., 1.),
780        ];
781        let ls = LineString::from(points);
782        // A Random point "inside" the LineString
783        let p = Point::new(5.5, 2.1);
784        let dist = p.euclidean_distance(&ls);
785        assert_relative_eq!(dist, 1.1313708498984762);
786    }
787    #[test]
788    // Point to LineString, point lies on the LineString
789    fn point_linestring_contains_test() {
790        // like an octagon, but missing the lowest horizontal segment
791        let points = vec![
792            (5., 1.),
793            (4., 2.),
794            (4., 3.),
795            (5., 4.),
796            (6., 4.),
797            (7., 3.),
798            (7., 2.),
799            (6., 1.),
800        ];
801        let ls = LineString::from(points);
802        // A point which lies on the LineString
803        let p = Point::new(5.0, 4.0);
804        let dist = p.euclidean_distance(&ls);
805        assert_relative_eq!(dist, 0.0);
806    }
807    #[test]
808    // Point to LineString, closed triangle
809    fn point_linestring_triangle_test() {
810        let points = vec![(3.5, 3.5), (4.4, 2.0), (2.6, 2.0), (3.5, 3.5)];
811        let ls = LineString::from(points);
812        let p = Point::new(3.5, 2.5);
813        let dist = p.euclidean_distance(&ls);
814        assert_relative_eq!(dist, 0.5);
815    }
816    #[test]
817    // Point to LineString, empty LineString
818    fn point_linestring_empty_test() {
819        let points = vec![];
820        let ls = LineString::new(points);
821        let p = Point::new(5.0, 4.0);
822        let dist = p.euclidean_distance(&ls);
823        assert_relative_eq!(dist, 0.0);
824    }
825    #[test]
826    fn distance_multilinestring_test() {
827        let v1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0)]);
828        let v2 = LineString::from(vec![(1.0, 10.0), (2.0, 0.0), (3.0, 1.0)]);
829        let mls = MultiLineString::new(vec![v1, v2]);
830        let p = Point::new(50.0, 50.0);
831        assert_relative_eq!(p.euclidean_distance(&mls), 63.25345840347388);
832    }
833    #[test]
834    fn distance1_test() {
835        assert_relative_eq!(
836            Point::new(0., 0.).euclidean_distance(&Point::new(1., 0.)),
837            1.
838        );
839    }
840    #[test]
841    fn distance2_test() {
842        let dist = Point::new(-72.1235, 42.3521).euclidean_distance(&Point::new(72.1260, 70.612));
843        assert_relative_eq!(dist, 146.99163308930207);
844    }
845    #[test]
846    fn distance_multipoint_test() {
847        let v = vec![
848            Point::new(0.0, 10.0),
849            Point::new(1.0, 1.0),
850            Point::new(10.0, 0.0),
851            Point::new(1.0, -1.0),
852            Point::new(0.0, -10.0),
853            Point::new(-1.0, -1.0),
854            Point::new(-10.0, 0.0),
855            Point::new(-1.0, 1.0),
856            Point::new(0.0, 10.0),
857        ];
858        let mp = MultiPoint::new(v);
859        let p = Point::new(50.0, 50.0);
860        assert_relative_eq!(p.euclidean_distance(&mp), 64.03124237432849)
861    }
862    #[test]
863    fn distance_line_test() {
864        let line0 = Line::from([(0., 0.), (5., 0.)]);
865        let p0 = Point::new(2., 3.);
866        let p1 = Point::new(3., 0.);
867        let p2 = Point::new(6., 0.);
868        assert_relative_eq!(line0.euclidean_distance(&p0), 3.);
869        assert_relative_eq!(p0.euclidean_distance(&line0), 3.);
870
871        assert_relative_eq!(line0.euclidean_distance(&p1), 0.);
872        assert_relative_eq!(p1.euclidean_distance(&line0), 0.);
873
874        assert_relative_eq!(line0.euclidean_distance(&p2), 1.);
875        assert_relative_eq!(p2.euclidean_distance(&line0), 1.);
876    }
877    #[test]
878    fn distance_line_line_test() {
879        let line0 = Line::from([(0., 0.), (5., 0.)]);
880        let line1 = Line::from([(2., 1.), (7., 2.)]);
881        assert_relative_eq!(line0.euclidean_distance(&line1), 1.);
882        assert_relative_eq!(line1.euclidean_distance(&line0), 1.);
883    }
884    #[test]
885    // See https://github.com/georust/geo/issues/476
886    fn distance_line_polygon_test() {
887        let line = Line::new(
888            coord! {
889                x: -0.17084137691985102,
890                y: 0.8748085493016657,
891            },
892            coord! {
893                x: -0.17084137691985102,
894                y: 0.09858870312437906,
895            },
896        );
897        let poly: Polygon<f64> = polygon![
898            coord! {
899                x: -0.10781391405721802,
900                y: -0.15433610862574643,
901            },
902            coord! {
903                x: -0.7855276236615211,
904                y: 0.23694208404779793,
905            },
906            coord! {
907                x: -0.7855276236615214,
908                y: -0.5456143012992907,
909            },
910            coord! {
911                x: -0.10781391405721802,
912                y: -0.15433610862574643,
913            },
914        ];
915        assert_eq!(line.euclidean_distance(&poly), 0.18752558079168907);
916    }
917    #[test]
918    // test edge-vertex minimum distance
919    fn test_minimum_polygon_distance() {
920        let points_raw = vec![
921            (126., 232.),
922            (126., 212.),
923            (112., 202.),
924            (97., 204.),
925            (87., 215.),
926            (87., 232.),
927            (100., 246.),
928            (118., 247.),
929        ];
930        let points = points_raw
931            .iter()
932            .map(|e| Point::new(e.0, e.1))
933            .collect::<Vec<_>>();
934        let poly1 = Polygon::new(LineString::from(points), vec![]);
935
936        let points_raw_2 = vec![
937            (188., 231.),
938            (189., 207.),
939            (174., 196.),
940            (164., 196.),
941            (147., 220.),
942            (158., 242.),
943            (177., 242.),
944        ];
945        let points2 = points_raw_2
946            .iter()
947            .map(|e| Point::new(e.0, e.1))
948            .collect::<Vec<_>>();
949        let poly2 = Polygon::new(LineString::from(points2), vec![]);
950        let dist = nearest_neighbour_distance(poly1.exterior(), poly2.exterior());
951        assert_relative_eq!(dist, 21.0);
952    }
953    #[test]
954    // test vertex-vertex minimum distance
955    fn test_minimum_polygon_distance_2() {
956        let points_raw = vec![
957            (118., 200.),
958            (153., 179.),
959            (106., 155.),
960            (88., 190.),
961            (118., 200.),
962        ];
963        let points = points_raw
964            .iter()
965            .map(|e| Point::new(e.0, e.1))
966            .collect::<Vec<_>>();
967        let poly1 = Polygon::new(LineString::from(points), vec![]);
968
969        let points_raw_2 = vec![
970            (242., 186.),
971            (260., 146.),
972            (182., 175.),
973            (216., 193.),
974            (242., 186.),
975        ];
976        let points2 = points_raw_2
977            .iter()
978            .map(|e| Point::new(e.0, e.1))
979            .collect::<Vec<_>>();
980        let poly2 = Polygon::new(LineString::from(points2), vec![]);
981        let dist = nearest_neighbour_distance(poly1.exterior(), poly2.exterior());
982        assert_relative_eq!(dist, 29.274562336608895);
983    }
984    #[test]
985    // test edge-edge minimum distance
986    fn test_minimum_polygon_distance_3() {
987        let points_raw = vec![
988            (182., 182.),
989            (182., 168.),
990            (138., 160.),
991            (136., 193.),
992            (182., 182.),
993        ];
994        let points = points_raw
995            .iter()
996            .map(|e| Point::new(e.0, e.1))
997            .collect::<Vec<_>>();
998        let poly1 = Polygon::new(LineString::from(points), vec![]);
999
1000        let points_raw_2 = vec![
1001            (232., 196.),
1002            (234., 150.),
1003            (194., 165.),
1004            (194., 191.),
1005            (232., 196.),
1006        ];
1007        let points2 = points_raw_2
1008            .iter()
1009            .map(|e| Point::new(e.0, e.1))
1010            .collect::<Vec<_>>();
1011        let poly2 = Polygon::new(LineString::from(points2), vec![]);
1012        let dist = nearest_neighbour_distance(poly1.exterior(), poly2.exterior());
1013        assert_relative_eq!(dist, 12.0);
1014    }
1015    #[test]
1016    fn test_large_polygon_distance() {
1017        let ls = geo_test_fixtures::norway_main::<f64>();
1018        let poly1 = Polygon::new(ls, vec![]);
1019        let vec2 = vec![
1020            (4.921875, 66.33750501996518),
1021            (3.69140625, 65.21989393613207),
1022            (6.15234375, 65.07213008560697),
1023            (4.921875, 66.33750501996518),
1024        ];
1025        let poly2 = Polygon::new(vec2.into(), vec![]);
1026        let distance = poly1.euclidean_distance(&poly2);
1027        // GEOS says 2.2864896295566055
1028        assert_relative_eq!(distance, 2.2864896295566055);
1029    }
1030    #[test]
1031    // A polygon inside another polygon's ring; they're disjoint in the DE-9IM sense:
1032    // FF2FF1212
1033    fn test_poly_in_ring() {
1034        let shell = geo_test_fixtures::shell::<f64>();
1035        let ring = geo_test_fixtures::ring::<f64>();
1036        let poly_in_ring = geo_test_fixtures::poly_in_ring::<f64>();
1037        // inside is "inside" outside's ring, but they are disjoint
1038        let outside = Polygon::new(shell, vec![ring]);
1039        let inside = Polygon::new(poly_in_ring, vec![]);
1040        assert_relative_eq!(outside.euclidean_distance(&inside), 5.992772737231033);
1041    }
1042    #[test]
1043    // two ring LineStrings; one encloses the other but they neither touch nor intersect
1044    fn test_linestring_distance() {
1045        let ring = geo_test_fixtures::ring::<f64>();
1046        let poly_in_ring = geo_test_fixtures::poly_in_ring::<f64>();
1047        assert_relative_eq!(ring.euclidean_distance(&poly_in_ring), 5.992772737231033);
1048    }
1049    #[test]
1050    // Line-Polygon test: closest point on Polygon is NOT nearest to a Line end-point
1051    fn test_line_polygon_simple() {
1052        let line = Line::from([(0.0, 0.0), (0.0, 3.0)]);
1053        let v = vec![(5.0, 1.0), (5.0, 2.0), (0.25, 1.5), (5.0, 1.0)];
1054        let poly = Polygon::new(v.into(), vec![]);
1055        assert_relative_eq!(line.euclidean_distance(&poly), 0.25);
1056    }
1057    #[test]
1058    // Line-Polygon test: Line intersects Polygon
1059    fn test_line_polygon_intersects() {
1060        let line = Line::from([(0.5, 0.0), (0.0, 3.0)]);
1061        let v = vec![(5.0, 1.0), (5.0, 2.0), (0.25, 1.5), (5.0, 1.0)];
1062        let poly = Polygon::new(v.into(), vec![]);
1063        assert_relative_eq!(line.euclidean_distance(&poly), 0.0);
1064    }
1065    #[test]
1066    // Line-Polygon test: Line contained by interior ring
1067    fn test_line_polygon_inside_ring() {
1068        let line = Line::from([(4.4, 1.5), (4.45, 1.5)]);
1069        let v = vec![(5.0, 1.0), (5.0, 2.0), (0.25, 1.0), (5.0, 1.0)];
1070        let v2 = vec![(4.5, 1.2), (4.5, 1.8), (3.5, 1.2), (4.5, 1.2)];
1071        let poly = Polygon::new(v.into(), vec![v2.into()]);
1072        assert_relative_eq!(line.euclidean_distance(&poly), 0.04999999999999982);
1073    }
1074    #[test]
1075    // LineString-Line test
1076    fn test_linestring_line_distance() {
1077        let line = Line::from([(0.0, 0.0), (0.0, 2.0)]);
1078        let ls: LineString<_> = vec![(3.0, 0.0), (1.0, 1.0), (3.0, 2.0)].into();
1079        assert_relative_eq!(ls.euclidean_distance(&line), 1.0);
1080    }
1081
1082    #[test]
1083    // Triangle-Point test: point on vertex
1084    fn test_triangle_point_on_vertex_distance() {
1085        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
1086        let point = Point::new(0.0, 0.0);
1087        assert_relative_eq!(triangle.euclidean_distance(&point), 0.0);
1088    }
1089
1090    #[test]
1091    // Triangle-Point test: point on edge
1092    fn test_triangle_point_on_edge_distance() {
1093        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
1094        let point = Point::new(1.5, 0.0);
1095        assert_relative_eq!(triangle.euclidean_distance(&point), 0.0);
1096    }
1097
1098    #[test]
1099    // Triangle-Point test
1100    fn test_triangle_point_distance() {
1101        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
1102        let point = Point::new(2.0, 3.0);
1103        assert_relative_eq!(triangle.euclidean_distance(&point), 1.0);
1104    }
1105
1106    #[test]
1107    // Triangle-Point test: point within triangle
1108    fn test_triangle_point_inside_distance() {
1109        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
1110        let point = Point::new(1.0, 0.5);
1111        assert_relative_eq!(triangle.euclidean_distance(&point), 0.0);
1112    }
1113
1114    #[test]
1115    fn convex_and_nearest_neighbour_comparison() {
1116        let ls1: LineString<f64> = vec![
1117            Coord::from((57.39453770777941, 307.60533608924663)),
1118            Coord::from((67.1800355576469, 309.6654408997451)),
1119            Coord::from((84.89693692793338, 225.5101593908847)),
1120            Coord::from((75.1114390780659, 223.45005458038628)),
1121            Coord::from((57.39453770777941, 307.60533608924663)),
1122        ]
1123        .into();
1124        let first_polygon: Polygon<f64> = Polygon::new(ls1, vec![]);
1125        let ls2: LineString<f64> = vec![
1126            Coord::from((138.11769866645008, -45.75134112915392)),
1127            Coord::from((130.50230476949187, -39.270154833870336)),
1128            Coord::from((184.94426964987397, 24.699153900578573)),
1129            Coord::from((192.55966354683218, 18.217967605294987)),
1130            Coord::from((138.11769866645008, -45.75134112915392)),
1131        ]
1132        .into();
1133        let second_polygon = Polygon::new(ls2, vec![]);
1134
1135        assert_relative_eq!(
1136            first_polygon.euclidean_distance(&second_polygon),
1137            224.35357967013238
1138        );
1139    }
1140    #[test]
1141    fn fast_path_regression() {
1142        // this test will fail if the fast path algorithm is reintroduced without being fixed
1143        let p1 = polygon!(
1144            (x: 0_f64, y: 0_f64),
1145            (x: 300_f64, y: 0_f64),
1146            (x: 300_f64, y: 100_f64),
1147            (x: 0_f64, y: 100_f64),
1148        )
1149        .orient(Direction::Default);
1150        let p2 = polygon!(
1151            (x: 100_f64, y: 150_f64),
1152            (x: 150_f64, y: 200_f64),
1153            (x: 50_f64, y: 200_f64),
1154        )
1155        .orient(Direction::Default);
1156        let p3 = polygon!(
1157            (x: 0_f64, y: 0_f64),
1158            (x: 300_f64, y: 0_f64),
1159            (x: 300_f64, y: 100_f64),
1160            (x: 0_f64, y: 100_f64),
1161        )
1162        .orient(Direction::Reversed);
1163        let p4 = polygon!(
1164            (x: 100_f64, y: 150_f64),
1165            (x: 150_f64, y: 200_f64),
1166            (x: 50_f64, y: 200_f64),
1167        )
1168        .orient(Direction::Reversed);
1169        assert_eq!(p1.euclidean_distance(&p2), 50.0f64);
1170        assert_eq!(p3.euclidean_distance(&p4), 50.0f64);
1171        assert_eq!(p1.euclidean_distance(&p4), 50.0f64);
1172        assert_eq!(p2.euclidean_distance(&p3), 50.0f64);
1173    }
1174}
geo/algorithm/euclidean_distance.rs

geo/algorithm/
euclidean_distance.rs
