diff --git a/geo-types/src/private_utils.rs b/geo-types/src/private_utils.rs index 001fcffbb8..67e619ae15 100644 --- a/geo-types/src/private_utils.rs +++ b/geo-types/src/private_utils.rs @@ -61,6 +61,14 @@ fn get_min_max(p: T, min: T, max: T) -> (T, T) { } pub fn line_segment_distance(point: C, start: C, end: C) -> T +where + T: CoordFloat, + C: Into>, +{ + line_segment_distance_squared(point, start, end).sqrt() +} + +pub fn line_segment_distance_squared(point: C, start: C, end: C) -> T where T: CoordFloat, C: Into>, @@ -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(line: Line) -> T where T: CoordFloat, { - line.dx().hypot(line.dy()) + line_euclidean_length_squared(line).sqrt() +} + +pub fn line_euclidean_length_squared(line: Line) -> T +where + T: CoordFloat, +{ + line.dx().powi(2) + line.dy().powi(2) } pub fn point_line_string_euclidean_distance(p: Point, l: &LineString) -> T @@ -114,6 +135,14 @@ where line_segment_distance(p.into(), l.start, l.end) } +pub fn point_line_euclidean_distance_squared(p: C, l: Line) -> T +where + T: CoordFloat, + C: Into>, +{ + line_segment_distance_squared(p.into(), l.start, l.end) +} + pub fn point_contains_point(p1: Point, p2: Point) -> bool where T: CoordFloat, diff --git a/geo/Cargo.toml b/geo/Cargo.toml index 814b119790..9e101687c6 100644 --- a/geo/Cargo.toml +++ b/geo/Cargo.toml @@ -126,3 +126,7 @@ harness = false [[bench]] name = "stitch" harness = false + +[[bench]] +name = "line_locate_point" +harness = false diff --git a/geo/benches/line_locate_point.rs b/geo/benches/line_locate_point.rs new file mode 100644 index 0000000000..5baaa2b5e7 --- /dev/null +++ b/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); \ No newline at end of file diff --git a/geo/src/algorithm/line_locate_point.rs b/geo/src/algorithm/line_locate_point.rs index 5d34c2d940..4f6a699236 100644 --- a/geo/src/algorithm/line_locate_point.rs +++ b/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 { diff --git a/geo/src/algorithm/line_measures/comparable_distance.rs b/geo/src/algorithm/line_measures/comparable_distance.rs new file mode 100644 index 0000000000..f45d22ba38 --- /dev/null +++ b/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 { + /// 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; +} diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs index 868fff6679..648e2a2048 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs +++ b/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 Distance, Coord> for Euclidean { fn distance(&self, origin: Coord, destination: Coord) -> F { + Self.comparable_distance(origin, destination).sqrt() + } +} +impl ComparableDistance, Coord> for Euclidean { + fn comparable_distance(&self, origin: Coord, destination: Coord) -> F { let delta = origin - destination; - delta.x.hypot(delta.y) + delta.x.powi(2) + delta.y.powi(2) } } impl Distance, &Line> for Euclidean { fn distance(&self, coord: Coord, line: &Line) -> F { - ::geo_types::private_utils::point_line_euclidean_distance(Point(coord), *line) + Self::comparable_distance(&self, coord, line).sqrt() + } +} +impl ComparableDistance, &Line> for Euclidean { + fn comparable_distance(&self, coord: Coord, line: &Line) -> F { + ::geo_types::private_utils::point_line_euclidean_distance_squared(Point(coord), *line) } } @@ -75,15 +85,33 @@ impl Distance, Point> for Euclidean { } } +impl ComparableDistance, Point> for Euclidean { + fn comparable_distance(&self, origin: Point, destination: Point) -> F { + self.comparable_distance(origin.0, destination.0) + } +} + impl Distance, &Point> for Euclidean { fn distance(&self, origin: &Point, destination: &Point) -> F { self.distance(*origin, *destination) } } +impl ComparableDistance, &Point> for Euclidean { + fn comparable_distance(&self, origin: &Point, destination: &Point) -> F { + self.comparable_distance(*origin, *destination) + } +} + impl Distance, &Line> for Euclidean { fn distance(&self, origin: &Point, destination: &Line) -> F { - geo_types::private_utils::point_line_euclidean_distance(*origin, *destination) + self.comparable_distance(origin, destination).sqrt() + } +} + +impl ComparableDistance, &Line> for Euclidean { + fn comparable_distance(&self, origin: &Point, destination: &Line) -> F { + geo_types::private_utils::point_line_euclidean_distance_squared(*origin, *destination) } } diff --git a/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs b/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs index 4312a308b5..d6e441bc6d 100644 --- a/geo/src/algorithm/line_measures/metric_spaces/euclidean/mod.rs +++ b/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; diff --git a/geo/src/algorithm/line_measures/mod.rs b/geo/src/algorithm/line_measures/mod.rs index ff41a45f39..2727225c48 100644 --- a/geo/src/algorithm/line_measures/mod.rs +++ b/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;