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//! \
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::{Vector3d, calculate_great_circle_atd, normalise, normalise_centroid, sq_distance};
45use angle_sc::{Angle, Radians};
46
47/// Calculate an intersection point between the poles of two Great Circles.
48/// See: <http://www.movable-type.co.uk/scripts/latlong-vectors.html#intersection>
49/// * `pole_0`, `pole_1` the poles.
50/// * `min_sq_value` the minimum square of a vector length to normalize.
51///
52/// return an intersection point or None if the poles represent coincident Great Circles.
53#[must_use]
54pub fn calculate_intersection(
55 pole_0: &Vector3d,
56 pole_1: &Vector3d,
57 min_sq_value: f64,
58) -> Option<Vector3d> {
59 normalise(&pole_0.cross(pole_1), min_sq_value)
60}
61
62/// Determine whether the antipodal point is closer to the centroid of the
63/// `Arc`s.
64///
65/// * `point` a great-circle intersection point.
66/// * `centroid` the centroid of the `Arc`s mid points.
67///
68/// returns true if the antipodal intersection is closer to the `centroid`
69/// of the `Arc`s otherwise returns false.
70#[must_use]
71pub fn use_antipodal_point(point: &Vector3d, centroid: &Vector3d) -> bool {
72 sq_distance(centroid, &(-*point)) < sq_distance(centroid, point)
73}
74
75/// Return the closer intersection point to the centroid of the `Arc`s.
76///
77/// * `point` a great-circle intersection point.
78/// * `centroid` the centroid of the `Arc`s mid points.
79///
80/// returns the antipodal point if it is closer to the `centroid`,
81/// otherwise returns the point.
82#[must_use]
83pub fn closest_intersection_point(point: &Vector3d, centroid: &Vector3d) -> Vector3d {
84 if use_antipodal_point(point, centroid) {
85 -*point
86 } else {
87 *point
88 }
89}
90
91/// Determine the reference point of a pair of arcs.
92/// I.e. the closest intersection point if they intersect or the
93/// centroid normalized to lie on the unit sphere if they don't.
94///
95/// * `mid_point_0`, `mid_point_1` the mid points of the `Arc`s.
96/// * `pole_0`, `pole_1` the poles of the `Arc` great circles.
97/// * `sq_sin_max_coincident_angle` the square of the sine of the
98/// maximum angle between coincident great circles.
99///
100/// returns the closest intersection point or normalized centroid and the
101/// sine of the angle between the arcs, zero if the arcs are coincident.
102/// And the absolute relative angle at the intersection point or centroid.
103#[must_use]
104pub fn calculate_reference_point_and_angle(
105 mid_point_0: &Vector3d,
106 pole_0: &Vector3d,
107 mid_point_1: &Vector3d,
108 pole_1: &Vector3d,
109 sq_sin_max_coincident_angle: f64,
110) -> (Vector3d, Angle) {
111 let centroid = mid_point_0 + mid_point_1;
112 // calculate the intersection point between the great circles
113 let point = pole_0.cross(pole_1);
114 normalise(&point, sq_sin_max_coincident_angle).map_or_else(
115 || {
116 // the great circles are coincident
117
118 let c = normalise_centroid(¢roid, mid_point_0, pole_0);
119 (c, Angle::from_y_x(point.norm(), pole_0.dot(pole_1)))
120 },
121 |p| {
122 // the great circles intersect
123
124 let x = closest_intersection_point(&p, ¢roid);
125 (x, Angle::from_y_x(point.norm(), pole_0.dot(pole_1)))
126 },
127 )
128}
129
130/// Calculate signed great circle distances from two arc mid points to their
131/// closest intersection point or normalized centroid if the arcs are on coincident
132/// great circles.
133///
134/// * `mid_point_0`, `mid_point_1` the mid points of the arcs.
135/// * `pole_0`, `pole_1` the poles of the arc great circles.
136/// * `sq_sin_max_coincident_angle` the square of the sine of the
137/// maximum angle between coincident great circles.
138///
139/// returns the signed great circle distances of the closest intersection
140/// point or centroid from the arc mid points in `Radians`,
141/// and the relative angle between the arc great circles.
142#[must_use]
143pub fn calculate_arc_reference_distances_and_angle(
144 mid_point_0: &Vector3d,
145 pole_0: &Vector3d,
146 mid_point_1: &Vector3d,
147 pole_1: &Vector3d,
148 sq_sin_max_coincident_angle: f64,
149) -> (Radians, Radians, Angle) {
150 let (point, angle) = calculate_reference_point_and_angle(
151 mid_point_0,
152 pole_0,
153 mid_point_1,
154 pole_1,
155 sq_sin_max_coincident_angle,
156 );
157 let distance_0 = calculate_great_circle_atd(mid_point_0, pole_0, &point);
158 let distance_1 = calculate_great_circle_atd(mid_point_1, pole_1, &point);
159
160 (distance_0, distance_1, angle)
161}
162
163#[cfg(test)]
164mod tests {
165 use std::f64;
166
167 use super::*;
168 use crate::{LatLong, vector};
169 use angle_sc::{Angle, Degrees, is_within_tolerance};
170
171 #[test]
172 fn test_calculate_intersection() {
173 let lat_lon_south = LatLong::new(Degrees(-90.0), Degrees(0.0));
174 let south_pole = Vector3d::from(&lat_lon_south);
175
176 let lat_lon_north = LatLong::new(Degrees(90.0), Degrees(0.0));
177 let north_pole = Vector3d::from(&lat_lon_north);
178
179 let lat_lon_idl = LatLong::new(Degrees(0.0), Degrees(180.0));
180 let idl = Vector3d::from(&lat_lon_idl);
181
182 let equator_intersection =
183 calculate_intersection(&south_pole, &north_pole, vector::MIN_SQ_NORM);
184 assert!(equator_intersection.is_none());
185
186 let gc_intersection1 =
187 calculate_intersection(&idl, &north_pole, vector::MIN_SQ_NORM).unwrap();
188 let gc_intersection2 =
189 calculate_intersection(&idl, &south_pole, vector::MIN_SQ_NORM).unwrap();
190
191 assert_eq!(gc_intersection1, -gc_intersection2);
192 }
193
194 #[test]
195 fn test_calculate_arc_reference_distances_and_angle_coincident_great_circles() {
196 let point_1 = Vector3d::new(1.0, 0.0, 0.0);
197 let pole_1 = Vector3d::new(0.0, 0.0, 1.0);
198
199 // same mid points and great circles
200 let result = calculate_arc_reference_distances_and_angle(
201 &point_1,
202 &pole_1,
203 &point_1,
204 &pole_1,
205 vector::MIN_SQ_NORM,
206 );
207 assert_eq!(Radians(0.0), result.0);
208 assert_eq!(Radians(0.0), result.1);
209 assert_eq!(Degrees(0.0), Degrees::from(result.2));
210
211 // opposite mid points and same great circles
212 let point_m1 = -point_1;
213 let result = calculate_arc_reference_distances_and_angle(
214 &point_1,
215 &pole_1,
216 &point_m1,
217 &pole_1,
218 vector::MIN_SQ_NORM,
219 );
220 assert!(is_within_tolerance(
221 -f64::consts::FRAC_PI_2,
222 result.0.0,
223 f64::EPSILON
224 ));
225 assert!(is_within_tolerance(
226 f64::consts::FRAC_PI_2,
227 result.1.0,
228 f64::EPSILON
229 ));
230 assert_eq!(Degrees(0.0), Degrees::from(result.2));
231
232 // opposite mid points and great circles
233 let pole_m1 = -pole_1;
234 let result = calculate_arc_reference_distances_and_angle(
235 &point_1,
236 &pole_1,
237 &point_m1,
238 &pole_m1,
239 vector::MIN_SQ_NORM,
240 );
241 assert!(is_within_tolerance(
242 -f64::consts::FRAC_PI_2,
243 result.0.0,
244 f64::EPSILON
245 ));
246 assert!(is_within_tolerance(
247 -f64::consts::FRAC_PI_2,
248 result.1.0,
249 f64::EPSILON
250 ));
251 assert_eq!(Degrees(180.0), Degrees::from(result.2));
252 }
253
254 #[test]
255 fn test_calculate_arc_reference_distances_and_angle_intersecting_great_circles() {
256 let point_1 = Vector3d::new(1.0, 0.0, 0.0);
257 let pole_1 = Vector3d::new(0.0, 0.0, 1.0);
258 let pole_2 = Vector3d::new(0.0, 1.0, 0.0);
259
260 // intersection, same mid points
261 let result = calculate_arc_reference_distances_and_angle(
262 &point_1,
263 &pole_1,
264 &point_1,
265 &pole_2,
266 vector::MIN_SQ_NORM,
267 );
268 assert_eq!(Radians(0.0), result.0);
269 assert_eq!(Radians(0.0), result.1);
270 assert_eq!(Degrees(90.0), Degrees::from(result.2));
271
272 // intersection, same mid points, acute angle
273 let pole_3 = (pole_1 + pole_2).normalize();
274 let result = calculate_arc_reference_distances_and_angle(
275 &point_1,
276 &pole_1,
277 &point_1,
278 &pole_3,
279 vector::MIN_SQ_NORM,
280 );
281 assert_eq!(Radians(0.0), result.0);
282 assert_eq!(Radians(0.0), result.1);
283 assert_eq!(Degrees(45.0), Degrees::from(result.2));
284
285 // intersection, same mid points, obtuse angle
286 let pole_m3 = -pole_3;
287 let result = calculate_arc_reference_distances_and_angle(
288 &point_1,
289 &pole_1,
290 &point_1,
291 &pole_m3,
292 vector::MIN_SQ_NORM,
293 );
294 assert_eq!(Radians(0.0), result.0);
295 assert_eq!(Radians(0.0), result.1);
296 assert_eq!(Degrees(135.0), Degrees::from(result.2));
297
298 // intersection, different mid points, acute angle
299 let point_2 = vector::position(
300 &point_1,
301 &vector::direction(&point_1, &pole_3),
302 Angle::default().quarter_turn_cw(),
303 );
304 let result = calculate_arc_reference_distances_and_angle(
305 &point_1,
306 &pole_1,
307 &point_2,
308 &pole_3,
309 vector::MIN_SQ_NORM,
310 );
311 assert_eq!(Radians(0.0), result.0);
312 assert!(is_within_tolerance(
313 -f64::consts::FRAC_PI_2,
314 result.1.0,
315 f64::EPSILON
316 ));
317 assert_eq!(Degrees(45.0), Degrees::from(result.2));
318
319 // intersection, different mid points, obtuse angle
320 let result = calculate_arc_reference_distances_and_angle(
321 &point_1,
322 &pole_1,
323 &point_2,
324 &pole_m3,
325 vector::MIN_SQ_NORM,
326 );
327 assert_eq!(Radians(0.0), result.0);
328 assert!(is_within_tolerance(
329 f64::consts::FRAC_PI_2,
330 result.1.0,
331 f64::EPSILON
332 ));
333 assert_eq!(Degrees(135.0), Degrees::from(result.2));
334 }
335}