gistools/geometry/tools/lines/
point_to_line_distance.rs

1use crate::{
2    geometry::{euclidean_distance, haversine_distance},
3    proj::Coords,
4};
5use libm::{fmax, fmin, hypot};
6use s2json::{GetXY, GetZ};
7
8/// The method to use to calculate the distance
9#[derive(Debug, Default, Copy, Clone, PartialEq)]
10pub enum DistanceMethod {
11    /// Euclidean
12    #[default]
13    Euclidean,
14    /// Haversine
15    Haversine,
16}
17
18#[derive(Debug, Default, Copy, Clone, PartialEq)]
19struct ClosestIndex {
20    index: usize,
21    dist: f64,
22}
23
24/// Check to see how far away the point is from the line. Supports both Euclidean and Haversine methods
25///
26/// ## Parameters
27/// - `line`: the line to check against
28/// - `point`: the point to check if it is on the line
29/// - `method`: the method to use, either 'euclidean' or 'haversine'. Defaults to [`DistanceMethod::Euclidean`]
30///
31/// ## Returns
32/// The shortest distance between the point and a line. Returns -1 if line is empty
33pub fn point_to_line_distance<P: GetXY + GetZ, Q: GetXY + GetZ>(
34    line: &[P],
35    point: &Q,
36    method: Option<DistanceMethod>,
37) -> f64 {
38    let method = method.unwrap_or(DistanceMethod::Euclidean);
39    let haversine = method == DistanceMethod::Haversine;
40
41    let mut closest_index: Option<ClosestIndex> = None;
42    for (i, line) in line.iter().enumerate() {
43        // get the distance between the point and the line's point at index
44        let dist = if haversine {
45            haversine_distance(point, line)
46        } else {
47            euclidean_distance(point, line)
48        };
49        if dist == 0. {
50            return 0.;
51        }
52        if closest_index.is_none() || dist < closest_index.unwrap().dist {
53            closest_index = Some(ClosestIndex { index: i, dist });
54        }
55    }
56
57    // If there is no closest point, return -1
58    if closest_index.is_none() {
59        return -1.;
60    }
61    let closest_index = closest_index.unwrap();
62
63    // If line is a single point, return distance to that point
64    if line.len() == 1 {
65        return closest_index.dist;
66    }
67
68    let curr = &line[closest_index.index];
69    // If the point is the start or end of the line, return distance to that point and next/prev
70    if closest_index.index == 0 {
71        return distance_point_to_segment(curr, &line[closest_index.index + 1], point, method);
72    }
73    if closest_index.index == line.len() - 1 {
74        return distance_point_to_segment(curr, &line[closest_index.index - 1], point, method);
75    }
76    let prev = &line[closest_index.index - 1];
77    let next = &line[closest_index.index + 1];
78
79    // Check against both sides of the line's closest point
80    let dist1 = distance_point_to_segment(curr, prev, point, method);
81    let dist2 = distance_point_to_segment(curr, next, point, method);
82
83    if dist1 < dist2 { dist1 } else { dist2 }
84}
85
86/// Get the distance between a point and a segment
87///
88/// @param a - the segment start point
89/// @param b - the segment end point
90/// @param p - the point to check
91/// @param method - the method to use, either 'euclidean' or 'haversine'. Defaults to 'euclidean'
92/// @returns - the distance
93fn distance_point_to_segment<A: GetXY, B: GetXY, P: GetXY>(
94    a: &A,
95    b: &B,
96    p: &P,
97    method: DistanceMethod,
98) -> f64 {
99    if method == DistanceMethod::Haversine {
100        // approximate by sampling along the great-circle segment
101        // project p onto AB using Euclidean math in lat/lon degrees
102        // but compute distances with haversine_distance
103        let abx = b.x() - a.x();
104        let aby = b.y() - a.y();
105        let apx = p.x() - a.x();
106        let apy = p.y() - a.y();
107        let ab_len_sq = abx * abx + aby * aby;
108
109        if ab_len_sq == 0. {
110            return haversine_distance(p, a);
111        }
112
113        let mut t = (apx * abx + apy * aby) / ab_len_sq;
114        t = fmax(0., fmin(1., t));
115
116        let proj = Coords::new_xy(a.x() + t * abx, a.y() + t * aby);
117        return haversine_distance(p, &proj);
118    }
119
120    // Euclidean fallback
121    let ab = Coords::new_xy(b.x() - a.x(), b.y() - a.y());
122    let ap = Coords::new_xy(p.x() - a.x(), p.y() - a.y());
123    let ab_len_sq = ab.x() * ab.x() + ab.y() * ab.y();
124
125    if ab_len_sq == 0. {
126        return hypot(p.x() - a.x(), p.y() - a.y());
127    }
128
129    let mut t = (ap.x() * ab.x() + ap.y() * ab.y()) / ab_len_sq;
130    t = fmax(0., fmin(1., t));
131
132    let closest = Coords::new_xy(a.x() + t * ab.x(), a.y() + t * ab.y());
133    hypot(p.x() - closest.x(), p.y() - closest.y())
134}
gistools/geometry/tools/lines/point_to_line_distance.rs

gistools/geometry/tools/lines/
point_to_line_distance.rs
