From 44dd23143bff184f22cc539330a28ed2138e4fbb Mon Sep 17 00:00:00 2001
From: Michael Kirk <michael.code@endoftheworl.de>
Date: Wed, 25 Sep 2024 16:55:28 -0700
Subject: [PATCH] Unify `Densify` trait across metric spaces

Just like #1228, but for Densify rather than Length.

This enables Geodesic, Haversine, Euclidean, and Rhumb Densification.

This required implementing InterpolatePoint for Euclidean.

Adjacent work:

deprecated legacy DensifyHaversine
linestring_segment is now implemented in terms of the new Densify

NOTE: linestring_segment would be a good future candidate for
a similar unification across all the metric spaces
---
 geo/CHANGES.md                                |  14 +
 geo/src/algorithm/densify.rs                  | 263 -----------
 geo/src/algorithm/densify_haversine.rs        |  80 +---
 geo/src/algorithm/line_measures/densify.rs    | 437 ++++++++++++++++++
 .../line_measures/metric_spaces/euclidean.rs  |  28 +-
 .../line_measures/metric_spaces/geodesic.rs   |   2 +-
 .../line_measures/metric_spaces/haversine.rs  |   2 +-
 .../line_measures/metric_spaces/rhumb.rs      |   2 +
 geo/src/algorithm/line_measures/mod.rs        |   3 +
 geo/src/algorithm/linestring_segment.rs       |  18 +-
 geo/src/algorithm/mod.rs                      |   7 +-
 11 files changed, 513 insertions(+), 343 deletions(-)
 delete mode 100644 geo/src/algorithm/densify.rs
 create mode 100644 geo/src/algorithm/line_measures/densify.rs

diff --git a/geo/CHANGES.md b/geo/CHANGES.md
index 62c48fa33..4b7c7895d 100644
--- a/geo/CHANGES.md
+++ b/geo/CHANGES.md
@@ -53,6 +53,20 @@
   line_string.length::<Haversine>();
   ```
   * <https://github.com/georust/geo/pull/1228>
+* Deprecated `DensifyHaversine`
+* BREAKING: `Densify::densify` is no longer strictly Euclidean, and now accepts a generic line measure parameter.
+   ```
+  // Before
+  line_string.densify();
+  line_string.densify_haversine();
+  // After
+  line_string.densify::<Euclidean>();
+  line_string.densify::<Haversine>();
+
+  // Additional measures are now supported
+  line_string.densify::<Geodesic>();
+  line_string.densify::<Rhumb>();
+   ```
 * Change IntersectionMatrix::is_equal_topo to now consider empty geometries as equal.
   * <https://github.com/georust/geo/pull/1223>
 * Fix `(LINESTRING EMPTY).contains(LINESTRING EMPTY)` and `(MULTIPOLYGON EMPTY).contains(MULTIPOINT EMPTY)` which previously
diff --git a/geo/src/algorithm/densify.rs b/geo/src/algorithm/densify.rs
deleted file mode 100644
index 3a930f904..000000000
--- a/geo/src/algorithm/densify.rs
+++ /dev/null
@@ -1,263 +0,0 @@
-use crate::{
-    CoordFloat, Line, LineInterpolatePoint, LineString, MultiLineString, MultiPolygon, Point,
-    Polygon, Rect, Triangle,
-};
-
-// This is still used in the trait constraints - but Densify too will soon be replaced with a
-// generic version, at which point this implementation detail can be removed.
-#[allow(deprecated)]
-use crate::EuclideanLength;
-
-/// Return a new linear geometry containing both existing and new interpolated coordinates with
-/// a maximum distance of `max_distance` between them.
-///
-/// Note: `max_distance` must be greater than 0.
-///
-/// # Examples
-/// ```
-/// use geo::{coord, Line, LineString};
-/// use geo::Densify;
-///
-/// let line: Line<f64> = Line::new(coord! {x: 0.0, y: 6.0}, coord! {x: 1.0, y: 8.0});
-/// let correct: LineString<f64> = vec![[0.0, 6.0], [0.5, 7.0], [1.0, 8.0]].into();
-/// let max_dist = 2.0;
-/// let densified = line.densify(max_dist);
-/// assert_eq!(densified, correct);
-///```
-pub trait Densify<F: CoordFloat> {
-    type Output;
-
-    fn densify(&self, max_distance: F) -> Self::Output;
-}
-
-// Helper for densification trait
-#[allow(deprecated)]
-fn densify_line<T: CoordFloat>(line: Line<T>, container: &mut Vec<Point<T>>, max_distance: T) {
-    assert!(max_distance > T::zero());
-    container.push(line.start_point());
-    let num_segments = (line.euclidean_length() / max_distance)
-        .ceil()
-        .to_u64()
-        .unwrap();
-    // distance "unit" for this line segment
-    let frac = T::one() / T::from(num_segments).unwrap();
-    for segment_idx in 1..num_segments {
-        let ratio = frac * T::from(segment_idx).unwrap();
-        let interpolated_point = line
-            .line_interpolate_point(ratio)
-            .expect("ratio should be between 0..1");
-        container.push(interpolated_point);
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for MultiPolygon<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = MultiPolygon<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        MultiPolygon::new(
-            self.iter()
-                .map(|polygon| polygon.densify(max_distance))
-                .collect(),
-        )
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for Polygon<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = Polygon<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        let densified_exterior = self.exterior().densify(max_distance);
-        let densified_interiors = self
-            .interiors()
-            .iter()
-            .map(|ring| ring.densify(max_distance))
-            .collect();
-        Polygon::new(densified_exterior, densified_interiors)
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for MultiLineString<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = MultiLineString<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        MultiLineString::new(
-            self.iter()
-                .map(|linestring| linestring.densify(max_distance))
-                .collect(),
-        )
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for LineString<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = LineString<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        if self.0.is_empty() {
-            return LineString::new(vec![]);
-        }
-
-        let mut new_line = vec![];
-
-        self.lines()
-            .for_each(|line| densify_line(line, &mut new_line, max_distance));
-        // we're done, push the last coordinate on to finish
-        new_line.push(self.points().last().unwrap());
-        LineString::from(new_line)
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for Line<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = LineString<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        let mut new_line = vec![];
-        densify_line(*self, &mut new_line, max_distance);
-        // we're done, push the last coordinate on to finish
-        new_line.push(self.end_point());
-        LineString::from(new_line)
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for Triangle<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = Polygon<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        self.to_polygon().densify(max_distance)
-    }
-}
-
-#[allow(deprecated)]
-impl<T> Densify<T> for Rect<T>
-where
-    T: CoordFloat,
-    Line<T>: EuclideanLength<T>,
-    LineString<T>: EuclideanLength<T>,
-{
-    type Output = Polygon<T>;
-
-    fn densify(&self, max_distance: T) -> Self::Output {
-        self.to_polygon().densify(max_distance)
-    }
-}
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::{coord, Coord};
-
-    #[test]
-    fn test_polygon_densify() {
-        let linestring: LineString<f64> =
-            vec![[-5.0, 0.0], [0.0, 5.0], [5.0, 0.0], [-5.0, 0.0]].into();
-        let interior: LineString<f64> =
-            vec![[-3.0, 0.0], [0.0, 3.0], [3.0, 0.0], [-3.0, 0.0]].into();
-        let polygon = Polygon::new(linestring, vec![interior]);
-        let correct_ext: LineString<f64> = LineString(vec![
-            Coord { x: -5.0, y: 0.0 },
-            Coord { x: -3.75, y: 1.25 },
-            Coord { x: -2.5, y: 2.5 },
-            Coord { x: -1.25, y: 3.75 },
-            Coord { x: 0.0, y: 5.0 },
-            Coord { x: 1.25, y: 3.75 },
-            Coord { x: 2.5, y: 2.5 },
-            Coord { x: 3.75, y: 1.25 },
-            Coord { x: 5.0, y: 0.0 },
-            Coord { x: 3.0, y: 0.0 },
-            Coord { x: 1.0, y: 0.0 },
-            Coord {
-                x: -1.0000000000000009,
-                y: 0.0,
-            },
-            Coord { x: -3.0, y: 0.0 },
-            Coord { x: -5.0, y: 0.0 },
-        ]);
-        let correct_int: LineString<f64> = LineString(vec![
-            Coord { x: -3.0, y: 0.0 },
-            Coord { x: -2.0, y: 1.0 },
-            Coord { x: -1.0, y: 2.0 },
-            Coord { x: 0.0, y: 3.0 },
-            Coord { x: 1.0, y: 2.0 },
-            Coord { x: 2.0, y: 1.0 },
-            Coord { x: 3.0, y: 0.0 },
-            Coord { x: 1.0, y: 0.0 },
-            Coord { x: -1.0, y: 0.0 },
-            Coord { x: -3.0, y: 0.0 },
-        ]);
-        let correct_polygon = Polygon::new(correct_ext, vec![correct_int]);
-        let max_dist = 2.0;
-        let densified = polygon.densify(max_dist);
-        assert_eq!(densified, correct_polygon);
-    }
-
-    #[test]
-    fn test_empty_linestring_densify() {
-        let linestring = LineString::<f64>::new(vec![]);
-        let max_dist = 2.0;
-        let densified = linestring.densify(max_dist);
-        assert!(densified.0.is_empty());
-    }
-
-    #[test]
-    fn test_linestring_densify() {
-        let linestring: LineString<f64> =
-            vec![[-1.0, 0.0], [0.0, 0.0], [0.0, 6.0], [1.0, 8.0]].into();
-        let correct: LineString<f64> = vec![
-            [-1.0, 0.0],
-            [0.0, 0.0],
-            [0.0, 2.0],
-            [0.0, 4.0],
-            [0.0, 6.0],
-            [0.5, 7.0],
-            [1.0, 8.0],
-        ]
-        .into();
-        let max_dist = 2.0;
-        let densified = linestring.densify(max_dist);
-        assert_eq!(densified, correct);
-    }
-
-    #[test]
-    fn test_line_densify() {
-        let line: Line<f64> = Line::new(coord! {x: 0.0, y: 6.0}, coord! {x: 1.0, y: 8.0});
-        let correct: LineString<f64> = vec![[0.0, 6.0], [0.5, 7.0], [1.0, 8.0]].into();
-        let max_dist = 2.0;
-        let densified = line.densify(max_dist);
-        assert_eq!(densified, correct);
-    }
-}
diff --git a/geo/src/algorithm/densify_haversine.rs b/geo/src/algorithm/densify_haversine.rs
index e8b2f0dbc..754779d83 100644
--- a/geo/src/algorithm/densify_haversine.rs
+++ b/geo/src/algorithm/densify_haversine.rs
@@ -1,14 +1,17 @@
 use num_traits::FromPrimitive;
 
-use crate::line_measures::{Haversine, InterpolatePoint, Length};
+use crate::line_measures::Haversine;
 // Densify will soon be deprecated too, so let's just allow deprecated for now
 #[allow(deprecated)]
 use crate::HaversineLength;
 use crate::{
-    CoordFloat, CoordsIter, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Rect,
-    Triangle,
+    CoordFloat, Densify, Line, LineString, MultiLineString, MultiPolygon, Polygon, Rect, Triangle,
 };
 
+#[deprecated(
+    since = "0.29.0",
+    note = "Please use the `line.densify::<Haversine>()` via the `Densify` trait instead."
+)]
 /// Returns a new spherical geometry containing both existing and new interpolated coordinates with
 /// a maximum distance of `max_distance` between them.
 ///
@@ -21,6 +24,7 @@ use crate::{
 /// # Examples
 /// ```
 /// use geo::{coord, Line, LineString};
+/// #[allow(deprecated)]
 /// use geo::DensifyHaversine;
 ///
 /// let line = Line::new(coord! {x: 0.0, y: 0.0}, coord! { x: 0.0, y: 1.0 });
@@ -36,29 +40,6 @@ pub trait DensifyHaversine<F: CoordFloat> {
     fn densify_haversine(&self, max_distance: F) -> Self::Output;
 }
 
-// Helper for densification trait
-fn densify_line<T: CoordFloat + FromPrimitive>(
-    line: Line<T>,
-    container: &mut Vec<Point<T>>,
-    max_distance: T,
-) {
-    assert!(max_distance > T::zero());
-    container.push(line.start_point());
-    let num_segments = (line.length::<Haversine>() / max_distance)
-        .ceil()
-        .to_u64()
-        .unwrap();
-    // distance "unit" for this line segment
-    let frac = T::one() / T::from(num_segments).unwrap();
-    for segment_idx in 1..num_segments {
-        let ratio = frac * T::from(segment_idx).unwrap();
-        let start = line.start;
-        let end = line.end;
-        let interpolated_point = Haversine::point_at_ratio_between(Point(start), Point(end), ratio);
-        container.push(interpolated_point);
-    }
-}
-
 #[allow(deprecated)]
 impl<T> DensifyHaversine<T> for MultiPolygon<T>
 where
@@ -69,11 +50,7 @@ where
     type Output = MultiPolygon<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        MultiPolygon::new(
-            self.iter()
-                .map(|polygon| polygon.densify_haversine(max_distance))
-                .collect(),
-        )
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -87,13 +64,7 @@ where
     type Output = Polygon<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        let densified_exterior = self.exterior().densify_haversine(max_distance);
-        let densified_interiors = self
-            .interiors()
-            .iter()
-            .map(|ring| ring.densify_haversine(max_distance))
-            .collect();
-        Polygon::new(densified_exterior, densified_interiors)
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -107,11 +78,7 @@ where
     type Output = MultiLineString<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        MultiLineString::new(
-            self.iter()
-                .map(|linestring| linestring.densify_haversine(max_distance))
-                .collect(),
-        )
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -125,16 +92,7 @@ where
     type Output = LineString<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        if self.coords_count() == 0 {
-            return LineString::new(vec![]);
-        }
-
-        let mut new_line = vec![];
-        self.lines()
-            .for_each(|line| densify_line(line, &mut new_line, max_distance));
-        // we're done, push the last coordinate on to finish
-        new_line.push(self.points().last().unwrap());
-        LineString::from(new_line)
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -148,11 +106,7 @@ where
     type Output = LineString<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        let mut new_line = vec![];
-        densify_line(*self, &mut new_line, max_distance);
-        // we're done, push the last coordinate on to finish
-        new_line.push(self.end_point());
-        LineString::from(new_line)
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -166,7 +120,7 @@ where
     type Output = Polygon<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        self.to_polygon().densify_haversine(max_distance)
+        self.densify::<Haversine>(max_distance)
     }
 }
 
@@ -180,14 +134,14 @@ where
     type Output = Polygon<T>;
 
     fn densify_haversine(&self, max_distance: T) -> Self::Output {
-        self.to_polygon().densify_haversine(max_distance)
+        self.densify::<Haversine>(max_distance)
     }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::coord;
+    use crate::{coord, CoordsIter};
 
     #[test]
     fn test_polygon_densify() {
@@ -218,6 +172,7 @@ mod tests {
         ]
         .into();
 
+        #[allow(deprecated)]
         let dense = polygon.densify_haversine(50000.0);
         assert_relative_eq!(dense.exterior(), &output_exterior);
     }
@@ -250,6 +205,7 @@ mod tests {
         ]
         .into();
 
+        #[allow(deprecated)]
         let dense = linestring.densify_haversine(110.0);
         assert_relative_eq!(dense, output);
     }
@@ -258,6 +214,7 @@ mod tests {
     fn test_line_densify() {
         let output: LineString = vec![[0.0, 0.0], [0.0, 0.5], [0.0, 1.0]].into();
         let line = Line::new(coord! {x: 0.0, y: 0.0}, coord! { x: 0.0, y: 1.0 });
+        #[allow(deprecated)]
         let dense = line.densify_haversine(100000.0);
         assert_relative_eq!(dense, output);
     }
@@ -265,6 +222,7 @@ mod tests {
     #[test]
     fn test_empty_linestring() {
         let linestring: LineString<f64> = LineString::new(vec![]);
+        #[allow(deprecated)]
         let dense = linestring.densify_haversine(10.0);
         assert_eq!(0, dense.coords_count());
     }
diff --git a/geo/src/algorithm/line_measures/densify.rs b/geo/src/algorithm/line_measures/densify.rs
new file mode 100644
index 000000000..b0a3eb08b
--- /dev/null
+++ b/geo/src/algorithm/line_measures/densify.rs
@@ -0,0 +1,437 @@
+use super::{Distance, InterpolatePoint};
+use crate::{
+    CoordFloat, CoordsIter, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Rect,
+    Triangle,
+};
+use num_traits::FromPrimitive;
+
+/// Creates a copy of the geometry with additional points inserted as necessary to ensure there
+/// is never more than `max_segment_length` between points.
+///
+/// ## Units
+/// - `max_segment_length` units depend on the implementing [metric space]. It must be greater than 0.
+///
+/// # Examples
+/// ```
+/// # use approx::assert_relative_eq;
+/// use geo::{wkt, Densify};
+/// use geo::line_measures::Euclidean;
+///
+/// let line_string = wkt!(LINESTRING(0.0 0.0,0.0 6.0,1.0 7.0));
+///
+/// // For Euclidean calculations, the unit of distance is the same as the units
+/// // of your coordinates.
+/// let max_dist = 2.0;
+/// let densified = line_string.densify::<Euclidean>(max_dist);
+/// let expected_output = wkt!(LINESTRING(
+///     0.0 0.0,
+///     0.0 2.0,
+///     0.0 4.0,
+///     0.0 6.0,
+///     1.0 7.0
+/// ));
+/// assert_relative_eq!(densified, expected_output);
+///```
+///
+/// For lng/lat geometries, consider using a different [metric space] like [`Haversine`](crate::Haversine) or [`Geodesic`](crate::Geodesic).
+///
+/// ```
+/// # use approx::assert_relative_eq;
+/// use geo::{wkt, Densify};
+/// use geo::line_measures::Haversine;
+/// let line_string = wkt!(LINESTRING(0.0 0.0,0.0 6.0,1.0 7.0));
+///
+/// // For Haversine, the unit of distance is in meters
+/// let max_dist = 200_000.0;
+/// let densified = line_string.densify::<Haversine>(max_dist);
+/// // Haversine interprets coordinate points as lng/lat
+/// let expected_output = wkt!(LINESTRING(
+///     0.0 0.0,
+///     0.0 1.5,
+///     0.0 3.0,
+///     0.0 4.5,
+///     0.0 6.0,
+///     1.0 7.0
+/// ));
+/// assert_relative_eq!(densified, expected_output, epsilon = 1e-14);
+/// ```
+/// [metric space]: crate::line_measures::metric_spaces
+pub trait Densify<F: CoordFloat> {
+    type Output;
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>;
+}
+
+pub(crate) fn densify_between<F, MetricSpace>(
+    line_start: Point<F>,
+    line_end: Point<F>,
+    container: &mut Vec<Point<F>>,
+    max_segment_length: F,
+) where
+    F: CoordFloat + FromPrimitive,
+    MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+{
+    assert!(max_segment_length > F::zero());
+    let num_segments = (MetricSpace::distance(line_start, line_end) / max_segment_length)
+        .ceil()
+        .to_u64()
+        .expect("unreasonable number of segments");
+
+    // distance "unit" for this line segment
+    let frac = F::one() / F::from(num_segments).unwrap();
+
+    for segment_num in 1..num_segments {
+        let ratio = frac * F::from(segment_num).unwrap();
+
+        // PERF TODO: We recompute "total_distance" every step of this loop.
+        // If we impl point_at_distance_between, we could compute it once and use it here.
+        // At that point, I think this function could be a good candidate to be *the single* basis
+        // for a unified generic of points_along_line for all metric spaces.
+        let interpolated_point = MetricSpace::point_at_ratio_between(line_start, line_end, ratio);
+        container.push(interpolated_point);
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for Line<F> {
+    type Output = LineString<F>;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        let mut points = vec![self.start_point()];
+        densify_between::<F, MetricSpace>(
+            self.start_point(),
+            self.end_point(),
+            &mut points,
+            max_segment_length,
+        );
+        points.push(self.end_point());
+        LineString::from(points)
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for LineString<F> {
+    type Output = Self;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> LineString<F>
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        if self.coords_count() == 0 {
+            return LineString::new(vec![]);
+        }
+
+        let mut points = vec![];
+        self.lines().for_each(|line| {
+            points.push(line.start_point());
+            densify_between::<F, MetricSpace>(
+                line.start_point(),
+                line.end_point(),
+                &mut points,
+                max_segment_length,
+            )
+        });
+
+        // we're done, push the last coordinate on to finish
+        let final_coord = *self
+            .0
+            .last()
+            .expect("we already asserted the line string is not empty");
+        points.push(final_coord.into());
+
+        LineString::from(points)
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for MultiLineString<F> {
+    type Output = Self;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        MultiLineString::new(
+            self.iter()
+                .map(|line_string| line_string.densify::<MetricSpace>(max_segment_length))
+                .collect(),
+        )
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for Polygon<F> {
+    type Output = Self;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        Polygon::new(
+            self.exterior().densify::<MetricSpace>(max_segment_length),
+            self.interiors()
+                .iter()
+                .map(|interior| interior.densify::<MetricSpace>(max_segment_length))
+                .collect(),
+        )
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for MultiPolygon<F> {
+    type Output = Self;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        MultiPolygon::new(
+            self.iter()
+                .map(|polygon| polygon.densify::<MetricSpace>(max_segment_length))
+                .collect(),
+        )
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for Rect<F> {
+    type Output = Polygon<F>;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        self.to_polygon().densify::<MetricSpace>(max_segment_length)
+    }
+}
+
+impl<F: CoordFloat + FromPrimitive> Densify<F> for Triangle<F> {
+    type Output = Polygon<F>;
+
+    fn densify<MetricSpace>(&self, max_segment_length: F) -> Self::Output
+    where
+        MetricSpace: Distance<F, Point<F>, Point<F>> + InterpolatePoint<F>,
+    {
+        self.to_polygon().densify::<MetricSpace>(max_segment_length)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::{coord, polygon, wkt, Euclidean, Geodesic, Haversine, Rhumb};
+
+    #[test]
+    fn densify_line() {
+        // London to Paris
+        let line = Line::new(
+            coord!(x: -0.1278f64, y: 51.5074),
+            coord!(x: 2.3522, y: 48.8566),
+        );
+
+        let densified_line = line.densify::<Geodesic>(100_000.0); // max segment length 100km
+        assert!(densified_line.coords_count() > 2);
+
+        let densified_rhumb = line.densify::<Rhumb>(100_000.0);
+        assert!(densified_rhumb.coords_count() > 2);
+
+        let densified_haversine = line.densify::<Haversine>(100_000.0);
+        assert!(densified_haversine.coords_count() > 2);
+    }
+
+    #[test]
+    fn densify_line_string() {
+        let line_string = LineString::new(vec![
+            coord!(x: -58.3816f64, y: -34.6037), // Buenos Aires, Argentina
+            coord!(x: -77.0428, y: -12.0464),    // Lima, Peru
+            coord!(x: -47.9292, y: -15.7801),    // Brasília, Brazil
+        ]);
+
+        let densified_ls = line_string.densify::<Geodesic>(500_000.0); // 500 km max segment length
+        assert!(densified_ls.coords_count() > line_string.coords_count());
+
+        let densified_rhumb_ls = line_string.densify::<Rhumb>(500_000.0);
+        assert!(densified_rhumb_ls.coords_count() > line_string.coords_count());
+
+        let densified_haversine_ls = line_string.densify::<Haversine>(500_000.0);
+        assert!(densified_haversine_ls.coords_count() > line_string.coords_count());
+    }
+
+    #[test]
+    fn densify_polygon() {
+        let polygon = polygon![
+            (x: -58.3816f64, y: -34.6037), // Buenos Aires
+            (x: -77.0428, y: -12.0464),    // Lima
+            (x: -47.9292, y: -15.7801),    // Brasília
+        ];
+
+        let densified_polygon = polygon.densify::<Geodesic>(500_000.0); // 500 km max segment length
+        assert!(densified_polygon.exterior().coords_count() > polygon.exterior().coords_count());
+    }
+
+    // ported from the old Deprecated trait, which only worked with Euclidean measures
+    mod euclidean {
+        use super::*;
+
+        #[test]
+        fn test_polygon_densify() {
+            let polygon = wkt!(POLYGON(
+                (-5.0 0.0,0.0 5.0,5.0 0.0,-5.0 0.0),
+                (-3.0 0.0,0.0 3.0,3.0 0.0,-3.0 0.0)
+            ));
+
+            let expected = wkt!(POLYGON(
+                (-5.0 0.0,-3.75 1.25,-2.5 2.5,-1.25 3.75,0.0 5.0,1.25 3.75,2.5 2.5,3.75 1.25,5.0 0.0,3.0 0.0,1.0 0.0,-1.0000000000000009 0.0,-3.0 0.0, -5.0 0.0),
+                (-3.0 0.0,-2.0 1.0,-1.0 2.0,0.0 3.0,1.0 2.0,2.0 1.0,3.0 0.0,1.0 0.0,-1.0 0.0,-3.0 0.0)
+            ));
+
+            let max_dist = 2.0;
+            let densified = polygon.densify::<Euclidean>(max_dist);
+            assert_eq!(densified, expected);
+        }
+
+        #[test]
+        fn test_empty_linestring_densify() {
+            let linestring = LineString::<f64>::new(vec![]);
+            let max_dist = 2.0;
+            let densified = linestring.densify::<Euclidean>(max_dist);
+            assert!(densified.0.is_empty());
+        }
+
+        #[test]
+        fn test_linestring_densify() {
+            let linestring = wkt!(LINESTRING(
+               -1.0 0.0,
+                0.0 0.0,
+                0.0 6.0,
+                1.0 8.0
+            ));
+            let expected = wkt!(LINESTRING(
+               -1.0 0.0,
+                0.0 0.0,
+                0.0 2.0,
+                0.0 4.0,
+                0.0 6.0,
+                0.5 7.0,
+                1.0 8.0
+            ));
+            let max_dist = 2.0;
+            let densified = linestring.densify::<Euclidean>(max_dist);
+            assert_eq!(densified, expected);
+        }
+
+        #[test]
+        fn test_line_densify() {
+            let line: Line<f64> = Line::new(coord! {x: 0.0, y: 6.0}, coord! {x: 1.0, y: 8.0});
+            let correct: LineString<f64> = vec![[0.0, 6.0], [0.5, 7.0], [1.0, 8.0]].into();
+            let max_dist = 2.0;
+            let densified = line.densify::<Euclidean>(max_dist);
+            assert_eq!(densified, correct);
+        }
+    }
+
+    // ported from the now deprecated DensifyHaversine
+    mod lon_lat_tests {
+        use super::*;
+
+        #[test]
+        fn test_polygon_densify() {
+            let polygon = wkt!(POLYGON((
+                4.925 45.804,
+                4.732 45.941,
+                4.935 46.513,
+                5.821 46.103,
+                5.627 45.611,
+                5.355 45.883,
+                4.925 45.804
+            )));
+
+            let exepcted_haversine = wkt!(POLYGON((
+                4.925 45.804,
+                4.732 45.941,
+                4.8329711649985505 46.2270449096239,
+                4.935 46.513,
+                5.379659133344039 46.30885540136222,
+                5.821 46.103,
+                5.723570877658867 45.85704103535437,
+                5.627 45.611,
+                5.355 45.883,
+                4.925 45.804
+            )));
+
+            let actual_haversine = polygon.densify::<Haversine>(50000.0);
+            assert_relative_eq!(actual_haversine, exepcted_haversine);
+
+            let expected_geodesic = wkt!(POLYGON((
+                4.925 45.804,
+                4.732 45.941,
+                4.832972865149862 46.22705224065524,
+                4.935 46.513,
+                5.379653814979939 46.30886184400083,
+                5.821 46.103,
+                5.723572275808633 45.85704648840237,
+                5.627 45.611,
+                5.355 45.883,
+                4.925 45.804
+            )));
+            let actual_geodesic = polygon.densify::<Geodesic>(50000.0);
+            assert_relative_eq!(actual_geodesic, expected_geodesic);
+        }
+
+        #[test]
+        fn test_linestring_densify() {
+            let linestring = wkt!(LINESTRING(
+                -3.202 55.9471,
+                -3.2012 55.9476,
+                -3.1994 55.9476,
+                -3.1977 55.9481,
+                -3.196 55.9483,
+                -3.1947 55.9487,
+                -3.1944 55.9488,
+                -3.1944 55.949
+            ));
+
+            let expected = wkt!(LINESTRING(
+                -3.202 55.9471,
+                -3.2012 55.9476,
+                -3.2002999999999995 55.94760000327935,
+                -3.1994 55.9476,
+                -3.1985500054877773 55.94785000292509,
+                -3.1977 55.9481,
+                -3.196 55.9483,
+                -3.1947 55.9487,
+                -3.1944 55.9488,
+                -3.1944 55.949
+            ));
+
+            let dense = linestring.densify::<Haversine>(110.0);
+            assert_relative_eq!(dense, expected);
+        }
+
+        #[test]
+        fn test_line_densify() {
+            let output = wkt!(LINESTRING(0.0 0.0, 0.0 0.5, 0.0 1.0));
+            let line = Line::new(coord! {x: 0.0, y: 0.0}, coord! { x: 0.0, y: 1.0 });
+            let dense = line.densify::<Haversine>(100000.0);
+            assert_relative_eq!(dense, output);
+        }
+    }
+
+    mod degenerate {
+        use super::*;
+
+        #[test]
+        fn test_empty_linestring() {
+            let input = wkt!(LINESTRING EMPTY);
+            let dense = input.densify::<Euclidean>(1.0);
+            assert_eq!(0, dense.coords_count());
+            assert_eq!(input, dense);
+        }
+
+        #[test]
+        fn test_one_point_linestring() {
+            let input = wkt!(LINESTRING(1.0 1.0));
+            let dense = input.densify::<Euclidean>(1.0);
+            assert_eq!(1, dense.coords_count());
+            assert_eq!(input, dense);
+        }
+    }
+}
diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs
index 74dcb81fe..75e2c90f9 100644
--- a/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs
+++ b/geo/src/algorithm/line_measures/metric_spaces/euclidean.rs
@@ -1,5 +1,7 @@
-use super::super::Distance;
+use super::super::{Distance, InterpolatePoint};
+use crate::line_measures::densify::densify_between;
 use crate::{CoordFloat, Point};
+use num_traits::FromPrimitive;
 
 /// Operations on the [Euclidean plane] measure distance with the pythagorean formula -
 /// what you'd measure with a ruler.
@@ -53,6 +55,30 @@ impl<F: CoordFloat> Distance<F, Point<F>, Point<F>> for Euclidean {
     }
 }
 
+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;
+        start + diff * ratio_from_start
+    }
+
+    fn points_along_line(
+        start: Point<F>,
+        end: Point<F>,
+        max_distance: F,
+        include_ends: bool,
+    ) -> impl Iterator<Item = Point<F>> {
+        let mut container = vec![];
+        if include_ends {
+            container.push(start);
+        }
+        densify_between::<F, Self>(start, end, &mut container, max_distance);
+        if include_ends {
+            container.push(end);
+        }
+        container.into_iter()
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
diff --git a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs
index fbcfc7804..1821f28de 100644
--- a/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs
+++ b/geo/src/algorithm/line_measures/metric_spaces/geodesic.rs
@@ -4,7 +4,7 @@ use geographiclib_rs::{DirectGeodesic, InverseGeodesic};
 
 /// An ellipsoidal model of the earth, using methods given by [Karney (2013)].
 ///
-/// Distances are computed using [geodesic lines].
+/// Distances are computed using [geodesic lines] and are measured in meters.
 ///
 /// [geodesic lines]: https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid
 /// [Karney (2013)]:  https://arxiv.org/pdf/1109.4448.pdf
diff --git a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs
index 6813bbe77..dee062839 100644
--- a/geo/src/algorithm/line_measures/metric_spaces/haversine.rs
+++ b/geo/src/algorithm/line_measures/metric_spaces/haversine.rs
@@ -6,7 +6,7 @@ use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
 
 /// A spherical model of the earth using the [haversine formula].
 ///
-/// Distances are considered [great circle] lengths.
+/// Distances are considered [great circle] lengths and are measured in meters.
 ///
 /// # References
 ///
diff --git a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs
index 6c550311d..9a8b0fa94 100644
--- a/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs
+++ b/geo/src/algorithm/line_measures/metric_spaces/rhumb.rs
@@ -7,6 +7,8 @@ use crate::{CoordFloat, Point, MEAN_EARTH_RADIUS};
 /// Provides [rhumb line] (a.k.a. loxodrome) geometry operations. A rhumb line appears as a straight
 /// line on a Mercator projection map.
 ///
+/// Rhumb distance is measured in meters.
+///
 /// # References
 ///
 /// The distance, destination, and bearing implementations are adapted in part
diff --git a/geo/src/algorithm/line_measures/mod.rs b/geo/src/algorithm/line_measures/mod.rs
index 3c2bc41e2..cb9208489 100644
--- a/geo/src/algorithm/line_measures/mod.rs
+++ b/geo/src/algorithm/line_measures/mod.rs
@@ -15,5 +15,8 @@ pub use interpolate_point::InterpolatePoint;
 mod length;
 pub use length::Length;
 
+mod densify;
+pub use densify::Densify;
+
 pub mod metric_spaces;
 pub use metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb};
diff --git a/geo/src/algorithm/linestring_segment.rs b/geo/src/algorithm/linestring_segment.rs
index b361c1069..7ee5cd3e3 100644
--- a/geo/src/algorithm/linestring_segment.rs
+++ b/geo/src/algorithm/linestring_segment.rs
@@ -1,8 +1,6 @@
-use crate::line_interpolate_point::LineInterpolatePoint;
-use crate::{
-    Coord, Densify, DensifyHaversine, Euclidean, Haversine, Length, LineString, LinesIter,
-    MultiLineString,
-};
+use crate::algorithm::{Densify, Length, LineInterpolatePoint, LinesIter};
+use crate::geometry::{Coord, LineString, MultiLineString};
+use crate::line_measures::{Euclidean, Haversine};
 
 /// Segments a LineString into `segment_count` equal length LineStrings as a MultiLineString
 /// using Euclidean distance calculations.  See `LineStringSegmentizeHaversine`
@@ -47,7 +45,7 @@ pub trait LineStringSegmentizeHaversine {
 }
 
 macro_rules! implement_segmentize {
-    ($trait_name:ident, $method_name:ident, $metric_space:ty, $densify_method:ident) => {
+    ($trait_name:ident, $method_name:ident, $metric_space:ty) => {
         impl $trait_name for LineString {
             fn $method_name(&self, n: usize) -> Option<MultiLineString> {
                 if (n == usize::MIN) || (n == usize::MAX) {
@@ -62,7 +60,7 @@ macro_rules! implement_segmentize {
                 let mut cum_length = 0_f64;
                 let segment_prop = (1_f64) / (n as f64);
                 let segment_length = total_length * segment_prop;
-                let densified = self.$densify_method(segment_length - f64::EPSILON);
+                let densified = self.densify::<$metric_space>(segment_length - f64::EPSILON);
 
                 if densified.lines().count() == n {
                     let linestrings = densified
@@ -112,13 +110,11 @@ macro_rules! implement_segmentize {
     };
 }
 
-implement_segmentize!(LineStringSegmentize, line_segmentize, Euclidean, densify);
-
+implement_segmentize!(LineStringSegmentize, line_segmentize, Euclidean);
 implement_segmentize!(
     LineStringSegmentizeHaversine,
     line_segmentize_haversine,
-    Haversine,
-    densify_haversine
+    Haversine
 );
 
 #[cfg(test)]
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 9ced466c7..32e140c7a 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -66,12 +66,9 @@ pub use coordinate_position::CoordinatePosition;
 pub mod coords_iter;
 pub use coords_iter::CoordsIter;
 
-/// Densify linear geometry components
-pub mod densify;
-pub use densify::Densify;
-
 /// Densify spherical geometry components
 pub mod densify_haversine;
+#[allow(deprecated)]
 pub use densify_haversine::DensifyHaversine;
 
 /// Dimensionality of a geometry and its boundary, based on OGC-SFA.
@@ -190,7 +187,7 @@ pub use lines_iter::LinesIter;
 
 pub mod line_measures;
 pub use line_measures::metric_spaces::{Euclidean, Geodesic, Haversine, Rhumb};
-pub use line_measures::{Bearing, Destination, Distance, InterpolatePoint, Length};
+pub use line_measures::{Bearing, Densify, Destination, Distance, InterpolatePoint, Length};
 
 /// Split a LineString into n segments
 pub mod linestring_segment;