unit_sphere/vector/
intersection.rs

1// Copyright (c) 2024-2026 Ken Barker
2
3// Permission is hereby granted, free of charge, to any person obtaining a copy
4// of this software and associated documentation files (the "Software"),
5// to deal in the Software without restriction, including without limitation the
6// rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
7// sell copies of the Software, and to permit persons to whom the Software is
8// furnished to do so, subject to the following conditions:
9
10// The above copyright notice and this permission notice shall be included in
11// all copies or substantial portions of the Software.
12
13// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
18// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
19// THE SOFTWARE.
20
21//! The `intersection` module contains functions for calculating great-circle
22//! intersections using vectors.
23//!
24//! A pair of great circles intersect at two points unless they are coincident.\
25//! For example, points `u` and `v` in *Figure1*.
26//!
27//! ![great circle path](https://upload.wikimedia.org/wikipedia/commons/thumb/c/cb/Illustration_of_great-circle_distance.svg/220px-Illustration_of_great-circle_distance.svg.png)\
28//! *Figure 1 A pair of intersecting great circles*
29//!
30//! A great circle intersection point can simply be calculated by normalizing
31//! the [cross product](https://en.wikipedia.org/wiki/Cross_product) of their
32//! pole vectors.\
33//! If the resulting vector is too small to normalize, then the great circles
34//! are coincident, in which case they effectively *intersect* everywhere.
35//!
36//! If a pair of `Arc`s are on coincident great circles,
37//! the mormalized centroid of the arc midpoints is used instead of the
38//! intersection point.
39//!
40//! Otherwise `closest_intersection_point` calls `use_antipodal_point`to determine
41//! which intersection point is closer to the
42//! [centroid](https://en.wikipedia.org/wiki/Centroid) of the `Arc`s midpoints.
43
44use super::{Vector3, calculate_great_circle_atd, normalise, normalise_centroid, sq_distance};
45use angle_sc::{Angle, Radians};
46use num_traits::{Float, float::FloatConst};
47
48/// Calculate an intersection point between the poles of two Great Circles.
49/// See: <http://www.movable-type.co.uk/scripts/latlong-vectors.html#intersection>
50/// * `pole_0`, `pole_1` the poles.
51/// * `min_sq_value` the minimum square of a vector length to normalize.
52///
53/// return an intersection point or None if the poles represent coincident Great Circles.
54#[must_use]
55pub fn calculate_intersection<T>(
56    pole_0: &Vector3<T>,
57    pole_1: &Vector3<T>,
58    min_sq_value: T,
59) -> Option<Vector3<T>>
60where
61    T: Float + na::Scalar + na::ComplexField<RealField = T>,
62{
63    normalise(&pole_0.cross(pole_1), min_sq_value)
64}
65
66/// Determine whether the antipodal point is closer to the centroid of the
67/// `Arc`s.
68///
69/// * `point` a great-circle intersection point.
70/// * `centroid` the centroid of the `Arc`s mid points.
71///
72/// returns true if the antipodal intersection is closer to the `centroid`
73/// of the `Arc`s otherwise returns false.
74#[must_use]
75pub fn use_antipodal_point<T>(point: &Vector3<T>, centroid: &Vector3<T>) -> bool
76where
77    T: Float + na::Scalar + na::ComplexField<RealField = T>,
78{
79    sq_distance(centroid, &(-*point)) < sq_distance(centroid, point)
80}
81
82/// Return the closer intersection point to the centroid of the `Arc`s.
83///
84/// * `point` a great-circle intersection point.
85/// * `centroid` the centroid of the `Arc`s mid points.
86///
87/// returns the antipodal point if it is closer to the `centroid`,
88/// otherwise returns the point.
89#[must_use]
90pub fn closest_intersection_point<T>(point: &Vector3<T>, centroid: &Vector3<T>) -> Vector3<T>
91where
92    T: Float + na::Scalar + na::ComplexField<RealField = T>,
93{
94    if use_antipodal_point(point, centroid) {
95        -*point
96    } else {
97        *point
98    }
99}
100
101/// Determine the reference point of a pair of arcs.
102/// I.e. the closest intersection point if they intersect or the
103/// centroid normalized to lie on the unit sphere if they don't.
104///
105/// * `mid_point_0`, `mid_point_1` the mid points of the `Arc`s.
106/// * `pole_0`, `pole_1` the poles of the `Arc` great circles.
107/// * `sq_sin_max_coincident_angle` the square of the sine of the
108///   maximum angle between coincident great circles.
109///
110/// returns the closest intersection point or normalized centroid and the
111/// sine of the angle between the arcs, zero if the arcs are coincident.
112/// And the absolute relative angle at the intersection point or centroid.
113#[must_use]
114pub fn calculate_reference_point_and_angle<T>(
115    mid_point_0: &Vector3<T>,
116    pole_0: &Vector3<T>,
117    mid_point_1: &Vector3<T>,
118    pole_1: &Vector3<T>,
119    sq_sin_max_coincident_angle: T,
120) -> (Vector3<T>, Angle<T>)
121where
122    T: Float + na::Scalar + na::ComplexField<RealField = T>,
123    f64: From<T>,
124{
125    let centroid = mid_point_0 + mid_point_1;
126    // calculate the intersection point between the great circles
127    let point = pole_0.cross(pole_1);
128    normalise(&point, sq_sin_max_coincident_angle).map_or_else(
129        || {
130            // the great circles are coincident
131
132            let c = normalise_centroid(&centroid, mid_point_0, pole_0);
133            (c, Angle::from_y_x(point.norm(), pole_0.dot(pole_1)))
134        },
135        |p| {
136            // the great circles intersect
137
138            let x = closest_intersection_point(&p, &centroid);
139            (x, Angle::from_y_x(point.norm(), pole_0.dot(pole_1)))
140        },
141    )
142}
143
144/// Calculate signed great circle distances from two arc mid points to their
145/// closest intersection point or normalized centroid if the arcs are on coincident
146/// great circles.
147///
148/// * `mid_point_0`, `mid_point_1` the mid points of the arcs.
149/// * `pole_0`, `pole_1` the poles of the arc great circles.
150/// * `sq_sin_max_coincident_angle` the square of the sine of the
151///   maximum angle between coincident great circles.
152///
153/// returns the signed great circle distances of the closest intersection
154/// point or centroid  from the arc mid points in `Radians`,
155/// and the relative angle between the arc great circles.
156#[must_use]
157pub fn calculate_arc_reference_distances_and_angle<T>(
158    mid_point_0: &Vector3<T>,
159    pole_0: &Vector3<T>,
160    mid_point_1: &Vector3<T>,
161    pole_1: &Vector3<T>,
162    sq_sin_max_coincident_angle: T,
163) -> (Radians<T>, Radians<T>, Angle<T>)
164where
165    T: Float + FloatConst + na::Scalar + na::ComplexField<RealField = T>,
166    f64: From<T>,
167{
168    let (point, angle) = calculate_reference_point_and_angle(
169        mid_point_0,
170        pole_0,
171        mid_point_1,
172        pole_1,
173        sq_sin_max_coincident_angle,
174    );
175    let distance_0 = calculate_great_circle_atd(mid_point_0, pole_0, &point);
176    let distance_1 = calculate_great_circle_atd(mid_point_1, pole_1, &point);
177
178    (distance_0, distance_1, angle)
179}
180
181#[cfg(test)]
182mod tests {
183    use std::f64;
184
185    use super::*;
186    use crate::{LatLong, vector};
187    use angle_sc::{Angle, Degrees, is_within_tolerance};
188
189    pub const MIN_SQ_NORM: f64 = vector::tests::MIN_SQ_NORM;
190
191    #[test]
192    fn test_calculate_intersection() {
193        let lat_lon_south = LatLong::new(Degrees(-90.0), Degrees(0.0));
194        let south_pole = Vector3::from(&lat_lon_south);
195
196        let lat_lon_north = LatLong::new(Degrees(90.0), Degrees(0.0));
197        let north_pole = Vector3::from(&lat_lon_north);
198
199        let lat_lon_idl = LatLong::new(Degrees(0.0), Degrees(180.0));
200        let idl = Vector3::from(&lat_lon_idl);
201
202        let equator_intersection = calculate_intersection(&south_pole, &north_pole, MIN_SQ_NORM);
203        assert!(equator_intersection.is_none());
204
205        let gc_intersection1 = calculate_intersection(&idl, &north_pole, MIN_SQ_NORM).unwrap();
206        let gc_intersection2 = calculate_intersection(&idl, &south_pole, MIN_SQ_NORM).unwrap();
207
208        assert_eq!(gc_intersection1, -gc_intersection2);
209    }
210
211    #[test]
212    fn test_calculate_arc_reference_distances_and_angle_coincident_great_circles() {
213        let point_1 = Vector3::new(1.0, 0.0, 0.0);
214        let pole_1 = Vector3::new(0.0, 0.0, 1.0);
215
216        // same mid points and great circles
217        let result = calculate_arc_reference_distances_and_angle(
218            &point_1,
219            &pole_1,
220            &point_1,
221            &pole_1,
222            MIN_SQ_NORM,
223        );
224        assert_eq!(Radians(0.0), result.0);
225        assert_eq!(Radians(0.0), result.1);
226        assert_eq!(Degrees(0.0), Degrees::from(result.2));
227
228        // opposite mid points and same great circles
229        let point_m1 = -point_1;
230        let result = calculate_arc_reference_distances_and_angle(
231            &point_1,
232            &pole_1,
233            &point_m1,
234            &pole_1,
235            MIN_SQ_NORM,
236        );
237        assert!(is_within_tolerance(
238            -f64::consts::FRAC_PI_2,
239            result.0.0,
240            f64::EPSILON
241        ));
242        assert!(is_within_tolerance(
243            f64::consts::FRAC_PI_2,
244            result.1.0,
245            f64::EPSILON
246        ));
247        assert_eq!(Degrees(0.0), Degrees::from(result.2));
248
249        // opposite mid points and great circles
250        let pole_m1 = -pole_1;
251        let result = calculate_arc_reference_distances_and_angle(
252            &point_1,
253            &pole_1,
254            &point_m1,
255            &pole_m1,
256            MIN_SQ_NORM,
257        );
258        assert!(is_within_tolerance(
259            -f64::consts::FRAC_PI_2,
260            result.0.0,
261            f64::EPSILON
262        ));
263        assert!(is_within_tolerance(
264            -f64::consts::FRAC_PI_2,
265            result.1.0,
266            f64::EPSILON
267        ));
268        assert_eq!(Degrees(180.0), Degrees::from(result.2));
269    }
270
271    #[test]
272    fn test_calculate_arc_reference_distances_and_angle_intersecting_great_circles() {
273        let point_1 = Vector3::new(1.0, 0.0, 0.0);
274        let pole_1 = Vector3::new(0.0, 0.0, 1.0);
275        let pole_2 = Vector3::new(0.0, 1.0, 0.0);
276
277        // intersection, same mid points
278        let result = calculate_arc_reference_distances_and_angle(
279            &point_1,
280            &pole_1,
281            &point_1,
282            &pole_2,
283            MIN_SQ_NORM,
284        );
285        assert_eq!(Radians(0.0), result.0);
286        assert_eq!(Radians(0.0), result.1);
287        assert_eq!(Degrees(90.0), Degrees::from(result.2));
288
289        // intersection, same mid points, acute angle
290        let pole_3 = (pole_1 + pole_2).normalize();
291        let result = calculate_arc_reference_distances_and_angle(
292            &point_1,
293            &pole_1,
294            &point_1,
295            &pole_3,
296            MIN_SQ_NORM,
297        );
298        assert_eq!(Radians(0.0), result.0);
299        assert_eq!(Radians(0.0), result.1);
300        assert_eq!(Degrees(45.0), Degrees::from(result.2));
301
302        // intersection, same mid points, obtuse angle
303        let pole_m3 = -pole_3;
304        let result = calculate_arc_reference_distances_and_angle(
305            &point_1,
306            &pole_1,
307            &point_1,
308            &pole_m3,
309            MIN_SQ_NORM,
310        );
311        assert_eq!(Radians(0.0), result.0);
312        assert_eq!(Radians(0.0), result.1);
313        assert_eq!(Degrees(135.0), Degrees::from(result.2));
314
315        // intersection, different mid points, acute angle
316        let point_2 = vector::position(
317            &point_1,
318            &vector::direction(&point_1, &pole_3),
319            Angle::default().quarter_turn_cw(),
320        );
321        let result = calculate_arc_reference_distances_and_angle(
322            &point_1,
323            &pole_1,
324            &point_2,
325            &pole_3,
326            MIN_SQ_NORM,
327        );
328        assert_eq!(Radians(0.0), result.0);
329        assert!(is_within_tolerance(
330            -f64::consts::FRAC_PI_2,
331            result.1.0,
332            f64::EPSILON
333        ));
334        assert_eq!(Degrees(45.0), Degrees::from(result.2));
335
336        // intersection, different mid points, obtuse angle
337        let result = calculate_arc_reference_distances_and_angle(
338            &point_1,
339            &pole_1,
340            &point_2,
341            &pole_m3,
342            MIN_SQ_NORM,
343        );
344        assert_eq!(Radians(0.0), result.0);
345        assert!(is_within_tolerance(
346            f64::consts::FRAC_PI_2,
347            result.1.0,
348            f64::EPSILON
349        ));
350        assert_eq!(Degrees(135.0), Degrees::from(result.2));
351    }
352}
unit_sphere/vector/intersection.rs

unit_sphere/vector/
intersection.rs
