Optimise line locate point

geo-types/src/private_utils.rs

@@ -61,6 +61,14 @@ fn get_min_max<T: PartialOrd>(p: T, min: T, max: T) -> (T, T) {
 }
 
 pub fn line_segment_distance<T, C>(point: C, start: C, end: C) -> T
+where
+    T: CoordFloat,
+    C: Into<Coord<T>>,
+{
+    line_segment_distance_squared(point, start, end).sqrt()
+}
+
+pub fn line_segment_distance_squared<T, C>(point: C, start: C, end: C) -> T
 where
     T: CoordFloat,
     C: Into<Coord<T>>,
@@ -69,28 +77,41 @@ where
     let start = start.into();
     let end = end.into();
 
+    // Degenerate case for line with length 0 - treat as a point
     if start == end {
-        return line_euclidean_length(Line::new(point, start));
+        return line_euclidean_length_squared(Line::new(point, start));
     }
     let dx = end.x - start.x;
     let dy = end.y - start.y;
     let d_squared = dx * dx + dy * dy;
+
+    // Projection of point onto the line segment
     let r = ((point.x - start.x) * dx + (point.y - start.y) * dy) / d_squared;
+    // Projection lies beyond start - start point is closest
     if r <= T::zero() {
-        return line_euclidean_length(Line::new(point, start));
+        return line_euclidean_length_squared(Line::new(point, start));
     }
+    // Projection lies beyond end - end point is closest
     if r >= T::one() {
-        return line_euclidean_length(Line::new(point, end));
+        return line_euclidean_length_squared(Line::new(point, end));
     }
+    // Projection lies on midpoint between start-end
     let s = ((start.y - point.y) * dx - (start.x - point.x) * dy) / d_squared;
-    s.abs() * dx.hypot(dy)
+    s.powi(2) * d_squared
 }
 
 pub fn line_euclidean_length<T>(line: Line<T>) -> T
 where
     T: CoordFloat,
 {
-    line.dx().hypot(line.dy())
+    line_euclidean_length_squared(line).sqrt()
+}
+
+pub fn line_euclidean_length_squared<T>(line: Line<T>) -> T
+where
+    T: CoordFloat,
+{
+    line.dx().powi(2) + line.dy().powi(2)
 }
 
 pub fn point_line_string_euclidean_distance<T>(p: Point<T>, l: &LineString<T>) -> T
@@ -114,6 +135,14 @@ where
     line_segment_distance(p.into(), l.start, l.end)
 }
 
+pub fn point_line_euclidean_distance_squared<C, T>(p: C, l: Line<T>) -> T
+where
+    T: CoordFloat,
+    C: Into<Coord<T>>,
+{
+    line_segment_distance_squared(p.into(), l.start, l.end)
+}
+
 pub fn point_contains_point<T>(p1: Point<T>, p2: Point<T>) -> bool
 where
     T: CoordFloat,

geo/Cargo.toml

@@ -126,3 +126,7 @@ harness = false
 [[bench]]
 name = "stitch"
 harness = false
+
+[[bench]]
+name = "line_locate_point"
+harness = false

geo/benches/line_locate_point.rs

@@ -0,0 +1,21 @@
+use criterion::{criterion_group, criterion_main, Criterion};
+use geo::{InterpolatableLine, Euclidean, LineLocatePoint};
+use geo::LineString;
+
+fn criterion_benchmark(c: &mut Criterion) {
+    let test_line: LineString = geo_test_fixtures::norway_main();
+
+    let percentage = 0.45;
+
+    c.bench_function("LineString line_locate_point", |bencher| {
+        bencher.iter(|| {
+            let point = test_line.point_at_ratio_from_start(&Euclidean, percentage).expect("Did not interpolate out point");
+            let percentage_calced = criterion::black_box(test_line.line_locate_point(&point)).expect("Didn't interpolate point back to percentage");
+
+            approx::assert_abs_diff_eq!(percentage, percentage_calced, epsilon = 0.0000001);
+        });
+    });
+}
+
+criterion_group!(benches, criterion_benchmark);
+criterion_main!(benches);

geo/src/algorithm/line_locate_point.rs

@@ -8,6 +8,8 @@ use crate::{
 };
 use std::ops::AddAssign;
 
+use super::{line_measures::ComparableDistance, Euclidean};
+
 /// Returns a (option of the) fraction of the line's total length
 /// representing the location of the closest point on the line to
 /// the given point.
@@ -98,7 +100,8 @@ where
         let mut closest_dist_to_point = T::infinity();
         let mut fraction = T::zero();
         for segment in self.lines() {
-            let segment_distance_to_point = segment.euclidean_distance(p);
+            let segment_distance_to_point = geo_types::private_utils::point_line_euclidean_distance_squared(p.clone(), segment.clone());
+            // let segment_distance_to_point = Euclidean::comparable_distance(&Euclidean, p, &segment);
             let segment_length = segment.euclidean_length();
             let segment_fraction = segment.line_locate_point(p)?; // if any segment has a None fraction, return None
             if segment_distance_to_point < closest_dist_to_point {

geo/src/algorithm/line_measures/comparable_distance.rs

@@ -0,0 +1,34 @@
+/// Calculate a minimum distance between two geometries in a way that is useful for sorting operations
+pub trait ComparableDistance<F, Origin, Destination> {
+    /// This trait differs from [Distance](geo::Distance) in that the value returned is not a true distance 
+    /// but a value that is indicative of the true distance that can be used for sorting. It is generally
+    /// faster to compute this value rather than true distance if all you need is a comparable figure between two 
+    /// geometries 
+    /// 
+    /// Note that not all implementations support all geometry combinations, but at least `Point` to `Point`
+    /// is supported.
+    /// See [specific implementations](#implementors) for details.
+    ///
+    /// # Units
+    ///
+    /// - `origin`, `destination`: geometry where the units of x/y depend on the trait implementation.
+    /// - returns: depends on the trait implementation.
+    ///
+    /// # Examples
+    ///
+    /// ```
+    /// use geo::{Euclidean, ComparableDistance, Point};
+    /// let p1: Point = Point::new(0.0, 0.0);
+    /// let p2: Point = Point::new(0.0, 2.0);
+    /// let p3: Point = Point::new(0.0, 5.0);
+    /// 
+    /// let comparable_1_2 = Euclidean.comparable_distance(p1, p2);
+    /// let comparable_1_3 = Euclidean.comparable_distance(p1, p3);
+    ///
+    /// assert_eq!(comparable_1_2, 4.0);
+    /// assert_eq!(comparable_1_3, 25.0);
+    /// assert_lt!(comparable_1_2, comparable_1_3);
+    ///
+    /// ```
+    fn comparable_distance(&self, origin: Origin, destination: Destination) -> F;
+}

geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs

@@ -1,4 +1,4 @@
-use super::{Distance, Euclidean};
+use super::{ComparableDistance, Distance, Euclidean};
 use crate::algorithm::Intersects;
 use crate::coordinate_position::{coord_pos_relative_to_ring, CoordPos};
 use crate::geometry::*;
@@ -27,13 +27,23 @@ macro_rules! symmetric_distance_impl {
 
 impl<F: CoordFloat> Distance<F, Coord<F>, Coord<F>> for Euclidean {
     fn distance(&self, origin: Coord<F>, destination: Coord<F>) -> F {
+        Self.comparable_distance(origin, destination).sqrt()
+    }
+}
+impl<F: CoordFloat> ComparableDistance<F, Coord<F>, Coord<F>> for Euclidean {
+    fn comparable_distance(&self, origin: Coord<F>, destination: Coord<F>) -> F {
         let delta = origin - destination;
-        delta.x.hypot(delta.y)
+        delta.x.powi(2) + delta.y.powi(2)
     }
 }
 impl<F: CoordFloat> Distance<F, Coord<F>, &Line<F>> for Euclidean {
     fn distance(&self, coord: Coord<F>, line: &Line<F>) -> F {
-        ::geo_types::private_utils::point_line_euclidean_distance(Point(coord), *line)
+        Self::comparable_distance(&self, coord, line).sqrt()
+    }
+}
+impl<F: CoordFloat> ComparableDistance<F, Coord<F>, &Line<F>> for Euclidean {
+    fn comparable_distance(&self, coord: Coord<F>, line: &Line<F>) -> F {
+        ::geo_types::private_utils::point_line_euclidean_distance_squared(Point(coord), *line)
     }
 }
 
@@ -75,15 +85,33 @@ impl<F: CoordFloat> Distance<F, Point<F>, Point<F>> for Euclidean {
     }
 }
 
+impl<F: CoordFloat> ComparableDistance<F, Point<F>, Point<F>> for Euclidean {
+    fn comparable_distance(&self, origin: Point<F>, destination: Point<F>) -> F {
+        self.comparable_distance(origin.0, destination.0)
+    }
+}
+
 impl<F: CoordFloat> Distance<F, &Point<F>, &Point<F>> for Euclidean {
     fn distance(&self, origin: &Point<F>, destination: &Point<F>) -> F {
         self.distance(*origin, *destination)
     }
 }
 
+impl<F: CoordFloat> ComparableDistance<F, &Point<F>, &Point<F>> for Euclidean {
+    fn comparable_distance(&self, origin: &Point<F>, destination: &Point<F>) -> F {
+        self.comparable_distance(*origin, *destination)
+    }
+}
+
 impl<F: CoordFloat> Distance<F, &Point<F>, &Line<F>> for Euclidean {
     fn distance(&self, origin: &Point<F>, destination: &Line<F>) -> F {
-        geo_types::private_utils::point_line_euclidean_distance(*origin, *destination)
+        self.comparable_distance(origin, destination).sqrt()
+    }
+}
+
+impl<F: CoordFloat> ComparableDistance<F, &Point<F>, &Line<F>> for Euclidean {
+    fn comparable_distance(&self, origin: &Point<F>, destination: &Line<F>) -> F {
+        geo_types::private_utils::point_line_euclidean_distance_squared(*origin, *destination)
     }
 }
 

geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs

@@ -1,6 +1,6 @@
 mod distance;
 
-use super::super::{Distance, InterpolatePoint};
+use super::super::{ComparableDistance, Distance, InterpolatePoint};
 use crate::line_measures::densify::densify_between;
 use crate::{CoordFloat, Point};
 use num_traits::FromPrimitive;

geo/src/algorithm/line_measures/mod.rs

@@ -22,6 +22,9 @@
 mod bearing;
 pub use bearing::Bearing;
 
+mod comparable_distance;
+pub use comparable_distance::ComparableDistance;
+
 mod destination;
 pub use destination::Destination;