From 813d33edc6e2646bb185a4596510d5eb5e6d32bd Mon Sep 17 00:00:00 2001
From: Michael Kirk <michael.code@endoftheworl.de>
Date: Tue, 22 Oct 2024 11:26:51 -0700
Subject: [PATCH] Move all geometries impls into `impl Distance for Euclidean`

As part of aligning all our line measure traits (see
https://github.com/georust/geo/issues/1181), all the Point x Point
distance calculations had been moved into `line_measures::metric_spaces`.

Euclidean space is the only space that implements distance calculations
on *non* point geometries. This commit moves all those calculations onto
the `line_measures::metric_spaces::Euclidean` so that we can fully
deprecate the existing EuclideanDistance trait.
---
 geo/benches/euclidean_distance.rs             |   10 +-
 geo/src/algorithm/concave_hull.rs             |   16 +-
 geo/src/algorithm/euclidean_distance.rs       |  222 +---
 geo/src/algorithm/frechet_distance.rs         |   12 +-
 geo/src/algorithm/hausdorff_distance.rs       |    6 +-
 geo/src/algorithm/interior_point.rs           |   22 +-
 .../metric_spaces/euclidean/distance.rs       | 1150 +++++++++++++++++
 .../{euclidean.rs => euclidean/mod.rs}        |   37 +-
 geo/src/algorithm/mod.rs                      |    1 +
 geo/src/algorithm/relate/geomgraph/mod.rs     |    2 +-
 geo/src/algorithm/simplify.rs                 |   12 +-
 geo/src/algorithm/triangulate_spade.rs        |   10 +-
 geo/src/types.rs                              |    4 +-
 geo/src/utils.rs                              |    3 -
 14 files changed, 1271 insertions(+), 236 deletions(-)
 create mode 100644 geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs
 rename geo/src/algorithm/line_measures/metric_spaces/{euclidean.rs => euclidean/mod.rs} (63%)

diff --git a/geo/benches/euclidean_distance.rs b/geo/benches/euclidean_distance.rs
index 6a1c85d60a..b4e2e27c31 100644
--- a/geo/benches/euclidean_distance.rs
+++ b/geo/benches/euclidean_distance.rs
@@ -1,5 +1,5 @@
 use criterion::{criterion_group, criterion_main};
-use geo::algorithm::{ConvexHull, EuclideanDistance};
+use geo::algorithm::{ConvexHull, Distance, Euclidean};
 use geo::{polygon, Polygon};
 
 fn criterion_benchmark(c: &mut criterion::Criterion) {
@@ -37,9 +37,7 @@ fn criterion_benchmark(c: &mut criterion::Criterion) {
             (x: -6.064453, y: 68.49604),
         ];
         bencher.iter(|| {
-            criterion::black_box(
-                criterion::black_box(&poly1).euclidean_distance(criterion::black_box(&poly2)),
-            );
+            criterion::black_box(Euclidean::distance(&poly1, &poly2));
         });
     });
 
@@ -81,9 +79,7 @@ fn criterion_benchmark(c: &mut criterion::Criterion) {
             ]
             .convex_hull();
             bencher.iter(|| {
-                criterion::black_box(
-                    criterion::black_box(&poly1).euclidean_distance(criterion::black_box(&poly2)),
-                );
+                criterion::black_box(Euclidean::distance(&poly1, &poly2));
             });
         },
     );
diff --git a/geo/src/algorithm/concave_hull.rs b/geo/src/algorithm/concave_hull.rs
index e4a288adf9..321c9eb4f6 100644
--- a/geo/src/algorithm/concave_hull.rs
+++ b/geo/src/algorithm/concave_hull.rs
@@ -1,8 +1,8 @@
 use crate::convex_hull::qhull;
 use crate::utils::partial_min;
 use crate::{
-    coord, Centroid, Coord, CoordNum, Euclidean, EuclideanDistance, GeoFloat, Length, Line,
-    LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon,
+    coord, Centroid, Coord, CoordNum, Distance, Euclidean, GeoFloat, Length, Line, LineString,
+    MultiLineString, MultiPoint, MultiPolygon, Point, Polygon,
 };
 use rstar::{RTree, RTreeNum};
 use std::collections::VecDeque;
@@ -134,8 +134,8 @@ where
             let closest_point =
                 candidates.fold(Point::new(point.x, point.y), |acc_point, candidate| {
                     let candidate_point = Point::new(candidate.x, candidate.y);
-                    if line.euclidean_distance(&acc_point)
-                        > line.euclidean_distance(&candidate_point)
+                    if Euclidean::distance(&line, &acc_point)
+                        > Euclidean::distance(&line, &candidate_point)
                     {
                         candidate_point
                     } else {
@@ -154,8 +154,8 @@ where
             let closest_edge_option = match peeked_edge {
                 None => None,
                 Some(&edge) => Some(edges_nearby_point.fold(*edge, |acc, candidate| {
-                    if closest_point.euclidean_distance(&acc)
-                        > closest_point.euclidean_distance(candidate)
+                    if Euclidean::distance(&closest_point, &acc)
+                        > Euclidean::distance(&closest_point, candidate)
                     {
                         *candidate
                     } else {
@@ -164,8 +164,8 @@ where
                 })),
             };
             let decision_distance = partial_min(
-                closest_point.euclidean_distance(&line.start_point()),
-                closest_point.euclidean_distance(&line.end_point()),
+                Euclidean::distance(&closest_point, &line.start_point()),
+                Euclidean::distance(&closest_point, &line.end_point()),
             );
             if let Some(closest_edge) = closest_edge_option {
                 let far_enough = edge_length / decision_distance > concavity;
diff --git a/geo/src/algorithm/euclidean_distance.rs b/geo/src/algorithm/euclidean_distance.rs
index b0ee3c5c84..b66b33ef4f 100644
--- a/geo/src/algorithm/euclidean_distance.rs
+++ b/geo/src/algorithm/euclidean_distance.rs
@@ -1,17 +1,19 @@
-use crate::utils::{coord_pos_relative_to_ring, CoordPos};
 use crate::{
-    Coord, GeoFloat, GeoNum, Geometry, GeometryCollection, Line, LineString, MultiLineString,
-    MultiPoint, MultiPolygon, Point, Polygon, Rect, Triangle,
+    Coord, GeoFloat, Geometry, GeometryCollection, Line, LineString, MultiLineString, MultiPoint,
+    MultiPolygon, Point, Polygon, Rect, Triangle,
 };
-use crate::{Distance, Euclidean, Intersects};
-use num_traits::{float::FloatConst, Bounded, Float, Signed};
+use crate::{Distance, Euclidean};
+use num_traits::{float::FloatConst, Bounded, Signed};
 
 use rstar::primitives::CachedEnvelope;
 use rstar::RTree;
 use rstar::RTreeNum;
 
+#[deprecated(
+    since = "0.29.0",
+    note = "Please use the `Euclidean::distance` method from the `Distance` trait instead"
+)]
 /// Returns the distance between two geometries.
-
 pub trait EuclideanDistance<T, Rhs = Self> {
     /// Returns the distance between two geometries
     ///
@@ -96,24 +98,25 @@ pub trait EuclideanDistance<T, Rhs = Self> {
 // ┌───────────────────────────┐
 // │ Implementations for Coord │
 // └───────────────────────────┘
-
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Coord<T>> for Coord<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance between two `Coord`s
     fn euclidean_distance(&self, c: &Coord<T>) -> T {
-        Euclidean::distance((*self).into(), (*c).into())
+        Euclidean::distance(Point(*self), Point(*c))
     }
 }
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Line<T>> for Coord<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a `Coord` to a `Line`
     fn euclidean_distance(&self, line: &Line<T>) -> T {
-        line.euclidean_distance(self)
+        Euclidean::distance(&Point(*self), line)
     }
 }
 
@@ -121,64 +124,47 @@ where
 // │ Implementations for Point │
 // └───────────────────────────┘
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Point<T>> for Point<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance between two Points
     fn euclidean_distance(&self, p: &Point<T>) -> T {
-        self.0.euclidean_distance(&p.0)
+        Euclidean::distance(*self, *p)
     }
 }
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Line<T>> for Point<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a Line to a Point
     fn euclidean_distance(&self, line: &Line<T>) -> T {
-        self.0.euclidean_distance(line)
+        Euclidean::distance(self, line)
     }
 }
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, LineString<T>> for Point<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a Point to a LineString
-    fn euclidean_distance(&self, linestring: &LineString<T>) -> T {
-        ::geo_types::private_utils::point_line_string_euclidean_distance(*self, linestring)
+    fn euclidean_distance(&self, line_string: &LineString<T>) -> T {
+        Euclidean::distance(self, line_string)
     }
 }
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Polygon<T>> for Point<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a Point to a Polygon
     fn euclidean_distance(&self, polygon: &Polygon<T>) -> T {
-        // No need to continue if the polygon intersects the point, or is zero-length
-        if polygon.exterior().0.is_empty() || polygon.intersects(self) {
-            return T::zero();
-        }
-        // fold the minimum interior ring distance if any, followed by the exterior
-        // shell distance, returning the minimum of the two distances
-        polygon
-            .interiors()
-            .iter()
-            .map(|ring| self.euclidean_distance(ring))
-            .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
-            .min(
-                polygon
-                    .exterior()
-                    .lines()
-                    .map(|line| {
-                        ::geo_types::private_utils::line_segment_distance(
-                            self.0, line.start, line.end,
-                        )
-                    })
-                    .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val)),
-            )
+        Euclidean::distance(self, polygon)
     }
 }
 
@@ -186,85 +172,58 @@ where
 // │ Implementations for Line │
 // └──────────────────────────┘
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Coord<T>> for Line<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a `Line` to a `Coord`
     fn euclidean_distance(&self, coord: &Coord<T>) -> T {
-        ::geo_types::private_utils::point_line_euclidean_distance(Point::from(*coord), *self)
+        Euclidean::distance(self, *coord)
     }
 }
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Point<T>> for Line<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a Line to a Point
     fn euclidean_distance(&self, point: &Point<T>) -> T {
-        self.euclidean_distance(&point.0)
+        Euclidean::distance(self, point)
     }
 }
 
+#[allow(deprecated)]
 /// Line to Line distance
 impl<T> EuclideanDistance<T, Line<T>> for Line<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &Line<T>) -> T {
-        if self.intersects(other) {
-            return T::zero();
-        }
-        // minimum of all Point-Line distances
-        self.start_point()
-            .euclidean_distance(other)
-            .min(self.end_point().euclidean_distance(other))
-            .min(other.start_point().euclidean_distance(self))
-            .min(other.end_point().euclidean_distance(self))
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 /// Line to LineString
 impl<T> EuclideanDistance<T, LineString<T>> for Line<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &LineString<T>) -> T {
-        other.euclidean_distance(self)
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 // Line to Polygon distance
 impl<T> EuclideanDistance<T, Polygon<T>> for Line<T>
 where
     T: GeoFloat + Signed + RTreeNum + FloatConst,
 {
     fn euclidean_distance(&self, other: &Polygon<T>) -> T {
-        if self.intersects(other) {
-            return T::zero();
-        }
-        // line-line distance between each exterior polygon segment and the line
-        let exterior_min = other
-            .exterior()
-            .lines()
-            .fold(<T as Bounded>::max_value(), |acc, point| {
-                acc.min(self.euclidean_distance(&point))
-            });
-        // line-line distance between each interior ring segment 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 as Bounded>::max_value(), |acc, line| {
-                    acc.min(self.euclidean_distance(&line))
-                })
-            })
-            .fold(<T as Bounded>::max_value(), |acc, ring_min| {
-                acc.min(ring_min)
-            });
-        // return smaller of the two values
-        exterior_min.min(interior_min)
+        Euclidean::distance(self, other)
     }
 }
 
@@ -272,62 +231,47 @@ where
 // │ Implementations for LineString │
 // └────────────────────────────────┘
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Point<T>> for LineString<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a LineString to a Point
     fn euclidean_distance(&self, point: &Point<T>) -> T {
-        point.euclidean_distance(self)
+        Euclidean::distance(self, point)
     }
 }
 
+#[allow(deprecated)]
 /// LineString to Line
 impl<T> EuclideanDistance<T, Line<T>> for LineString<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &Line<T>) -> T {
-        self.lines().fold(Bounded::max_value(), |acc, line| {
-            acc.min(line.euclidean_distance(other))
-        })
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 /// LineString-LineString distance
 impl<T> EuclideanDistance<T, LineString<T>> for LineString<T>
 where
     T: GeoFloat + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &LineString<T>) -> T {
-        if self.intersects(other) {
-            T::zero()
-        } else {
-            nearest_neighbour_distance(self, other)
-        }
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 /// LineString to Polygon
 impl<T> EuclideanDistance<T, Polygon<T>> for LineString<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &Polygon<T>) -> T {
-        if self.intersects(other) {
-            T::zero()
-        } else if !other.interiors().is_empty()
-            && ring_contains_point(other, Point::from(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())
-        }
+        Euclidean::distance(self, other)
     }
 }
 
@@ -335,65 +279,47 @@ where
 // │ Implementations for Polygon │
 // └─────────────────────────────┘
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T, Point<T>> for Polygon<T>
 where
     T: GeoFloat,
 {
     /// Minimum distance from a Polygon to a Point
     fn euclidean_distance(&self, point: &Point<T>) -> T {
-        point.euclidean_distance(self)
+        Euclidean::distance(self, point)
     }
 }
 
+#[allow(deprecated)]
 // Polygon to Line distance
 impl<T> EuclideanDistance<T, Line<T>> for Polygon<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &Line<T>) -> T {
-        other.euclidean_distance(self)
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 /// Polygon to LineString distance
 impl<T> EuclideanDistance<T, LineString<T>> for Polygon<T>
 where
     T: GeoFloat + FloatConst + Signed + RTreeNum,
 {
     fn euclidean_distance(&self, other: &LineString<T>) -> T {
-        other.euclidean_distance(self)
+        Euclidean::distance(self, other)
     }
 }
 
+#[allow(deprecated)]
 // Polygon to Polygon distance
 impl<T> EuclideanDistance<T, Polygon<T>> for Polygon<T>
 where
     T: GeoFloat + FloatConst + RTreeNum,
 {
     fn euclidean_distance(&self, poly2: &Polygon<T>) -> T {
-        if self.intersects(poly2) {
-            return T::zero();
-        }
-        // Containment check
-        if !self.interiors().is_empty()
-            && ring_contains_point(self, Point::from(poly2.exterior().0[0]))
-        {
-            // check each ring distance, returning the minimum
-            let mut mindist: T = Float::max_value();
-            for ring in self.interiors() {
-                mindist = mindist.min(nearest_neighbour_distance(poly2.exterior(), ring))
-            }
-            return mindist;
-        } else if !poly2.interiors().is_empty()
-            && ring_contains_point(poly2, Point::from(self.exterior().0[0]))
-        {
-            let mut mindist: T = Float::max_value();
-            for ring in poly2.interiors() {
-                mindist = mindist.min(nearest_neighbour_distance(self.exterior(), ring))
-            }
-            return mindist;
-        }
-        nearest_neighbour_distance(self.exterior(), poly2.exterior())
+        Euclidean::distance(self, poly2)
     }
 }
 
@@ -405,12 +331,13 @@ where
 macro_rules! impl_euclidean_distance_for_polygonlike_geometry {
   ($for:ty,  [$($target:ty),*]) => {
       $(
+          #[allow(deprecated)]
           impl<T> EuclideanDistance<T, $target> for $for
           where
               T: GeoFloat + Signed + RTreeNum + FloatConst,
           {
               fn euclidean_distance(&self, other: &$target) -> T {
-                  other.euclidean_distance(&self.to_polygon())
+                  Euclidean::distance(self, other)
               }
           }
       )*
@@ -424,12 +351,13 @@ impl_euclidean_distance_for_polygonlike_geometry!(Rect<T>,      [Point<T>, Multi
 macro_rules! impl_euclidean_distance_to_polygonlike_geometry {
   ($for:ty,  [$($target:ty),*]) => {
       $(
+          #[allow(deprecated)]
           impl<T> EuclideanDistance<T, $target> for $for
           where
               T: GeoFloat + Signed + RTreeNum + FloatConst,
           {
               fn euclidean_distance(&self, other: &$target) -> T {
-                  other.to_polygon().euclidean_distance(self)
+                  Euclidean::distance(self, other)
               }
           }
       )*
@@ -453,15 +381,13 @@ impl_euclidean_distance_to_polygonlike_geometry!(GeometryCollection<T>,  [Rect<T
 macro_rules! impl_euclidean_distance_for_iter_geometry {
   ($for:ty,  [$($target:ty),*]) => {
       $(
+          #[allow(deprecated)]
           impl<T> EuclideanDistance<T, $target> for $for
           where
               T: GeoFloat + FloatConst + RTreeNum,
           {
               fn euclidean_distance(&self, target: &$target) -> T {
-                  self
-                    .iter()
-                    .map(|g| g.euclidean_distance(target))
-                    .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
+                  Euclidean::distance(self, target)
               }
           }
       )*
@@ -478,12 +404,13 @@ impl_euclidean_distance_for_iter_geometry!(GeometryCollection<T>, [Point<T>, Mul
 macro_rules! impl_euclidean_distance_from_iter_geometry {
   ($for:ty,  [$($target:ty),*]) => {
       $(
+        #[allow(deprecated)]
         impl<T> EuclideanDistance<T, $target> for $for
         where
             T: GeoFloat + FloatConst + RTreeNum
         {
           fn euclidean_distance(&self, target: &$target) -> T {
-            target.euclidean_distance(self)
+              Euclidean::distance(self, target)
           }
         }
       )*
@@ -505,23 +432,13 @@ impl_euclidean_distance_from_iter_geometry!(Polygon<T>,       [MultiPoint<T>, Mu
 macro_rules! impl_euclidean_distance_to_geometry_for_specific {
   ([$($for:ty),*]) => {
       $(
+          #[allow(deprecated)]
           impl<T> EuclideanDistance<T, Geometry<T>> for $for
           where
               T: GeoFloat + FloatConst + RTreeNum,
           {
               fn euclidean_distance(&self, geom: &Geometry<T>) -> T {
-                  match geom {
-                      Geometry::Point(p) => self.euclidean_distance(p),
-                      Geometry::Line(l) => self.euclidean_distance(l),
-                      Geometry::LineString(ls) => self.euclidean_distance(ls),
-                      Geometry::Polygon(p) => self.euclidean_distance(p),
-                      Geometry::MultiPoint(mp) => self.euclidean_distance(mp),
-                      Geometry::MultiLineString(mls) => self.euclidean_distance(mls),
-                      Geometry::MultiPolygon(mp) => self.euclidean_distance(mp),
-                      Geometry::GeometryCollection(gc) => self.euclidean_distance(gc),
-                      Geometry::Rect(r) => self.euclidean_distance(r),
-                      Geometry::Triangle(t) => self.euclidean_distance(t),
-                  }
+                Euclidean::distance(self, geom)
               }
           }
       )*
@@ -538,6 +455,7 @@ impl_euclidean_distance_to_geometry_for_specific!([Point<T>, MultiPoint<T>, Line
 macro_rules! impl_euclidean_distance_to_specific_for_geometry {
   ([$($for:ty),*]) => {
       $(
+          #[allow(deprecated)]
           impl<T> EuclideanDistance<T, $for> for Geometry<T>
           where
               T: GeoFloat + FloatConst + RTreeNum
@@ -552,6 +470,7 @@ macro_rules! impl_euclidean_distance_to_specific_for_geometry {
 
 impl_euclidean_distance_to_specific_for_geometry!([Point<T>, MultiPoint<T>, Line<T>, LineString<T>, MultiLineString<T>, Polygon<T>, MultiPolygon<T>, Triangle<T>, Rect<T>, GeometryCollection<T>]);
 
+#[allow(deprecated)]
 impl<T> EuclideanDistance<T> for Geometry<T>
 where
     T: GeoFloat + FloatConst,
@@ -565,20 +484,6 @@ where
 // │ Utilities │
 // └───────────┘
 
-/// This method handles a corner case in which a candidate polygon
-/// is disjoint because it's contained in the inner ring
-/// we work around this by checking that Polygons with inner rings don't
-/// contain a point from the candidate Polygon's outer shell in their simple representations
-fn ring_contains_point<T>(poly: &Polygon<T>, p: Point<T>) -> bool
-where
-    T: GeoNum,
-{
-    match coord_pos_relative_to_ring(p.0, poly.exterior()) {
-        CoordPos::Inside => true,
-        CoordPos::OnBoundary | CoordPos::Outside => false,
-    }
-}
-
 /// 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
 pub fn nearest_neighbour_distance<T>(geom1: &LineString<T>, geom2: &LineString<T>) -> T
@@ -592,19 +497,26 @@ where
         .points()
         .fold(<T as Bounded>::max_value(), |acc, point| {
             let nearest = tree_a.nearest_neighbor(&point).unwrap();
+            #[allow(deprecated)]
             acc.min(nearest.euclidean_distance(&point))
         })
         .min(geom1.points().fold(Bounded::max_value(), |acc, point| {
             let nearest = tree_b.nearest_neighbor(&point).unwrap();
+            #[allow(deprecated)]
             acc.min(nearest.euclidean_distance(&point))
         }))
 }
 
 #[cfg(test)]
 mod test {
+    // These tests have been ported to the new line_measures::euclidean::distance module.
+    // They'll get removed eventually, so lets leave the deprecated tests here until then
+    // as assurance that they old traits are properly delegating to the new implementation.
+    #![allow(deprecated)]
+
     use super::*;
     use crate::orient::Direction;
-    use crate::{EuclideanDistance, Orient};
+    use crate::Orient;
     use crate::{Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon};
     use geo_types::{coord, polygon, private_utils::line_segment_distance};
 
diff --git a/geo/src/algorithm/frechet_distance.rs b/geo/src/algorithm/frechet_distance.rs
index 5e1c90ee27..b8e8e3ea84 100644
--- a/geo/src/algorithm/frechet_distance.rs
+++ b/geo/src/algorithm/frechet_distance.rs
@@ -1,6 +1,6 @@
 use crate::coords_iter::CoordsIter;
-use crate::euclidean_distance::EuclideanDistance;
-use crate::{GeoFloat, LineString, Point};
+use crate::line_measures::{Distance, Euclidean};
+use crate::{GeoFloat, LineString};
 use num_traits::FromPrimitive;
 
 /// Determine the similarity between two `LineStrings` using the [Frechet distance].
@@ -78,7 +78,7 @@ where
 
         for (i, &a) in self.ls_a.coords().enumerate() {
             for (j, &b) in self.ls_b.coords().enumerate() {
-                let dist = Point::from(a).euclidean_distance(&Point::from(b));
+                let dist = Euclidean::distance(a, b);
 
                 self.cache[i * columns_count + j] = match (i, j) {
                     (0, 0) => dist,
@@ -98,16 +98,14 @@ where
 
 #[cfg(test)]
 mod test {
-    use crate::euclidean_distance::EuclideanDistance;
-    use crate::FrechetDistance;
-    use crate::LineString;
+    use super::*;
 
     #[test]
     fn test_single_point_in_linestring() {
         let ls_a = LineString::from(vec![(1., 1.)]);
         let ls_b = LineString::from(vec![(0., 2.)]);
         assert_relative_eq!(
-            (ls_a.clone().into_points())[0].euclidean_distance(&(&ls_b.clone().into_points())[0]),
+            Euclidean::distance(ls_a.0[0], ls_b.0[0]),
             ls_a.frechet_distance(&ls_b)
         );
     }
diff --git a/geo/src/algorithm/hausdorff_distance.rs b/geo/src/algorithm/hausdorff_distance.rs
index bdb0a35d9f..ec987a238e 100644
--- a/geo/src/algorithm/hausdorff_distance.rs
+++ b/geo/src/algorithm/hausdorff_distance.rs
@@ -1,4 +1,4 @@
-use crate::algorithm::EuclideanDistance;
+use crate::algorithm::{Distance, Euclidean};
 use crate::CoordsIter;
 use crate::GeoFloat;
 use geo_types::{Coord, Point};
@@ -35,7 +35,7 @@ where
             .coords_iter()
             .map(|c| {
                 rhs.coords_iter()
-                    .map(|c2| c.euclidean_distance(&c2))
+                    .map(|c2| Euclidean::distance(c, c2))
                     .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
             })
             .fold(<T as Bounded>::min_value(), |accum, val| accum.max(val));
@@ -45,7 +45,7 @@ where
             .coords_iter()
             .map(|c| {
                 self.coords_iter()
-                    .map(|c2| c.euclidean_distance(&c2))
+                    .map(|c2| Euclidean::distance(c, c2))
                     .fold(<T as Bounded>::max_value(), |accum, val| accum.min(val))
             })
             .fold(<T as Bounded>::min_value(), |accum, val| accum.max(val));
diff --git a/geo/src/algorithm/interior_point.rs b/geo/src/algorithm/interior_point.rs
index 7f34c905e7..7695d9c555 100644
--- a/geo/src/algorithm/interior_point.rs
+++ b/geo/src/algorithm/interior_point.rs
@@ -1,9 +1,14 @@
 use std::cmp::{Ordering, Reverse};
 
 use crate::algorithm::{
-    bounding_rect::BoundingRect, centroid::Centroid, coords_iter::CoordsIter,
-    dimensions::HasDimensions, euclidean_distance::EuclideanDistance,
-    line_intersection::LineIntersection, lines_iter::LinesIter, relate::Relate,
+    bounding_rect::BoundingRect,
+    centroid::Centroid,
+    coords_iter::CoordsIter,
+    dimensions::HasDimensions,
+    line_intersection::LineIntersection,
+    line_measures::{Distance, Euclidean},
+    lines_iter::LinesIter,
+    relate::Relate,
 };
 use crate::geometry::*;
 use crate::sweep::{Intersections, SweepPoint};
@@ -107,7 +112,7 @@ where
                     .iter()
                     .map(|coord| {
                         let pt = Point::from(*coord);
-                        (pt, pt.euclidean_distance(&centroid))
+                        (pt, Euclidean::distance(pt, centroid))
                     })
                     .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(Ordering::Less))
                     .map(|(pt, _distance)| pt)
@@ -130,7 +135,7 @@ where
                 .filter_map(|linestring| {
                     linestring
                         .interior_point()
-                        .map(|pt| (pt, pt.euclidean_distance(&centroid)))
+                        .map(|pt| (pt, Euclidean::distance(pt, centroid)))
                 })
                 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(Ordering::Less))
                 .map(|(pt, _distance)| pt)
@@ -308,7 +313,7 @@ where
     fn interior_point(&self) -> Self::Output {
         if let Some(centroid) = self.centroid() {
             self.iter()
-                .map(|pt| (pt, pt.euclidean_distance(&centroid)))
+                .map(|pt| (pt, Euclidean::distance(pt, &centroid)))
                 .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(Ordering::Less))
                 .map(|(pt, _distance)| *pt)
         } else {
@@ -342,7 +347,10 @@ where
                         (
                             pt,
                             // maximize dimensions, minimize distance
-                            (Reverse(geom.dimensions()), pt.euclidean_distance(&centroid)),
+                            (
+                                Reverse(geom.dimensions()),
+                                Euclidean::distance(pt, centroid),
+                            ),
                         )
                     })
                 })
diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs
new file mode 100644
index 0000000000..41a5b89d6a
--- /dev/null
+++ b/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs
@@ -0,0 +1,1150 @@
+use super::{Distance, Euclidean};
+use crate::algorithm::Intersects;
+use crate::coordinate_position::{coord_pos_relative_to_ring, CoordPos};
+use crate::geometry::*;
+use crate::{CoordFloat, GeoFloat, GeoNum};
+use num_traits::{Bounded, Float};
+use rstar::primitives::CachedEnvelope;
+use rstar::RTree;
+
+// Distance is a symmetric operation, so we can implement it once for both
+macro_rules! symmetric_distance_impl {
+    ($t:ident, $a:ty, $b:ty) => {
+        impl<F> $crate::Distance<F, $a, $b> for Euclidean
+        where
+            F: $t,
+        {
+            fn distance(a: $a, b: $b) -> F {
+                Self::distance(b, a)
+            }
+        }
+    };
+}
+
+// ┌───────────────────────────┐
+// │ Implementations for Coord │
+// └───────────────────────────┘
+
+impl<F: CoordFloat> Distance<F, Coord<F>, Coord<F>> for Euclidean {
+    fn distance(origin: Coord<F>, destination: Coord<F>) -> F {
+        let delta = origin - destination;
+        delta.x.hypot(delta.y)
+    }
+}
+impl<F: CoordFloat> Distance<F, Coord<F>, &Line<F>> for Euclidean {
+    fn distance(coord: Coord<F>, line: &Line<F>) -> F {
+        ::geo_types::private_utils::point_line_euclidean_distance(Point(coord), *line)
+    }
+}
+
+// ┌───────────────────────────┐
+// │ Implementations for Point │
+// └───────────────────────────┘
+
+/// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points
+impl<F: CoordFloat> Distance<F, Point<F>, Point<F>> for Euclidean {
+    /// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points
+    ///
+    /// # Units
+    /// - `origin`, `destination`: Point where the units of x/y represent non-angular units
+    ///    — e.g. meters or miles, not lon/lat. For lon/lat points, use the
+    ///    [`Haversine`] or [`Geodesic`] [metric spaces].
+    /// - returns: distance in the same units as the `origin` and `destination` points
+    ///
+    /// # Example
+    /// ```
+    /// use geo::{Euclidean, Distance};
+    /// use geo::Point;
+    /// // web mercator
+    /// let new_york_city = Point::new(-8238310.24, 4942194.78);
+    /// // web mercator
+    /// let london = Point::new(-14226.63, 6678077.70);
+    /// let distance: f64 = Euclidean::distance(new_york_city, london);
+    ///
+    /// assert_eq!(
+    ///     8_405_286., // meters in web mercator
+    ///     distance.round()
+    /// );
+    /// ```
+    ///
+    /// [`Haversine`]: super::Haversine
+    /// [`Geodesic`]: super::Geodesic
+    /// [metric spaces]: super
+    fn distance(origin: Point<F>, destination: Point<F>) -> F {
+        Self::distance(origin.0, destination.0)
+    }
+}
+
+impl<F: CoordFloat> Distance<F, &Point<F>, &Point<F>> for Euclidean {
+    fn distance(origin: &Point<F>, destination: &Point<F>) -> F {
+        Self::distance(*origin, *destination)
+    }
+}
+
+impl<F: CoordFloat> Distance<F, &Point<F>, &Line<F>> for Euclidean {
+    fn distance(origin: &Point<F>, destination: &Line<F>) -> F {
+        geo_types::private_utils::point_line_euclidean_distance(*origin, *destination)
+    }
+}
+
+impl<F: CoordFloat> Distance<F, &Point<F>, &LineString<F>> for Euclidean {
+    fn distance(origin: &Point<F>, destination: &LineString<F>) -> F {
+        geo_types::private_utils::point_line_string_euclidean_distance(*origin, destination)
+    }
+}
+
+impl<F: GeoFloat> Distance<F, &Point<F>, &Polygon<F>> for Euclidean {
+    fn distance(point: &Point<F>, polygon: &Polygon<F>) -> F {
+        // No need to continue if the polygon intersects the point, or is zero-length
+        if polygon.exterior().0.is_empty() || polygon.intersects(point) {
+            return F::zero();
+        }
+        // fold the minimum interior ring distance if any, followed by the exterior
+        // shell distance, returning the minimum of the two distances
+        polygon
+            .interiors()
+            .iter()
+            .map(|ring| Self::distance(point, ring))
+            .fold(Bounded::max_value(), |accum: F, val| accum.min(val))
+            .min(
+                polygon
+                    .exterior()
+                    .lines()
+                    .map(|line| {
+                        ::geo_types::private_utils::line_segment_distance(
+                            point.0, line.start, line.end,
+                        )
+                    })
+                    .fold(Bounded::max_value(), |accum, val| accum.min(val)),
+            )
+    }
+}
+
+// ┌──────────────────────────┐
+// │ Implementations for Line │
+// └──────────────────────────┘
+
+symmetric_distance_impl!(CoordFloat, &Line<F>, Coord<F>);
+symmetric_distance_impl!(CoordFloat, &Line<F>, &Point<F>);
+
+impl<F: GeoFloat> Distance<F, &Line<F>, &Line<F>> for Euclidean {
+    fn distance(line_a: &Line<F>, line_b: &Line<F>) -> F {
+        if line_a.intersects(line_b) {
+            return F::zero();
+        }
+        // minimum of all Point-Line distances
+        Self::distance(&line_a.start_point(), line_b)
+            .min(Self::distance(&line_a.end_point(), line_b))
+            .min(Self::distance(&line_b.start_point(), line_a))
+            .min(Self::distance(&line_b.end_point(), line_a))
+    }
+}
+
+impl<F: GeoFloat> Distance<F, &Line<F>, &LineString<F>> for Euclidean {
+    fn distance(line: &Line<F>, line_string: &LineString<F>) -> F {
+        line_string
+            .lines()
+            .fold(Bounded::max_value(), |acc, segment| {
+                acc.min(Self::distance(line, &segment))
+            })
+    }
+}
+
+impl<F: GeoFloat> Distance<F, &Line<F>, &Polygon<F>> for Euclidean {
+    fn distance(line: &Line<F>, polygon: &Polygon<F>) -> F {
+        if line.intersects(polygon) {
+            return F::zero();
+        }
+
+        // REVIEW: This impl changed slightly.
+        std::iter::once(polygon.exterior())
+            .chain(polygon.interiors().iter())
+            .fold(Bounded::max_value(), |acc, line_string| {
+                acc.min(Self::distance(line, line_string))
+            })
+    }
+}
+
+// ┌────────────────────────────────┐
+// │ Implementations for LineString │
+// └────────────────────────────────┘
+
+symmetric_distance_impl!(CoordFloat, &LineString<F>, &Point<F>);
+symmetric_distance_impl!(GeoFloat, &LineString<F>, &Line<F>);
+
+impl<F: GeoFloat> Distance<F, &LineString<F>, &LineString<F>> for Euclidean {
+    fn distance(line_string_a: &LineString<F>, line_string_b: &LineString<F>) -> F {
+        if line_string_a.intersects(line_string_b) {
+            F::zero()
+        } else {
+            nearest_neighbour_distance(line_string_a, line_string_b)
+        }
+    }
+}
+
+impl<F: GeoFloat> Distance<F, &LineString<F>, &Polygon<F>> for Euclidean {
+    fn distance(line_string: &LineString<F>, polygon: &Polygon<F>) -> F {
+        if line_string.intersects(polygon) {
+            F::zero()
+        } else if !polygon.interiors().is_empty()
+            // FIXME: Explodes on empty line_string
+            && ring_contains_coord(polygon.exterior(), line_string.0[0])
+        {
+            // check each ring distance, returning the minimum
+            let mut mindist: F = Float::max_value();
+            for ring in polygon.interiors() {
+                mindist = mindist.min(crate::euclidean_distance::nearest_neighbour_distance(
+                    line_string,
+                    ring,
+                ))
+            }
+            mindist
+        } else {
+            crate::euclidean_distance::nearest_neighbour_distance(line_string, polygon.exterior())
+        }
+    }
+}
+
+// ┌─────────────────────────────┐
+// │ Implementations for Polygon │
+// └─────────────────────────────┘
+
+symmetric_distance_impl!(GeoFloat, &Polygon<F>, &Point<F>);
+symmetric_distance_impl!(GeoFloat, &Polygon<F>, &Line<F>);
+symmetric_distance_impl!(GeoFloat, &Polygon<F>, &LineString<F>);
+impl<F: GeoFloat> Distance<F, &Polygon<F>, &Polygon<F>> for Euclidean {
+    fn distance(polygon_a: &Polygon<F>, polygon_b: &Polygon<F>) -> F {
+        if polygon_a.intersects(polygon_b) {
+            return F::zero();
+        }
+        // FIXME: explodes when polygon_b.exterior() is empty
+        // Containment check
+        if !polygon_a.interiors().is_empty()
+            && ring_contains_coord(polygon_a.exterior(), polygon_b.exterior().0[0])
+        {
+            // check each ring distance, returning the minimum
+            let mut mindist: F = Float::max_value();
+            for ring in polygon_a.interiors() {
+                mindist = mindist.min(nearest_neighbour_distance(polygon_b.exterior(), ring))
+            }
+            return mindist;
+        } else if !polygon_b.interiors().is_empty()
+            // FIXME: explodes when polygon_a.exterior() is empty
+            && ring_contains_coord(polygon_b.exterior(), polygon_a.exterior().0[0])
+        {
+            let mut mindist: F = Float::max_value();
+            for ring in polygon_b.interiors() {
+                mindist = mindist.min(nearest_neighbour_distance(polygon_a.exterior(), ring))
+            }
+            return mindist;
+        }
+        nearest_neighbour_distance(polygon_a.exterior(), polygon_b.exterior())
+    }
+}
+
+// ┌────────────────────────────────────────┐
+// │ Implementations for Rect and Triangle  │
+// └────────────────────────────────────────┘
+
+/// Implements Euclidean distance for Triangles and Rects by converting them to polygons.
+macro_rules! impl_euclidean_distance_from_polygonlike_geometry {
+  ($polygonlike:ty,  [$($geometry_b:ty),*]) => {
+      $(
+          impl<F: GeoFloat> Distance<F, $polygonlike, $geometry_b> for Euclidean
+          {
+              fn distance(polygonlike: $polygonlike, geometry_b: $geometry_b) -> F {
+                    Self::distance(&polygonlike.to_polygon(), geometry_b)
+              }
+          }
+      )*
+  };
+}
+
+impl_euclidean_distance_from_polygonlike_geometry!(&Triangle<F>,  [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>, &Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_from_polygonlike_geometry!(&Rect<F>,      [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>, &Rect<F>, &Triangle<F>]);
+
+/// Implements Euclidean distance for other geometry types to Triangles and Rects by converting the Triangle or Rect to a polygon.
+macro_rules! impl_euclidean_distance_to_polygonlike_geometry {
+  ($geometry_a:ty,  [$($polygonlike:ty),*]) => {
+      $(
+          impl<F: GeoFloat> Distance<F, $geometry_a, $polygonlike> for Euclidean
+          {
+              fn distance(geometry_a: $geometry_a, polygonlike: $polygonlike) -> F {
+                    Self::distance(geometry_a, &polygonlike.to_polygon())
+              }
+          }
+      )*
+  };
+}
+
+impl_euclidean_distance_to_polygonlike_geometry!(&Point<F>,               [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&MultiPoint<F>,          [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&Line<F>,                [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&LineString<F>,          [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&MultiLineString<F>,     [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&Polygon<F>,             [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&MultiPolygon<F>,        [&Rect<F>, &Triangle<F>]);
+impl_euclidean_distance_to_polygonlike_geometry!(&GeometryCollection<F>,  [&Rect<F>, &Triangle<F>]);
+
+// ┌───────────────────────────────────────────┐
+// │ Implementations for multi geometry types  │
+// └───────────────────────────────────────────┘
+
+/// Euclidean distance implementation for multi geometry types.
+macro_rules! impl_euclidean_distance_from_iter_geometry {
+  ($iter_geometry:ty,  [$($to_geometry:ty),*]) => {
+      $(
+          impl<F: GeoFloat> Distance<F, $iter_geometry, $to_geometry> for Euclidean
+          {
+              fn distance(iter_geometry: $iter_geometry, to_geometry: $to_geometry) -> F {
+                  iter_geometry
+                    .iter()
+                    .fold(Bounded::max_value(), |accum: F, member| {
+                      accum.min(Self::distance(member, to_geometry))
+                    })
+              }
+          }
+      )*
+  };
+}
+
+impl_euclidean_distance_from_iter_geometry!(&MultiPoint<F>,         [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_from_iter_geometry!(&MultiLineString<F>,    [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_from_iter_geometry!(&MultiPolygon<F>,       [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_from_iter_geometry!(&GeometryCollection<F>, [&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+
+/// Euclidean distance implementation for other geometry types to multi geometry types,
+/// using the multi geometry type's implementation.
+macro_rules! impl_euclidean_distance_to_iter_geometry {
+  ($from_geometry:ty,  [$($iter_geometry:ty),*]) => {
+      $(
+        // Can I inline this in  impl_euclidean_distance_from_iter_geometry and get rid of impl_euclidean_distance_to_iter_geometry? Except... if $to == $from?
+        symmetric_distance_impl!(GeoFloat, $from_geometry, $iter_geometry);
+      )*
+  };
+}
+
+// This macro is used to implement EuclideanDistance to multi geometry types for non-multi geometry types.
+// Rect and Triangle are omitted here because those implementations are included in the Rect and Triangle section above.
+impl_euclidean_distance_to_iter_geometry!(&Point<F>,         [&MultiPoint<F>, &MultiLineString<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_to_iter_geometry!(&Line<F>,          [&MultiPoint<F>, &MultiLineString<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_to_iter_geometry!(&LineString<F>,    [&MultiPoint<F>, &MultiLineString<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+impl_euclidean_distance_to_iter_geometry!(&Polygon<F>,       [&MultiPoint<F>, &MultiLineString<F>, &MultiPolygon<F>, &GeometryCollection<F>]);
+
+// ┌─────────────────────────────────────────────────────────┐
+// │ Implementation to Geometry<T> for every geometry type   │
+// └─────────────────────────────────────────────────────────┘
+
+/// Euclidean distance implementation for every specific Geometry type to Geometry<T>.
+macro_rules! impl_euclidean_distance_to_geometry_for_specific {
+  ([$($from_geometry:ty),*]) => {
+      $(
+          impl<F: GeoFloat> Distance<F, $from_geometry, &Geometry<F>> for Euclidean
+          {
+              fn distance(from_geometry: $from_geometry, geom: &Geometry<F>) -> F {
+                  match geom {
+                      Geometry::Point(p) => Self::distance(from_geometry, p),
+                      Geometry::Line(l) => Self::distance(from_geometry, l),
+                      Geometry::LineString(ls) => Self::distance(from_geometry, ls),
+                      Geometry::Polygon(p) => Self::distance(from_geometry, p),
+                      Geometry::MultiPoint(mp) => Self::distance(from_geometry, mp),
+                      Geometry::MultiLineString(mls) => Self::distance(from_geometry, mls),
+                      Geometry::MultiPolygon(mp) => Self::distance(from_geometry, mp),
+                      Geometry::GeometryCollection(gc) => Self::distance(from_geometry, gc),
+                      Geometry::Rect(r) => Self::distance(from_geometry, r),
+                      Geometry::Triangle(t) => Self::distance(from_geometry, t),
+                  }
+              }
+          }
+      )*
+  };
+}
+
+impl_euclidean_distance_to_geometry_for_specific!([&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &Triangle<F>, &Rect<F>, &GeometryCollection<F>]);
+// ┌──────────────────────────────┐
+// │ Implementation for Geometry  │
+// └──────────────────────────────┘
+
+/// Euclidean distance implementation for Geometry<T> to every specific Geometry type.
+macro_rules! impl_euclidean_distance_to_specific_for_geometry {
+    ([$($destination:ty),*]) => {
+        $(
+            impl<F: GeoFloat> Distance<F, &Geometry<F>, $destination> for Euclidean {
+                fn distance(origin: &Geometry<F>, destination: $destination) -> F {
+                    match origin {
+                        Geometry::Point(point) => Self::distance(point, destination),
+                        Geometry::Line(line) => Self::distance(line, destination),
+                        Geometry::LineString(line_string) => Self::distance(line_string, destination),
+                        Geometry::Polygon(polygon) => Self::distance(polygon, destination),
+                        Geometry::MultiPoint(multi_point) => Self::distance(multi_point, destination),
+                        Geometry::MultiLineString(multi_line_string) => Self::distance(multi_line_string, destination),
+                        Geometry::MultiPolygon(multi_polygon) => Self::distance(multi_polygon, destination),
+                        Geometry::GeometryCollection(geometry_collection) => Self::distance(geometry_collection, destination),
+                        Geometry::Rect(rect) => Self::distance(rect, destination),
+                        Geometry::Triangle(triangle) => Self::distance(triangle, destination),
+                    }
+                }
+            }
+        )*
+    };
+}
+
+impl_euclidean_distance_to_specific_for_geometry!([&Point<F>, &MultiPoint<F>, &Line<F>, &LineString<F>, &MultiLineString<F>, &Polygon<F>, &MultiPolygon<F>, &Triangle<F>, &Rect<F>, &GeometryCollection<F>]);
+
+impl<F: GeoFloat> Distance<F, &Geometry<F>, &Geometry<F>> for Euclidean {
+    fn distance(origin: &Geometry<F>, destination: &Geometry<F>) -> F {
+        match origin {
+            Geometry::Point(point) => Self::distance(point, destination),
+            Geometry::Line(line) => Self::distance(line, destination),
+            Geometry::LineString(line_string) => Self::distance(line_string, destination),
+            Geometry::Polygon(polygon) => Self::distance(polygon, destination),
+            Geometry::MultiPoint(multi_point) => Self::distance(multi_point, destination),
+            Geometry::MultiLineString(multi_line_string) => {
+                Self::distance(multi_line_string, destination)
+            }
+            Geometry::MultiPolygon(multi_polygon) => Self::distance(multi_polygon, destination),
+            Geometry::GeometryCollection(geometry_collection) => {
+                Self::distance(geometry_collection, destination)
+            }
+            Geometry::Rect(rect) => Self::distance(rect, destination),
+            Geometry::Triangle(triangle) => Self::distance(triangle, destination),
+        }
+    }
+}
+
+// ┌───────────────────────────┐
+// │ Implementations utilities │
+// └───────────────────────────┘
+
+/// 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<F: GeoFloat>(geom1: &LineString<F>, geom2: &LineString<F>) -> F {
+    let tree_a = RTree::bulk_load(geom1.lines().map(CachedEnvelope::new).collect());
+    let tree_b = RTree::bulk_load(geom2.lines().map(CachedEnvelope::new).collect());
+    // Return minimum distance between all geom a points and geom b lines, and all geom b points and geom a lines
+    geom2
+        .points()
+        .fold(Bounded::max_value(), |acc: F, point| {
+            let nearest = tree_a.nearest_neighbor(&point).unwrap();
+            acc.min(Euclidean::distance(nearest as &Line<F>, &point))
+        })
+        .min(geom1.points().fold(Bounded::max_value(), |acc, point| {
+            let nearest = tree_b.nearest_neighbor(&point).unwrap();
+            acc.min(Euclidean::distance(nearest as &Line<F>, &point))
+        }))
+}
+
+fn ring_contains_coord<T: GeoNum>(ring: &LineString<T>, c: Coord<T>) -> bool {
+    match coord_pos_relative_to_ring(c, ring) {
+        CoordPos::Inside => true,
+        CoordPos::OnBoundary | CoordPos::Outside => false,
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    use crate::orient::{Direction, Orient};
+    use crate::{Line, LineString, MultiLineString, MultiPoint, MultiPolygon, Point, Polygon};
+    use geo_types::{coord, polygon, private_utils::line_segment_distance};
+
+    #[test]
+    fn line_segment_distance_test() {
+        let o1 = Point::new(8.0, 0.0);
+        let o2 = Point::new(5.5, 0.0);
+        let o3 = Point::new(5.0, 0.0);
+        let o4 = Point::new(4.5, 1.5);
+
+        let p1 = Point::new(7.2, 2.0);
+        let p2 = Point::new(6.0, 1.0);
+
+        let dist = line_segment_distance(o1, p1, p2);
+        let dist2 = line_segment_distance(o2, p1, p2);
+        let dist3 = line_segment_distance(o3, p1, p2);
+        let dist4 = line_segment_distance(o4, p1, p2);
+        // Results agree with Shapely
+        assert_relative_eq!(dist, 2.0485900789263356);
+        assert_relative_eq!(dist2, 1.118033988749895);
+        assert_relative_eq!(dist3, std::f64::consts::SQRT_2); // workaround clippy::correctness error approx_constant (1.4142135623730951)
+        assert_relative_eq!(dist4, 1.5811388300841898);
+        // Point is on the line
+        let zero_dist = line_segment_distance(p1, p1, p2);
+        assert_relative_eq!(zero_dist, 0.0);
+    }
+    #[test]
+    // Point to Polygon, outside point
+    fn point_polygon_distance_outside_test() {
+        // an octagon
+        let points = vec![
+            (5., 1.),
+            (4., 2.),
+            (4., 3.),
+            (5., 4.),
+            (6., 4.),
+            (7., 3.),
+            (7., 2.),
+            (6., 1.),
+            (5., 1.),
+        ];
+        let ls = LineString::from(points);
+        let poly = Polygon::new(ls, vec![]);
+        // A Random point outside the octagon
+        let p = Point::new(2.5, 0.5);
+        let dist = Euclidean::distance(&p, &poly);
+        assert_relative_eq!(dist, 2.1213203435596424);
+    }
+    #[test]
+    // Point to Polygon, inside point
+    fn point_polygon_distance_inside_test() {
+        // an octagon
+        let points = vec![
+            (5., 1.),
+            (4., 2.),
+            (4., 3.),
+            (5., 4.),
+            (6., 4.),
+            (7., 3.),
+            (7., 2.),
+            (6., 1.),
+            (5., 1.),
+        ];
+        let ls = LineString::from(points);
+        let poly = Polygon::new(ls, vec![]);
+        // A Random point inside the octagon
+        let p = Point::new(5.5, 2.1);
+        let dist = Euclidean::distance(&p, &poly);
+        assert_relative_eq!(dist, 0.0);
+    }
+    #[test]
+    // Point to Polygon, on boundary
+    fn point_polygon_distance_boundary_test() {
+        // an octagon
+        let points = vec![
+            (5., 1.),
+            (4., 2.),
+            (4., 3.),
+            (5., 4.),
+            (6., 4.),
+            (7., 3.),
+            (7., 2.),
+            (6., 1.),
+            (5., 1.),
+        ];
+        let ls = LineString::from(points);
+        let poly = Polygon::new(ls, vec![]);
+        // A point on the octagon
+        let p = Point::new(5.0, 1.0);
+        let dist = Euclidean::distance(&p, &poly);
+        assert_relative_eq!(dist, 0.0);
+    }
+    #[test]
+    // Point to Polygon, on boundary
+    fn point_polygon_boundary_test2() {
+        let exterior = LineString::from(vec![
+            (0., 0.),
+            (0., 0.0004),
+            (0.0004, 0.0004),
+            (0.0004, 0.),
+            (0., 0.),
+        ]);
+
+        let poly = Polygon::new(exterior, vec![]);
+        let bugged_point = Point::new(0.0001, 0.);
+        assert_relative_eq!(Euclidean::distance(&poly, &bugged_point), 0.);
+    }
+    #[test]
+    // Point to Polygon, empty Polygon
+    fn point_polygon_empty_test() {
+        // an empty Polygon
+        let points = vec![];
+        let ls = LineString::new(points);
+        let poly = Polygon::new(ls, vec![]);
+        // A point on the octagon
+        let p = Point::new(2.5, 0.5);
+        let dist = Euclidean::distance(&p, &poly);
+        assert_relative_eq!(dist, 0.0);
+    }
+    #[test]
+    // Point to Polygon with an interior ring
+    fn point_polygon_interior_cutout_test() {
+        // an octagon
+        let ext_points = vec![
+            (4., 1.),
+            (5., 2.),
+            (5., 3.),
+            (4., 4.),
+            (3., 4.),
+            (2., 3.),
+            (2., 2.),
+            (3., 1.),
+            (4., 1.),
+        ];
+        // cut out a triangle inside octagon
+        let int_points = vec![(3.5, 3.5), (4.4, 1.5), (2.6, 1.5), (3.5, 3.5)];
+        let ls_ext = LineString::from(ext_points);
+        let ls_int = LineString::from(int_points);
+        let poly = Polygon::new(ls_ext, vec![ls_int]);
+        // A point inside the cutout triangle
+        let p = Point::new(3.5, 2.5);
+        let dist = Euclidean::distance(&p, &poly);
+
+        // 0.41036467732879783 <-- Shapely
+        assert_relative_eq!(dist, 0.41036467732879767);
+    }
+
+    #[test]
+    fn line_distance_multipolygon_do_not_intersect_test() {
+        // checks that the distance from the multipolygon
+        // is equal to the distance from the closest polygon
+        // taken in isolation, whatever that distance is
+        let ls1 = LineString::from(vec![
+            (0.0, 0.0),
+            (10.0, 0.0),
+            (10.0, 10.0),
+            (5.0, 15.0),
+            (0.0, 10.0),
+            (0.0, 0.0),
+        ]);
+        let ls2 = LineString::from(vec![
+            (0.0, 30.0),
+            (0.0, 25.0),
+            (10.0, 25.0),
+            (10.0, 30.0),
+            (0.0, 30.0),
+        ]);
+        let ls3 = LineString::from(vec![
+            (15.0, 30.0),
+            (15.0, 25.0),
+            (20.0, 25.0),
+            (20.0, 30.0),
+            (15.0, 30.0),
+        ]);
+        let pol1 = Polygon::new(ls1, vec![]);
+        let pol2 = Polygon::new(ls2, vec![]);
+        let pol3 = Polygon::new(ls3, vec![]);
+        let mp = MultiPolygon::new(vec![pol1.clone(), pol2, pol3]);
+        let pnt1 = Point::new(0.0, 15.0);
+        let pnt2 = Point::new(10.0, 20.0);
+        let ln = Line::new(pnt1.0, pnt2.0);
+        let dist_mp_ln = Euclidean::distance(&ln, &mp);
+        let dist_pol1_ln = Euclidean::distance(&ln, &pol1);
+        assert_relative_eq!(dist_mp_ln, dist_pol1_ln);
+    }
+
+    #[test]
+    fn point_distance_multipolygon_test() {
+        let ls1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0), (2.0, 0.0), (0.0, 0.0)]);
+        let ls2 = LineString::from(vec![(3.0, 0.0), (4.0, 10.0), (5.0, 0.0), (3.0, 0.0)]);
+        let p1 = Polygon::new(ls1, vec![]);
+        let p2 = Polygon::new(ls2, vec![]);
+        let mp = MultiPolygon::new(vec![p1, p2]);
+        let p = Point::new(50.0, 50.0);
+        assert_relative_eq!(Euclidean::distance(&p, &mp), 60.959002616512684);
+    }
+    #[test]
+    // Point to LineString
+    fn point_linestring_distance_test() {
+        // like an octagon, but missing the lowest horizontal segment
+        let points = vec![
+            (5., 1.),
+            (4., 2.),
+            (4., 3.),
+            (5., 4.),
+            (6., 4.),
+            (7., 3.),
+            (7., 2.),
+            (6., 1.),
+        ];
+        let ls = LineString::from(points);
+        // A Random point "inside" the LineString
+        let p = Point::new(5.5, 2.1);
+        let dist = Euclidean::distance(&p, &ls);
+        assert_relative_eq!(dist, 1.1313708498984762);
+    }
+    #[test]
+    // Point to LineString, point lies on the LineString
+    fn point_linestring_contains_test() {
+        // like an octagon, but missing the lowest horizontal segment
+        let points = vec![
+            (5., 1.),
+            (4., 2.),
+            (4., 3.),
+            (5., 4.),
+            (6., 4.),
+            (7., 3.),
+            (7., 2.),
+            (6., 1.),
+        ];
+        let ls = LineString::from(points);
+        // A point which lies on the LineString
+        let p = Point::new(5.0, 4.0);
+        let dist = Euclidean::distance(&p, &ls);
+        assert_relative_eq!(dist, 0.0);
+    }
+    #[test]
+    // Point to LineString, closed triangle
+    fn point_linestring_triangle_test() {
+        let points = vec![(3.5, 3.5), (4.4, 2.0), (2.6, 2.0), (3.5, 3.5)];
+        let ls = LineString::from(points);
+        let p = Point::new(3.5, 2.5);
+        let dist = Euclidean::distance(&p, &ls);
+        assert_relative_eq!(dist, 0.5);
+    }
+    #[test]
+    // Point to LineString, empty LineString
+    fn point_linestring_empty_test() {
+        let points = vec![];
+        let ls = LineString::new(points);
+        let p = Point::new(5.0, 4.0);
+        let dist = Euclidean::distance(&p, &ls);
+        assert_relative_eq!(dist, 0.0);
+    }
+    #[test]
+    fn distance_multilinestring_test() {
+        let v1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0)]);
+        let v2 = LineString::from(vec![(1.0, 10.0), (2.0, 0.0), (3.0, 1.0)]);
+        let mls = MultiLineString::new(vec![v1, v2]);
+        let p = Point::new(50.0, 50.0);
+        assert_relative_eq!(Euclidean::distance(&p, &mls), 63.25345840347388);
+    }
+    #[test]
+    fn distance1_test() {
+        assert_relative_eq!(
+            Euclidean::distance(&Point::new(0., 0.), &Point::new(1., 0.)),
+            1.
+        );
+    }
+    #[test]
+    fn distance2_test() {
+        let dist =
+            Euclidean::distance(&Point::new(-72.1235, 42.3521), &Point::new(72.1260, 70.612));
+        assert_relative_eq!(dist, 146.99163308930207);
+    }
+    #[test]
+    fn distance_multipoint_test() {
+        let v = vec![
+            Point::new(0.0, 10.0),
+            Point::new(1.0, 1.0),
+            Point::new(10.0, 0.0),
+            Point::new(1.0, -1.0),
+            Point::new(0.0, -10.0),
+            Point::new(-1.0, -1.0),
+            Point::new(-10.0, 0.0),
+            Point::new(-1.0, 1.0),
+            Point::new(0.0, 10.0),
+        ];
+        let mp = MultiPoint::new(v);
+        let p = Point::new(50.0, 50.0);
+        assert_relative_eq!(Euclidean::distance(&p, &mp), 64.03124237432849)
+    }
+    #[test]
+    fn distance_line_test() {
+        let line0 = Line::from([(0., 0.), (5., 0.)]);
+        let p0 = Point::new(2., 3.);
+        let p1 = Point::new(3., 0.);
+        let p2 = Point::new(6., 0.);
+        assert_relative_eq!(Euclidean::distance(&line0, &p0), 3.);
+        assert_relative_eq!(Euclidean::distance(&p0, &line0), 3.);
+
+        assert_relative_eq!(Euclidean::distance(&line0, &p1), 0.);
+        assert_relative_eq!(Euclidean::distance(&p1, &line0), 0.);
+
+        assert_relative_eq!(Euclidean::distance(&line0, &p2), 1.);
+        assert_relative_eq!(Euclidean::distance(&p2, &line0), 1.);
+    }
+    #[test]
+    fn distance_line_line_test() {
+        let line0 = Line::from([(0., 0.), (5., 0.)]);
+        let line1 = Line::from([(2., 1.), (7., 2.)]);
+        assert_relative_eq!(Euclidean::distance(&line0, &line1), 1.);
+        assert_relative_eq!(Euclidean::distance(&line1, &line0), 1.);
+    }
+    #[test]
+    // See https://github.com/georust/geo/issues/476
+    fn distance_line_polygon_test() {
+        let line = Line::new(
+            coord! {
+                x: -0.17084137691985102,
+                y: 0.8748085493016657,
+            },
+            coord! {
+                x: -0.17084137691985102,
+                y: 0.09858870312437906,
+            },
+        );
+        let poly: Polygon<f64> = polygon![
+            coord! {
+                x: -0.10781391405721802,
+                y: -0.15433610862574643,
+            },
+            coord! {
+                x: -0.7855276236615211,
+                y: 0.23694208404779793,
+            },
+            coord! {
+                x: -0.7855276236615214,
+                y: -0.5456143012992907,
+            },
+            coord! {
+                x: -0.10781391405721802,
+                y: -0.15433610862574643,
+            },
+        ];
+        assert_eq!(Euclidean::distance(&line, &poly), 0.18752558079168907);
+    }
+    #[test]
+    // test edge-vertex minimum distance
+    fn test_minimum_polygon_distance() {
+        let points_raw = [
+            (126., 232.),
+            (126., 212.),
+            (112., 202.),
+            (97., 204.),
+            (87., 215.),
+            (87., 232.),
+            (100., 246.),
+            (118., 247.),
+        ];
+        let points = points_raw
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly1 = Polygon::new(LineString::from(points), vec![]);
+
+        let points_raw_2 = [
+            (188., 231.),
+            (189., 207.),
+            (174., 196.),
+            (164., 196.),
+            (147., 220.),
+            (158., 242.),
+            (177., 242.),
+        ];
+        let points2 = points_raw_2
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly2 = Polygon::new(LineString::from(points2), vec![]);
+        let dist = nearest_neighbour_distance(poly1.exterior(), poly2.exterior());
+        assert_relative_eq!(dist, 21.0);
+    }
+    #[test]
+    // test vertex-vertex minimum distance
+    fn test_minimum_polygon_distance_2() {
+        let points_raw = [
+            (118., 200.),
+            (153., 179.),
+            (106., 155.),
+            (88., 190.),
+            (118., 200.),
+        ];
+        let points = points_raw
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly1 = Polygon::new(LineString::from(points), vec![]);
+
+        let points_raw_2 = [
+            (242., 186.),
+            (260., 146.),
+            (182., 175.),
+            (216., 193.),
+            (242., 186.),
+        ];
+        let points2 = points_raw_2
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly2 = Polygon::new(LineString::from(points2), vec![]);
+        let dist = crate::euclidean_distance::nearest_neighbour_distance(
+            poly1.exterior(),
+            poly2.exterior(),
+        );
+        assert_relative_eq!(dist, 29.274562336608895);
+    }
+    #[test]
+    // test edge-edge minimum distance
+    fn test_minimum_polygon_distance_3() {
+        let points_raw = [
+            (182., 182.),
+            (182., 168.),
+            (138., 160.),
+            (136., 193.),
+            (182., 182.),
+        ];
+        let points = points_raw
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly1 = Polygon::new(LineString::from(points), vec![]);
+
+        let points_raw_2 = [
+            (232., 196.),
+            (234., 150.),
+            (194., 165.),
+            (194., 191.),
+            (232., 196.),
+        ];
+        let points2 = points_raw_2
+            .iter()
+            .map(|e| Point::new(e.0, e.1))
+            .collect::<Vec<_>>();
+        let poly2 = Polygon::new(LineString::from(points2), vec![]);
+        let dist = crate::euclidean_distance::nearest_neighbour_distance(
+            poly1.exterior(),
+            poly2.exterior(),
+        );
+        assert_relative_eq!(dist, 12.0);
+    }
+    #[test]
+    fn test_large_polygon_distance() {
+        let ls = geo_test_fixtures::norway_main::<f64>();
+        let poly1 = Polygon::new(ls, vec![]);
+        let vec2 = vec![
+            (4.921875, 66.33750501996518),
+            (3.69140625, 65.21989393613207),
+            (6.15234375, 65.07213008560697),
+            (4.921875, 66.33750501996518),
+        ];
+        let poly2 = Polygon::new(vec2.into(), vec![]);
+        let distance = Euclidean::distance(&poly1, &poly2);
+        // GEOS says 2.2864896295566055
+        assert_relative_eq!(distance, 2.2864896295566055);
+    }
+    #[test]
+    // A polygon inside another polygon's ring; they're disjoint in the DE-9IM sense:
+    // FF2FF1212
+    fn test_poly_in_ring() {
+        let shell = geo_test_fixtures::shell::<f64>();
+        let ring = geo_test_fixtures::ring::<f64>();
+        let poly_in_ring = geo_test_fixtures::poly_in_ring::<f64>();
+        // inside is "inside" outside's ring, but they are disjoint
+        let outside = Polygon::new(shell, vec![ring]);
+        let inside = Polygon::new(poly_in_ring, vec![]);
+        assert_relative_eq!(Euclidean::distance(&outside, &inside), 5.992772737231033);
+    }
+    #[test]
+    // two ring LineStrings; one encloses the other but they neither touch nor intersect
+    fn test_linestring_distance() {
+        let ring = geo_test_fixtures::ring::<f64>();
+        let poly_in_ring = geo_test_fixtures::poly_in_ring::<f64>();
+        assert_relative_eq!(Euclidean::distance(&ring, &poly_in_ring), 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::from([(0.0, 0.0), (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_relative_eq!(Euclidean::distance(&line, &poly), 0.25);
+    }
+    #[test]
+    // Line-Polygon test: Line intersects Polygon
+    fn test_line_polygon_intersects() {
+        let line = Line::from([(0.5, 0.0), (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_relative_eq!(Euclidean::distance(&line, &poly), 0.0);
+    }
+    #[test]
+    // Line-Polygon test: Line contained by interior ring
+    fn test_line_polygon_inside_ring() {
+        let line = Line::from([(4.4, 1.5), (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_relative_eq!(Euclidean::distance(&line, &poly), 0.04999999999999982);
+    }
+    #[test]
+    // LineString-Line test
+    fn test_linestring_line_distance() {
+        let line = Line::from([(0.0, 0.0), (0.0, 2.0)]);
+        let ls: LineString<_> = vec![(3.0, 0.0), (1.0, 1.0), (3.0, 2.0)].into();
+        assert_relative_eq!(Euclidean::distance(&ls, &line), 1.0);
+    }
+
+    #[test]
+    // Triangle-Point test: point on vertex
+    fn test_triangle_point_on_vertex_distance() {
+        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
+        let point = Point::new(0.0, 0.0);
+        assert_relative_eq!(Euclidean::distance(&triangle, &point), 0.0);
+    }
+
+    #[test]
+    // Triangle-Point test: point on edge
+    fn test_triangle_point_on_edge_distance() {
+        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
+        let point = Point::new(1.5, 0.0);
+        assert_relative_eq!(Euclidean::distance(&triangle, &point), 0.0);
+    }
+
+    #[test]
+    // Triangle-Point test
+    fn test_triangle_point_distance() {
+        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
+        let point = Point::new(2.0, 3.0);
+        assert_relative_eq!(Euclidean::distance(&triangle, &point), 1.0);
+    }
+
+    #[test]
+    // Triangle-Point test: point within triangle
+    fn test_triangle_point_inside_distance() {
+        let triangle = Triangle::from([(0.0, 0.0), (2.0, 0.0), (2.0, 2.0)]);
+        let point = Point::new(1.0, 0.5);
+        assert_relative_eq!(Euclidean::distance(&triangle, &point), 0.0);
+    }
+
+    #[test]
+    fn convex_and_nearest_neighbour_comparison() {
+        let ls1: LineString<f64> = vec![
+            Coord::from((57.39453770777941, 307.60533608924663)),
+            Coord::from((67.1800355576469, 309.6654408997451)),
+            Coord::from((84.89693692793338, 225.5101593908847)),
+            Coord::from((75.1114390780659, 223.45005458038628)),
+            Coord::from((57.39453770777941, 307.60533608924663)),
+        ]
+        .into();
+        let first_polygon: Polygon<f64> = Polygon::new(ls1, vec![]);
+        let ls2: LineString<f64> = vec![
+            Coord::from((138.11769866645008, -45.75134112915392)),
+            Coord::from((130.50230476949187, -39.270154833870336)),
+            Coord::from((184.94426964987397, 24.699153900578573)),
+            Coord::from((192.55966354683218, 18.217967605294987)),
+            Coord::from((138.11769866645008, -45.75134112915392)),
+        ]
+        .into();
+        let second_polygon = Polygon::new(ls2, vec![]);
+
+        assert_relative_eq!(
+            Euclidean::distance(&first_polygon, &second_polygon),
+            224.35357967013238
+        );
+    }
+    #[test]
+    fn fast_path_regression() {
+        // this test will fail if the fast path algorithm is reintroduced without being fixed
+        let p1 = polygon!(
+            (x: 0_f64, y: 0_f64),
+            (x: 300_f64, y: 0_f64),
+            (x: 300_f64, y: 100_f64),
+            (x: 0_f64, y: 100_f64),
+        )
+        .orient(Direction::Default);
+        let p2 = polygon!(
+            (x: 100_f64, y: 150_f64),
+            (x: 150_f64, y: 200_f64),
+            (x: 50_f64, y: 200_f64),
+        )
+        .orient(Direction::Default);
+        let p3 = polygon!(
+            (x: 0_f64, y: 0_f64),
+            (x: 300_f64, y: 0_f64),
+            (x: 300_f64, y: 100_f64),
+            (x: 0_f64, y: 100_f64),
+        )
+        .orient(Direction::Reversed);
+        let p4 = polygon!(
+            (x: 100_f64, y: 150_f64),
+            (x: 150_f64, y: 200_f64),
+            (x: 50_f64, y: 200_f64),
+        )
+        .orient(Direction::Reversed);
+        assert_eq!(Euclidean::distance(&p1, &p2), 50.0f64);
+        assert_eq!(Euclidean::distance(&p3, &p4), 50.0f64);
+        assert_eq!(Euclidean::distance(&p1, &p4), 50.0f64);
+        assert_eq!(Euclidean::distance(&p2, &p3), 50.0f64);
+    }
+    #[test]
+    fn all_types_geometry_collection_test() {
+        let p = Point::new(0.0, 0.0);
+        let line = Line::from([(-1.0, -1.0), (-2.0, -2.0)]);
+        let ls = LineString::from(vec![(0.0, 0.0), (1.0, 10.0), (2.0, 0.0)]);
+        let poly = Polygon::new(
+            LineString::from(vec![(0.0, 0.0), (1.0, 10.0), (2.0, 0.0), (0.0, 0.0)]),
+            vec![],
+        );
+        let tri = Triangle::from([(0.0, 0.0), (1.0, 10.0), (2.0, 0.0)]);
+        let rect = Rect::new((0.0, 0.0), (-1.0, -1.0));
+
+        let ls1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0), (2.0, 0.0), (0.0, 0.0)]);
+        let ls2 = LineString::from(vec![(3.0, 0.0), (4.0, 10.0), (5.0, 0.0), (3.0, 0.0)]);
+        let p1 = Polygon::new(ls1, vec![]);
+        let p2 = Polygon::new(ls2, vec![]);
+        let mpoly = MultiPolygon::new(vec![p1, p2]);
+
+        let v = vec![
+            Point::new(0.0, 10.0),
+            Point::new(1.0, 1.0),
+            Point::new(10.0, 0.0),
+            Point::new(1.0, -1.0),
+            Point::new(0.0, -10.0),
+            Point::new(-1.0, -1.0),
+            Point::new(-10.0, 0.0),
+            Point::new(-1.0, 1.0),
+            Point::new(0.0, 10.0),
+        ];
+        let mpoint = MultiPoint::new(v);
+
+        let v1 = LineString::from(vec![(0.0, 0.0), (1.0, 10.0)]);
+        let v2 = LineString::from(vec![(1.0, 10.0), (2.0, 0.0), (3.0, 1.0)]);
+        let mls = MultiLineString::new(vec![v1, v2]);
+
+        let gc = GeometryCollection(vec![
+            Geometry::Point(p),
+            Geometry::Line(line),
+            Geometry::LineString(ls),
+            Geometry::Polygon(poly),
+            Geometry::MultiPoint(mpoint),
+            Geometry::MultiLineString(mls),
+            Geometry::MultiPolygon(mpoly),
+            Geometry::Triangle(tri),
+            Geometry::Rect(rect),
+        ]);
+
+        let test_p = Point::new(50., 50.);
+        assert_relative_eq!(Euclidean::distance(&test_p, &gc), 60.959002616512684);
+
+        let test_multipoint = MultiPoint::new(vec![test_p]);
+        assert_relative_eq!(
+            Euclidean::distance(&test_multipoint, &gc),
+            60.959002616512684
+        );
+
+        let test_line = Line::from([(50., 50.), (60., 60.)]);
+        assert_relative_eq!(Euclidean::distance(&test_line, &gc), 60.959002616512684);
+
+        let test_ls = LineString::from(vec![(50., 50.), (60., 60.), (70., 70.)]);
+        assert_relative_eq!(Euclidean::distance(&test_ls, &gc), 60.959002616512684);
+
+        let test_mls = MultiLineString::new(vec![test_ls]);
+        assert_relative_eq!(Euclidean::distance(&test_mls, &gc), 60.959002616512684);
+
+        let test_poly = Polygon::new(
+            LineString::from(vec![
+                (50., 50.),
+                (60., 50.),
+                (60., 60.),
+                (55., 55.),
+                (50., 50.),
+            ]),
+            vec![],
+        );
+        assert_relative_eq!(Euclidean::distance(&test_poly, &gc), 60.959002616512684);
+
+        let test_multipoly = MultiPolygon::new(vec![test_poly]);
+        assert_relative_eq!(
+            Euclidean::distance(&test_multipoly, &gc),
+            60.959002616512684
+        );
+
+        let test_tri = Triangle::from([(50., 50.), (60., 50.), (55., 55.)]);
+        assert_relative_eq!(Euclidean::distance(&test_tri, &gc), 60.959002616512684);
+
+        let test_rect = Rect::new(coord! { x: 50., y: 50. }, coord! { x: 60., y: 60. });
+        assert_relative_eq!(Euclidean::distance(&test_rect, &gc), 60.959002616512684);
+
+        let test_gc = GeometryCollection(vec![Geometry::Rect(test_rect)]);
+        assert_relative_eq!(Euclidean::distance(&test_gc, &gc), 60.959002616512684);
+    }
+}
diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs
similarity index 63%
rename from geo/src/algorithm/line_measures/metric_spaces/euclidean.rs
rename to geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs
index 75e2c90f92..715e211328 100644
--- a/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs
+++ b/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs
@@ -1,3 +1,5 @@
+mod distance;
+
 use super::super::{Distance, InterpolatePoint};
 use crate::line_measures::densify::densify_between;
 use crate::{CoordFloat, Point};
@@ -20,41 +22,6 @@ use num_traits::FromPrimitive;
 /// [metric spaces]: super
 pub struct Euclidean;
 
-/// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points
-impl<F: CoordFloat> Distance<F, Point<F>, Point<F>> for Euclidean {
-    /// Calculate the Euclidean distance (a.k.a. pythagorean distance) between two Points
-    ///
-    /// # Units
-    /// - `origin`, `destination`: Point where the units of x/y represent non-angular units
-    ///    — e.g. meters or miles, not lon/lat. For lon/lat points, use the
-    ///    [`Haversine`] or [`Geodesic`] [metric spaces].
-    /// - returns: distance in the same units as the `origin` and `destination` points
-    ///
-    /// # Example
-    /// ```
-    /// use geo::{Euclidean, Distance};
-    /// use geo::Point;
-    /// // web mercator
-    /// let new_york_city = Point::new(-8238310.24, 4942194.78);
-    /// // web mercator
-    /// let london = Point::new(-14226.63, 6678077.70);
-    /// let distance: f64 = Euclidean::distance(new_york_city, london);
-    ///
-    /// assert_eq!(
-    ///     8_405_286., // meters in web mercator
-    ///     distance.round()
-    /// );
-    /// ```
-    ///
-    /// [`Haversine`]: super::Haversine
-    /// [`Geodesic`]: super::Geodesic
-    /// [metric spaces]: super
-    fn distance(origin: Point<F>, destination: Point<F>) -> F {
-        let delta = origin - destination;
-        delta.x().hypot(delta.y())
-    }
-}
-
 impl<F: CoordFloat + FromPrimitive> InterpolatePoint<F> for Euclidean {
     fn point_at_ratio_between(start: Point<F>, end: Point<F>, ratio_from_start: F) -> Point<F> {
         let diff = end - start;
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 32e140c7a6..344084dc3a 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -77,6 +77,7 @@ pub use dimensions::HasDimensions;
 
 /// Calculate the minimum Euclidean distance between two `Geometries`.
 pub mod euclidean_distance;
+#[allow(deprecated)]
 pub use euclidean_distance::EuclideanDistance;
 
 /// Calculate the length of a planar line between two `Geometries`.
diff --git a/geo/src/algorithm/relate/geomgraph/mod.rs b/geo/src/algorithm/relate/geomgraph/mod.rs
index 62af0bee06..19906fd65d 100644
--- a/geo/src/algorithm/relate/geomgraph/mod.rs
+++ b/geo/src/algorithm/relate/geomgraph/mod.rs
@@ -18,8 +18,8 @@ pub(crate) use quadrant::Quadrant;
 pub(crate) use robust_line_intersector::RobustLineIntersector;
 use topology_position::TopologyPosition;
 
+pub use crate::coordinate_position::CoordPos;
 use crate::dimensions::Dimensions;
-pub use crate::utils::CoordPos;
 
 mod edge;
 mod edge_end;
diff --git a/geo/src/algorithm/simplify.rs b/geo/src/algorithm/simplify.rs
index c5cdb07dab..83e9ec41e4 100644
--- a/geo/src/algorithm/simplify.rs
+++ b/geo/src/algorithm/simplify.rs
@@ -1,5 +1,6 @@
-use crate::{Coord, GeoFloat, Line, LineString, MultiLineString, MultiPolygon, Polygon};
-use crate::{CoordsIter, EuclideanDistance};
+use crate::algorithm::{CoordsIter, Distance, Euclidean};
+use crate::geometry::{Coord, Line, LineString, MultiLineString, MultiPolygon, Polygon};
+use crate::GeoFloat;
 
 const LINE_STRING_INITIAL_MIN: usize = 2;
 const POLYGON_INITIAL_MIN: usize = 4;
@@ -96,7 +97,12 @@ where
         .enumerate()
         .take(rdp_indices.len() - 1) // Don't include the last index
         .skip(1) // Don't include the first index
-        .map(|(index, rdp_index)| (index, rdp_index.coord.euclidean_distance(&first_last_line)))
+        .map(|(index, rdp_index)| {
+            (
+                index,
+                Euclidean::distance(rdp_index.coord, &first_last_line),
+            )
+        })
         .fold(
             (0usize, T::zero()),
             |(farthest_index, farthest_distance), (index, distance)| {
diff --git a/geo/src/algorithm/triangulate_spade.rs b/geo/src/algorithm/triangulate_spade.rs
index 90bb05c40c..1f5fa88586 100644
--- a/geo/src/algorithm/triangulate_spade.rs
+++ b/geo/src/algorithm/triangulate_spade.rs
@@ -4,7 +4,7 @@ use spade::{
 };
 
 use crate::{
-    line_intersection::line_intersection, CoordsIter, EuclideanDistance, GeoFloat,
+    line_intersection::line_intersection, CoordsIter, Distance, Euclidean, GeoFloat,
     LineIntersection, LinesIter,
 };
 use crate::{Centroid, Contains};
@@ -530,12 +530,12 @@ fn snap_or_register_point<T: SpadeTriangulationFloat>(
         .iter()
         // find closest
         .min_by(|a, b| {
-            a.euclidean_distance(&point)
-                .partial_cmp(&b.euclidean_distance(&point))
+            Euclidean::distance(**a, point)
+                .partial_cmp(&Euclidean::distance(**b, point))
                 .expect("Couldn't compare coordinate distances")
         })
-        // only snap if closest is within epsilone range
-        .filter(|nearest_point| nearest_point.euclidean_distance(&point) < snap_radius)
+        // only snap if closest is within epsilon range
+        .filter(|nearest_point| Euclidean::distance(**nearest_point, point) < snap_radius)
         .cloned()
         // otherwise register and use input point
         .unwrap_or_else(|| {
diff --git a/geo/src/types.rs b/geo/src/types.rs
index 03f7cb4b20..33ab22f988 100644
--- a/geo/src/types.rs
+++ b/geo/src/types.rs
@@ -16,7 +16,7 @@ impl<F: GeoFloat> Closest<F> {
     /// Compare two `Closest`s relative to `p` and return a copy of the best
     /// one.
     pub fn best_of_two(&self, other: &Self, p: Point<F>) -> Self {
-        use crate::EuclideanDistance;
+        use crate::{Distance, Euclidean};
 
         let left = match *self {
             Closest::Indeterminate => return *other,
@@ -29,7 +29,7 @@ impl<F: GeoFloat> Closest<F> {
             Closest::SinglePoint(r) => r,
         };
 
-        if left.euclidean_distance(&p) <= right.euclidean_distance(&p) {
+        if Euclidean::distance(left, p) <= Euclidean::distance(right, p) {
             *self
         } else {
             *other
diff --git a/geo/src/utils.rs b/geo/src/utils.rs
index 9db9a3cd42..15720364ce 100644
--- a/geo/src/utils.rs
+++ b/geo/src/utils.rs
@@ -93,9 +93,6 @@ pub fn partial_min<T: PartialOrd>(a: T, b: T) -> T {
     }
 }
 
-// Moved to their own module, but we re-export to avoid breaking the API.
-pub use crate::coordinate_position::{coord_pos_relative_to_ring, CoordPos};
-
 use std::cmp::Ordering;
 
 /// Compare two coordinates lexicographically: first by the