unit_sphere/vector/intersection.rs
1// Copyright (c) 2024-2025 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//! \
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//! `calculate_coincident_arc_distances` calculates the distances between
38//! `Arc` ends, zero if the `Arc`s overlap.
39//!
40//! Otherwise `use_antipodal_point` determines which intersection point
41//! is closer to the [centroid](https://en.wikipedia.org/wiki/Centroid)
42//! of the `Arc`s midpoints.
43//! `calculate_intersection_distances` then calculates great-circle distances
44//! along the `Arc`s to the intersection point.
45
46use super::{
47 calculate_great_circle_atd, normalise, sq_distance, Vector3d, MIN_SQ_DISTANCE, MIN_SQ_NORM,
48};
49use angle_sc::{max, Radians};
50
51/// Calculate an intersection point between the poles of two Great Circles.
52/// See: <http://www.movable-type.co.uk/scripts/latlong-vectors.html#intersection>
53/// * `pole1`, `pole2` the poles.
54/// * `min_sq_value` the minimum square of a vector length to normalize.
55///
56/// return an intersection point or None if the poles represent coincident Great Circles.
57#[must_use]
58pub fn calculate_intersection(
59 pole1: &Vector3d,
60 pole2: &Vector3d,
61 min_sq_value: f64,
62) -> Option<Vector3d> {
63 normalise(&pole1.cross(pole2), min_sq_value)
64}
65
66/// Calculate the great circle distances to an intersection point from the
67/// start points of a pair of great circle arcs, on different great circles.
68///
69/// * `a1`, `a2` the start points of the great circle arcs
70/// * `pole1`, `pole2` the poles of the great circle arcs
71/// * `c` the intersection point
72///
73/// returns a pair of great circle distances along the arcs to the
74/// intersection point in `Radians`.
75#[must_use]
76pub fn calculate_intersection_distances(
77 a1: &Vector3d,
78 pole1: &Vector3d,
79 a2: &Vector3d,
80 pole2: &Vector3d,
81 c: &Vector3d,
82) -> (Radians, Radians) {
83 (
84 calculate_great_circle_atd(a1, pole1, c),
85 calculate_great_circle_atd(a2, pole2, c),
86 )
87}
88
89/// Whether an intersection point is within an `Arc`.
90///
91/// * `distance` - the along track distance to the point from the start of the `Arc`.
92/// * `length` the length of the `Arc`.
93///
94/// return true if the intersection point is within the `Arc`, false otherwise.
95#[must_use]
96pub fn is_within(distance: f64, length: f64) -> bool {
97 (-f64::EPSILON <= distance) && (distance <= length + (f64::EPSILON * (1.0 + length)))
98}
99
100/// Calculate the great-circle distances along a pair of `Arc`s on coincident
101/// Great Circles to their closest (reference) points.
102///
103/// * `gc_d` the great-circle distance between the arc start points.
104/// * `reciprocal` whether the arcs are in reciprocal directions.
105/// * `arc1_length`, `arc2_length` the `Arc` lengths in `Radians`.
106///
107/// returns the distances along the first `Arc` and second `Arc` to their closest
108/// (reference) points in `Radians`.
109#[must_use]
110pub fn calculate_coincident_arc_distances(
111 gc_d: Radians,
112 reciprocal: bool,
113 arc1_length: Radians,
114 arc2_length: Radians,
115) -> (Radians, Radians) {
116 if reciprocal {
117 // if the arcs intersect
118 if is_within(gc_d.0, max(arc1_length, arc2_length).0) {
119 if gc_d <= arc2_length {
120 // The start of the first `Arc` is within the second `Arc`
121 (Radians(0.0), gc_d)
122 } else {
123 // The start of the second `Arc` is within the first `Arc`
124 (gc_d, Radians(0.0))
125 }
126 } else {
127 let abs_d = gc_d.abs();
128
129 // The distance between the `Arc` b ends
130 let b_d = abs_d.0 - arc1_length.0 - arc2_length.0;
131 // The distance between the `Arc` b ends around the Great Circle
132 let b_gc_d = if Radians(0.0) < gc_d {
133 b_d
134 } else {
135 core::f64::consts::TAU - b_d
136 };
137 if b_gc_d < abs_d.0 {
138 // The end of the second `Arc` is beyond the end of first `Arc`
139 (Radians(b_gc_d) + arc1_length, arc2_length)
140 } else {
141 // The start of the second `Arc` is before the start of first `Arc`
142 (-abs_d, Radians(0.0))
143 }
144 }
145 } else {
146 // The distance to the start of arc2 from the end of arc1
147 let b1a2 = if Radians(0.0) < gc_d {
148 gc_d.0 - arc1_length.0
149 } else {
150 core::f64::consts::TAU + gc_d.0 - arc1_length.0
151 };
152 // The distance to the start of arc1 from the end of arc2
153 let b2a1 = if Radians(0.0) < gc_d {
154 core::f64::consts::TAU - gc_d.0 - arc2_length.0
155 } else {
156 -gc_d.0 - arc2_length.0
157 };
158 if b2a1 < b1a2 {
159 // The start of the first arc is within the second arc
160 (Radians(0.0), Radians(b2a1 + arc2_length.0))
161 } else {
162 // The start of the second arc relative to the start of first arc.
163 (Radians(b1a2 + arc1_length.0), Radians(0.0))
164 }
165 }
166}
167
168/// Determine whether the antipodal point is closer to the centroid of the
169/// `Arc`s.
170///
171/// * `point` a great-circle intersection point.
172/// * `centroid` the centroid (geometric mean) of the `Arc`s mid points.
173///
174/// returns true if the antipodal intersection is closer to the `centroid`
175/// of the `Arc`s otherwise returns false.
176#[must_use]
177pub fn use_antipodal_point(point: &Vector3d, centroid: &Vector3d) -> bool {
178 sq_distance(centroid, &(-*point)) < sq_distance(centroid, point)
179}
180
181/// Calculate the great-circle distances along a pair of arcs to their
182/// closest intersection point or their coincident arc distances if the
183/// `Arc`s are on coincident Great Circles.
184///
185/// * `a1`, `a2` the `Arc` start points.
186/// * `pole1`, `pole1` the `Arc` poles.
187/// * `length1`, `length2` the `Arc` lengths.
188/// * `centroid` the centroid (geometric mean) of the `Arc`s mid points.
189///
190/// returns the distances along the first arc and second arc to the intersection
191/// point or to their coincident arc distances if the arcs do not intersect.
192#[must_use]
193pub fn calculate_intersection_point_distances(
194 a1: &Vector3d,
195 pole1: &Vector3d,
196 length1: Radians,
197 a2: &Vector3d,
198 pole2: &Vector3d,
199 length2: Radians,
200 centroid: &Vector3d,
201) -> (Radians, Radians) {
202 // Calculate the square of the Euclidean distance between the start points.
203 let sq_d = sq_distance(a1, a2);
204 if sq_d < MIN_SQ_DISTANCE {
205 (Radians(0.0), Radians(0.0))
206 } else {
207 calculate_intersection(pole1, pole2, MIN_SQ_NORM).map_or_else(
208 || {
209 calculate_coincident_arc_distances(
210 calculate_great_circle_atd(a1, pole1, a2),
211 pole1.dot(pole2) < 0.0,
212 length1,
213 length2,
214 )
215 },
216 |c| {
217 // Find the closest intersection point
218 let c = if use_antipodal_point(&c, centroid) {
219 -c
220 } else {
221 c
222 };
223 calculate_intersection_distances(a1, pole1, a2, pole2, &c)
224 },
225 )
226 }
227}
228
229#[cfg(test)]
230mod tests {
231 use super::*;
232 use crate::{vector, LatLong};
233 use angle_sc::{is_within_tolerance, Angle, Degrees};
234
235 #[test]
236 fn test_calculate_intersection() {
237 let lat_lon_south = LatLong::new(Degrees(-90.0), Degrees(0.0));
238 let south_pole = Vector3d::from(&lat_lon_south);
239
240 let lat_lon_north = LatLong::new(Degrees(90.0), Degrees(0.0));
241 let north_pole = Vector3d::from(&lat_lon_north);
242
243 let lat_lon_idl = LatLong::new(Degrees(0.0), Degrees(180.0));
244 let idl = Vector3d::from(&lat_lon_idl);
245
246 let equator_intersection = calculate_intersection(&south_pole, &north_pole, MIN_SQ_NORM);
247 assert!(equator_intersection.is_none());
248
249 let gc_intersection1 = calculate_intersection(&idl, &north_pole, MIN_SQ_NORM).unwrap();
250 let gc_intersection2 = calculate_intersection(&idl, &south_pole, MIN_SQ_NORM).unwrap();
251
252 assert_eq!(gc_intersection1, -gc_intersection2);
253 }
254
255 #[test]
256 fn test_calculate_intersection_distances() {
257 let start1 = LatLong::new(Degrees(-1.0), Degrees(-1.0));
258 let a1 = Vector3d::from(&start1);
259 let azimuth1 = Angle::from(Degrees(45.0));
260 let pole1 = vector::calculate_pole(
261 Angle::from(start1.lat()),
262 Angle::from(start1.lon()),
263 azimuth1,
264 );
265
266 let start2 = LatLong::new(Degrees(1.0), Degrees(-1.0));
267 let a2 = Vector3d::from(&start2);
268 let azimuth2 = Angle::from(Degrees(135.0));
269 let pole2 = vector::calculate_pole(
270 Angle::from(start2.lat()),
271 Angle::from(start2.lon()),
272 azimuth2,
273 );
274
275 let c = calculate_intersection(&pole1, &pole2, MIN_SQ_NORM).unwrap();
276 let (c1, c2) = calculate_intersection_distances(&a1, &pole1, &a2, &pole2, &c);
277 assert!(is_within_tolerance(-3.1169124762478333, c1.0, f64::EPSILON));
278 assert!(is_within_tolerance(-3.1169124762478333, c2.0, f64::EPSILON));
279
280 // Calculate the centre of the arc start points
281 let centre_point = vector::normalise(&(a1 + a2), MIN_SQ_NORM).unwrap();
282 assert!(sq_distance(&c, ¢re_point) > 2.0);
283
284 // opposite intersection point
285 let d = -c;
286 assert!(sq_distance(&d, ¢re_point) <= 2.0);
287
288 let (d1, d2) = calculate_intersection_distances(&a1, &pole1, &a2, &pole2, &d);
289 assert!(is_within_tolerance(
290 0.024680177341956263,
291 d1.0,
292 f64::EPSILON
293 ));
294 assert!(is_within_tolerance(
295 0.024680177341956263,
296 d2.0,
297 f64::EPSILON
298 ));
299
300 // Same start points and intersection point
301 let (e1, e2) = calculate_intersection_distances(&a1, &pole1, &a1, &pole2, &a1);
302 assert_eq!(0.0, e1.0);
303 assert_eq!(0.0, e2.0);
304 }
305
306 #[test]
307 fn test_is_within() {
308 assert!(!is_within(-2.0 * f64::EPSILON, 2.0));
309 assert!(is_within(-f64::EPSILON, 2.0));
310 assert!(is_within(2.0 * (1.0 + f64::EPSILON), 2.0));
311 assert!(!is_within(2.0 * (1.0 + 3.0 * f64::EPSILON), 2.0));
312 }
313
314 #[test]
315 fn test_calculate_coincident_arc_distances() {
316 let zero = Radians(0.0);
317 let length1 = Radians(0.25);
318 let length2 = Radians(0.75);
319
320 let result0 = calculate_coincident_arc_distances(length2, true, length2, length1);
321 assert_eq!(length2, result0.0);
322 assert_eq!(zero, result0.1);
323
324 let result1 = calculate_coincident_arc_distances(length2, true, length1, length2);
325 assert_eq!(zero, result1.0);
326 assert_eq!(length2, result1.1);
327
328 let result2 = calculate_coincident_arc_distances(Radians(1.0), true, length1, length2);
329 assert_eq!(length1, result2.0);
330 assert_eq!(length2, result2.1);
331
332 let result3 = calculate_coincident_arc_distances(Radians(1.5), true, length1, length2);
333 assert_eq!(length2, result3.0);
334 assert_eq!(length2, result3.1);
335
336 let result4 = calculate_coincident_arc_distances(Radians(-1.5), true, length1, length2);
337 assert_eq!(Radians(-1.5), result4.0);
338 assert_eq!(zero, result4.1);
339
340 let result5 = calculate_coincident_arc_distances(Radians(-1.0), false, length1, length2);
341 assert_eq!(zero, result5.0);
342 assert_eq!(Radians(1.0), result5.1);
343
344 let result6 = calculate_coincident_arc_distances(Radians(1.0), false, length1, length2);
345 assert_eq!(Radians(1.0), result6.0);
346 assert_eq!(zero, result6.1);
347
348 let result7 = calculate_coincident_arc_distances(-length2, false, length1, length2);
349 assert_eq!(zero, result7.0);
350 assert_eq!(length2, result7.1);
351
352 let result8 = calculate_coincident_arc_distances(length1, false, length1, length2);
353 assert_eq!(length1, result8.0);
354 assert_eq!(zero, result8.1);
355 }
356}