diff --git a/geo/CHANGES.md b/geo/CHANGES.md
index fe3543e2c..7f6a1e8c8 100644
--- a/geo/CHANGES.md
+++ b/geo/CHANGES.md
@@ -1,5 +1,9 @@
 # Changes
 
+## Unreleased
+* Add densification algorithm for linear geometry components
+  * <https://github.com/georust/geo/pull/847>
+
 ## 0.21.0
 
 * Boolean operations for `Polygon`s and `MultiPolygon`s: intersect, union, xor,
diff --git a/geo/src/algorithm/densify.rs b/geo/src/algorithm/densify.rs
new file mode 100644
index 000000000..8bfd422b0
--- /dev/null
+++ b/geo/src/algorithm/densify.rs
@@ -0,0 +1,237 @@
+use crate::{
+    CoordFloat, EuclideanLength, Line, LineInterpolatePoint, LineString, MultiLineString,
+    MultiPolygon, Point, Polygon, Rect, Triangle,
+};
+
+/// 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
+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);
+    }
+}
+
+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(),
+        )
+    }
+}
+
+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)
+    }
+}
+
+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(),
+        )
+    }
+}
+
+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 {
+        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)
+    }
+}
+
+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)
+    }
+}
+
+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)
+    }
+}
+
+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, Coordinate};
+
+    #[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![
+            Coordinate { x: -5.0, y: 0.0 },
+            Coordinate { x: -3.75, y: 1.25 },
+            Coordinate { x: -2.5, y: 2.5 },
+            Coordinate { x: -1.25, y: 3.75 },
+            Coordinate { x: 0.0, y: 5.0 },
+            Coordinate { x: 1.25, y: 3.75 },
+            Coordinate { x: 2.5, y: 2.5 },
+            Coordinate { x: 3.75, y: 1.25 },
+            Coordinate { x: 5.0, y: 0.0 },
+            Coordinate { x: 3.0, y: 0.0 },
+            Coordinate { x: 1.0, y: 0.0 },
+            Coordinate {
+                x: -1.0000000000000009,
+                y: 0.0,
+            },
+            Coordinate { x: -3.0, y: 0.0 },
+            Coordinate { x: -5.0, y: 0.0 },
+        ]);
+        let correct_int: LineString<f64> = LineString(vec![
+            Coordinate { x: -3.0, y: 0.0 },
+            Coordinate { x: -2.0, y: 1.0 },
+            Coordinate { x: -1.0, y: 2.0 },
+            Coordinate { x: 0.0, y: 3.0 },
+            Coordinate { x: 1.0, y: 2.0 },
+            Coordinate { x: 2.0, y: 1.0 },
+            Coordinate { x: 3.0, y: 0.0 },
+            Coordinate { x: 1.0, y: 0.0 },
+            Coordinate { x: -1.0, y: 0.0 },
+            Coordinate { 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_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/line_interpolate_point.rs b/geo/src/algorithm/line_interpolate_point.rs
index 8f9cebc1c..fdf2cdd2a 100644
--- a/geo/src/algorithm/line_interpolate_point.rs
+++ b/geo/src/algorithm/line_interpolate_point.rs
@@ -1,7 +1,6 @@
 use crate::coords_iter::CoordsIter;
-use std::ops::AddAssign;
-
 use crate::{CoordFloat, EuclideanLength, Line, LineString, Point};
+use std::ops::AddAssign;
 
 /// Returns an option of the point that lies a given fraction along the line.
 ///
@@ -109,6 +108,7 @@ where
 
 #[cfg(test)]
 mod test {
+
     use super::*;
     use crate::{coord, point};
     use crate::{ClosestPoint, LineLocatePoint};
diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs
index 66c08f781..3bfd7ed55 100644
--- a/geo/src/algorithm/mod.rs
+++ b/geo/src/algorithm/mod.rs
@@ -191,3 +191,7 @@ pub mod sweep;
 /// Boolean Ops;
 pub mod bool_ops;
 pub use bool_ops::{BooleanOps, OpType};
+
+/// Densify linear geometry components
+pub mod densify;
+pub use densify::Densify;
diff --git a/geo/src/lib.rs b/geo/src/lib.rs
index 3803d262d..c15702934 100644
--- a/geo/src/lib.rs
+++ b/geo/src/lib.rs
@@ -136,6 +136,7 @@
 //! - **[`HaversineIntermediate`](HaversineIntermediate)**:
 //! - **[`Proj`](Proj)**: Project geometries with the `proj` crate (requires the `use-proj` feature)
 //! - **[`ChaikinSmoothing`](ChaikinSmoothing)**: Smoothen `LineString`, `Polygon`, `MultiLineString` and `MultiPolygon` using Chaikins algorithm.
+//! - **[`Densify`](Densify)**: Densify linear geometry components by interpolating points
 //!
 //! # Features
 //!