geo_core/
line_utils.rs

1use crate::{
2    constants::EARTH_RADIUS, latlng::LatLng, line::Line, point::Point, point3d::Point3D,
3    point_utils::point_to_point_distance_cosine, utils::deg_to_rad,
4};
5
6// PointToLineCartesianDistance - Get the min distance from a point to a line (Cartesian Coordinates)
7// https://en.m.wikipedia.org/wiki/Distance_from_a_point_to_a_line
8pub fn point_to_line_cartesian_distance(p: Point, l: Line) -> f64 {
9    let top = (((l.coords[1].y - l.coords[0].y) * p.x) - ((l.coords[1].x - l.coords[0].x) * p.y)
10        + (l.coords[1].x * l.coords[0].y)
11        - (l.coords[1].y * l.coords[0].x))
12        .abs();
13    let bottom =
14        (((l.coords[1].y - l.coords[1].y) * 2.0) + (l.coords[1].x - l.coords[0].x * 2.0)).sqrt();
15    let result = top / bottom;
16    return result;
17}
18
19// PointToLineGeographicDistance - Gets the min distnace from a point to a line (Geographic Cordinates)
20pub fn point_to_line_geographic_distance(p: LatLng, l: Line) -> f64 {
21    let a = LatLng {
22        lat: deg_to_rad(l.coords[0].y),
23        lng: deg_to_rad(l.coords[0].x),
24    };
25    let b = LatLng {
26        lat: deg_to_rad(l.coords[1].y),
27        lng: deg_to_rad(l.coords[1].x),
28    };
29    let c = p.convert_to_radian();
30
31    let t = nearest_point_great_circle(&a, &b, &c);
32    let a_point = a.convert_to_point();
33    let b_point = b.convert_to_point();
34    let c_point = c.convert_to_point();
35    let t_point = t.convert_to_point();
36
37    //If closest point is on the line use that.
38    if on_segment(a_point, b_point, t_point) {
39        return point_to_point_distance_cosine(t_point, c_point);
40    }
41
42    //Otherwise just use start or end point whichever is closer.
43    let distance_ac = point_to_point_distance_cosine(a_point, c_point);
44    let distance_bc = point_to_point_distance_cosine(b_point, c_point);
45
46    if distance_ac < distance_bc {
47        return distance_ac;
48    }
49    return distance_bc;
50}
51
52fn nearest_point_great_circle(a: &LatLng, b: &LatLng, c: &LatLng) -> LatLng {
53    let a_cartesian = a.convert_to_point3d();
54    let b_cartesian = b.convert_to_point3d();
55    let c_cartesian = c.convert_to_point3d();
56
57    let g = vector_product(&a_cartesian, &b_cartesian);
58    let f: Point3D = vector_product(&c_cartesian, &g);
59    let t = vector_product(&g, &f);
60    let norm = normalize(t);
61    let multi = multiply_by_scalar(norm, EARTH_RADIUS);
62    let cart = multi.convert_to_lat_lng();
63
64    return cart;
65}
66
67fn vector_product(a: &Point3D, b: &Point3D) -> Point3D {
68    let result = Point3D {
69        x: a.y * b.z - a.z * b.y,
70        y: a.z * b.x - a.x * b.z,
71        z: a.x * b.y - a.y * b.x,
72    };
73    return result;
74}
75
76fn normalize(t: Point3D) -> Point3D {
77    let length = ((t.x * t.x) + (t.y * t.y) + (t.z * t.z)).sqrt();
78    let result = Point3D {
79        x: t.x / length,
80        y: t.y / length,
81        z: t.z / length,
82    };
83    return result;
84}
85
86fn multiply_by_scalar(normalize: Point3D, k: f64) -> Point3D {
87    let result = Point3D {
88        x: normalize.x * k,
89        y: normalize.y * k,
90        z: normalize.z * k,
91    };
92    return result;
93}
94
95//Needs Radians -- Checks if point is on the line by substracting distances to check if total lenght - two sub lenghts are near enough to be zero.
96fn on_segment(a: Point, b: Point, t: Point) -> bool {
97    let diff = (point_to_point_distance_cosine(a, b)
98        - point_to_point_distance_cosine(a, t)
99        - point_to_point_distance_cosine(b, t))
100    .abs();
101    if diff < 0.1 {
102        return true;
103    }
104    return false;
105}
106
107#[cfg(test)]
108mod tests {
109
110    use super::*;
111
112    #[test]
113    fn test_point_to_line_geographic() {
114        let test_line = Line::new(
115            Point {
116                x: 1.070678,
117                y: 51.279336,
118            },
119            Point {
120                x: 1.071720,
121                y: 51.277559,
122            },
123        );
124
125        //Point which should return a distnace to a point on the line.
126        let start_point = LatLng {
127            lat: 51.27936,
128            lng: 1.07263,
129        };
130        let distance = point_to_line_geographic_distance(start_point, test_line);
131        assert_eq!(distance.round(), 128.0);
132
133        //Point which should return a distance to the line's start point.
134        let start_point2 = LatLng {
135            lat: 51.280345,
136            lng: 1.070411,
137        };
138        let distance2 = point_to_line_geographic_distance(start_point2, test_line);
139        assert_eq!(distance2.round(), 114.0);
140
141        //Point which should return a distance to the line's end point.
142        let start_point3 = LatLng {
143            lat: 51.277410,
144            lng: 1.072814,
145        };
146        let distance3 = point_to_line_geographic_distance(start_point3, test_line);
147
148        assert_eq!(distance3.round(), 78.0);
149    }
150}