jord/spherical/
cap.rs

1use std::f64::consts::PI;
2
3use crate::{Angle, LatLong, Mat33, NVector, Vec3};
4
5use super::{ChordLength, Sphere};
6
7/// A [spherical cap](https://en.wikipedia.org/wiki/Spherical_cap): a portion of a sphere cut off by a plane.
8/// This struct and implementation is very much based on [S2Cap](https://github.com/google/s2geometry/blob/master/src/s2/s2cap.h).
9#[derive(PartialEq, Clone, Copy, Debug, Default)]
10#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))] // codecov:ignore:this
11pub struct Cap {
12    centre: NVector,
13    radius: ChordLength,
14}
15
16impl Cap {
17    /// Empty spherical cap: contains no position.
18    pub const EMPTY: Cap = Self {
19        centre: NVector::new(Vec3::UNIT_Z),
20        radius: ChordLength::NEGATIVE,
21    };
22
23    /// Full spherical cap: contains all positions.
24    pub const FULL: Cap = Self {
25        centre: NVector::new(Vec3::UNIT_Z),
26        radius: ChordLength::MAX,
27    };
28
29    /// Constructs a new cap from the given centre and given radius expressed as the angle between the
30    /// centre and all positions on the boundary of the cap.
31    pub fn from_centre_and_radius(centre: NVector, radius: Angle) -> Self {
32        Self {
33            centre,
34            radius: ChordLength::from_angle(radius),
35        }
36    }
37
38    /// Constructs a new cap from the given centre and a given position on the boundary of the cap.
39    pub fn from_centre_and_boundary_position(centre: NVector, boundary_position: NVector) -> Self {
40        Self {
41            centre,
42            radius: ChordLength::new(centre, boundary_position),
43        }
44    }
45
46    /// Constructs a new cap whose boundary passes by the 3 given positions: the returned cap is the circumcircle of the
47    /// triangle defined by the 3 given positions.
48    pub fn from_triangle(a: NVector, b: NVector, c: NVector) -> Self {
49        // see STRIPACK: http://orion.math.iastate.edu/burkardt/f_src/stripack/stripack.f90
50        // 3 positions must be in anti-clockwise order
51        let clockwise = Sphere::side(a, b, c) < 0;
52        let v1 = a.as_vec3();
53        let v2 = if clockwise { c.as_vec3() } else { b.as_vec3() };
54        let v3 = if clockwise { b.as_vec3() } else { c.as_vec3() };
55        let e1 = v2 - v1;
56        let e2 = v3 - v1;
57        let centre = NVector::new(e1.orthogonal_to(e2));
58        // all chord length should be equal, still take maximum to account for floating point errors.
59        let radius: ChordLength = ChordLength::new(a, centre)
60            .max(ChordLength::new(b, centre).max(ChordLength::new(c, centre)));
61        Self { centre, radius }
62    }
63
64    /// Determines whether this cap is [full](crate::spherical::Cap::FULL).
65    pub fn is_full(&self) -> bool {
66        self.radius == ChordLength::MAX
67    }
68
69    /// Determines whether this cap is [empty](crate::spherical::Cap::EMPTY).
70    pub fn is_empty(&self) -> bool {
71        self.radius == ChordLength::NEGATIVE
72    }
73
74    /// Returns the complement of this cap. Both caps have the same boundary but
75    /// disjoint interiors (the union of both caps is [full](crate::spherical::Cap::FULL)).
76    pub fn complement(&self) -> Self {
77        if self.is_empty() {
78            Self::FULL
79        } else if self.is_full() {
80            Self::EMPTY
81        } else {
82            Self {
83                centre: self.centre.antipode(),
84                radius: ChordLength::from_squared_length(
85                    ChordLength::MAX.length2() - self.radius.length2(),
86                ),
87            }
88        }
89    }
90
91    /// Determines whether this cap contains the given position (including the boundary).
92    ///
93    /// # Examples
94    ///
95    /// ```
96    /// use jord::NVector;
97    /// use jord::spherical::Cap;
98    ///
99    /// let cap = Cap::from_centre_and_boundary_position(
100    ///     NVector::from_lat_long_degrees(90.0, 0.0),
101    ///     NVector::from_lat_long_degrees(0.0, 0.0)
102    /// );
103    ///
104    /// assert!(cap.contains_position(NVector::from_lat_long_degrees(0.0, 0.0)));
105    /// assert!(cap.contains_position(NVector::from_lat_long_degrees(45.0, 45.0)));
106    /// ```
107    pub fn contains_position(&self, p: NVector) -> bool {
108        ChordLength::new(self.centre, p) <= self.radius
109    }
110
111    /// Determines whether the interior of this cap contains the given position.
112    ///
113    /// # Examples
114    ///
115    /// ```
116    /// use jord::NVector;
117    /// use jord::spherical::Cap;
118    ///
119    /// let cap = Cap::from_centre_and_boundary_position(
120    ///     NVector::from_lat_long_degrees(90.0, 0.0),
121    ///     NVector::from_lat_long_degrees(0.0, 0.0)
122    /// );
123    ///
124    /// assert!(!cap.interior_contains_position(NVector::from_lat_long_degrees(0.0, 0.0)));
125    /// assert!(cap.interior_contains_position(NVector::from_lat_long_degrees(45.0, 45.0)));
126    /// ```
127    pub fn interior_contains_position(&self, p: NVector) -> bool {
128        ChordLength::new(self.centre, p) < self.radius
129    }
130
131    /// Determines whether this cap contains the given cap. The full cap contains all caps
132    /// and the empty cap is contained by all caps.
133    ///
134    /// # Examples
135    ///
136    /// ```
137    /// use jord::NVector;
138    /// use jord::spherical::Cap;
139    ///
140    /// let cap1 = Cap::from_centre_and_boundary_position(
141    ///     NVector::from_lat_long_degrees(90.0, 0.0),
142    ///     NVector::from_lat_long_degrees(0.0, 0.0)
143    /// );
144    ///
145    /// let cap2 = Cap::from_centre_and_boundary_position(
146    ///     NVector::from_lat_long_degrees(90.0, 0.0),
147    ///     NVector::from_lat_long_degrees(45.0, 0.0)
148    /// );
149    ///
150    /// assert!(cap1.contains_cap(cap2));
151    /// ```
152    pub fn contains_cap(&self, other: Self) -> bool {
153        if self.is_full() || other.is_empty() {
154            true
155        } else {
156            self.radius.length2()
157                >= ChordLength::new(self.centre, other.centre).length2() + other.radius.length2()
158        }
159    }
160
161    /// Returns the smallest cap which encloses this cap and the other given cap.
162    pub fn union(&self, other: Self) -> Self {
163        if self.radius < other.radius {
164            return other.union(*self);
165        }
166        if self.is_full() || other.is_empty() {
167            return *self;
168        }
169
170        let self_radius = self.radius();
171        let other_radius = other.radius();
172        let distance = Sphere::angle(self.centre, other.centre);
173        if self_radius >= distance + other_radius {
174            return *self;
175        }
176        let union_radius = 0.5 * (distance + self_radius + other_radius);
177        let ang = 0.5 * (distance - self_radius + other_radius);
178        let centre = Sphere::position_on_great_circle(self.centre, other.centre, ang);
179        Self {
180            centre,
181            radius: ChordLength::from_angle(union_radius),
182        }
183    }
184
185    /// Returns the centre of this cap.
186    pub fn centre(&self) -> NVector {
187        self.centre
188    }
189
190    /// Returns the radius of this cap: central angle between the centre of this cap and
191    /// any position on the boundary (negative for [empty](crate::spherical::Cap::EMPTY) caps).
192    /// The returned value may not exactly equal the value passed
193    /// to [from_centre_and_boundary_position](crate::spherical::Cap::from_centre_and_boundary_position).
194    ///
195    /// # Examples
196    ///
197    /// ```
198    /// use std::f64::consts::PI;
199    ///
200    /// use jord::{Angle, NVector};
201    /// use jord::spherical::Cap;
202    ///
203    /// let cap = Cap::from_centre_and_boundary_position(
204    ///      NVector::from_lat_long_degrees(90.0, 0.0),
205    ///      NVector::from_lat_long_degrees(45.0, 45.0)
206    /// );
207    ///
208    /// assert_eq!(Angle::from_radians(PI / 4.0), cap.radius().round_d7());
209    /// ```
210    pub fn radius(&self) -> Angle {
211        self.radius.to_angle()
212    }
213
214    /// Returns the list of vertices defining the boundary of this cap. If this cap is [empty](crate::spherical::Cap::EMPTY)
215    /// or [full](crate::spherical::Cap::FULL) the returned vector is empty, otherwise it contains `max(3, nb_vertices)` vertices.
216    ///
217    /// ```
218    /// use jord::{Angle, NVector};
219    /// use jord::spherical::{Cap, Sphere};
220    ///
221    /// let centre = NVector::from_lat_long_degrees(55.6050, 13.0038);
222    /// let radius = Angle::from_degrees(1.0);
223    /// let cap = Cap::from_centre_and_radius(centre, radius);
224    ///
225    /// let vs = cap.boundary(10);
226    /// for v in vs {
227    ///     assert_eq!(radius, Sphere::angle(centre, v).round_d7());
228    /// }
229    /// ```
230    pub fn boundary(&self, nb_vertices: usize) -> Vec<NVector> {
231        if self.is_empty() || self.is_full() {
232            return Vec::new();
233        }
234
235        let radius = self.radius().as_radians();
236        let rm = radius.sin();
237        let z = (1.0 - rm * rm).sqrt();
238
239        let ll = LatLong::from_nvector(self.centre);
240        let lat = ll.latitude().as_radians();
241        let lon = ll.longitude().as_radians();
242
243        let rya = PI / 2.0 - lat;
244        let cy = rya.cos();
245        let sy = rya.sin();
246        let ry = Mat33::new(
247            Vec3::new(cy, 0.0, sy),
248            Vec3::new(0.0, 1.0, 0.0),
249            Vec3::new(-sy, 0.0, cy),
250        );
251
252        let rza = lon;
253        let cz = rza.cos();
254        let sz = rza.sin();
255        let rz = Mat33::new(
256            Vec3::new(cz, -sz, 0.0),
257            Vec3::new(sz, cz, 0.0),
258            Vec3::new(0.0, 0.0, 1.0),
259        );
260
261        let n = nb_vertices.max(3);
262
263        let mut angles = Vec::with_capacity(n);
264        let mut r = 0.0;
265        let inc = (2.0 * PI) / (n as f64);
266        for _i in 0..n {
267            angles.push(r);
268            r += inc;
269        }
270
271        let mut res = Vec::with_capacity(n);
272        for a in angles {
273            // arc at north pole.
274            let a_np = Vec3::new(-rm * a.cos(), rm * a.sin(), z);
275            // rotate each position to arc centre.
276            let a_cen = (a_np * ry) * rz;
277
278            let p = NVector::new(a_cen.unit());
279            res.push(p);
280        }
281        res
282    }
283}
284
285#[cfg(test)]
286mod tests {
287    use crate::{positions::assert_nv_eq_d7, spherical::Cap, Angle, LatLong, NVector};
288    use std::f64::consts::PI;
289
290    #[test]
291    fn full() {
292        assert!(Cap::FULL.contains_position(NVector::from_lat_long_degrees(90.0, 0.0)));
293        assert!(Cap::FULL.contains_position(NVector::from_lat_long_degrees(-90.0, 0.0)));
294        assert_eq!(Angle::from_radians(PI), Cap::FULL.radius());
295        assert_eq!(Cap::EMPTY, Cap::FULL.complement());
296    }
297
298    #[test]
299    fn empty() {
300        assert!(!Cap::EMPTY.contains_position(NVector::from_lat_long_degrees(90.0, 0.0)));
301        assert!(!Cap::EMPTY.contains_position(NVector::from_lat_long_degrees(-90.0, 0.0)));
302        assert_eq!(Angle::from_radians(-1.0), Cap::EMPTY.radius());
303        assert_eq!(Cap::FULL, Cap::EMPTY.complement());
304    }
305
306    #[test]
307    fn from_triangle() {
308        let a = NVector::from_lat_long_degrees(0.0, 0.0);
309        let b = NVector::from_lat_long_degrees(20.0, 0.0);
310        let c = NVector::from_lat_long_degrees(10.0, 10.0);
311        let cap = Cap::from_triangle(a, b, c);
312        assert!(cap.contains_position(a));
313        assert!(cap.contains_position(b));
314        assert!(cap.contains_position(c));
315
316        let o = Cap::from_triangle(c, b, a);
317        assert_nv_eq_d7(o.centre, cap.centre);
318        assert!((o.radius.length2() - cap.radius.length2()).abs() < 1e-16);
319    }
320
321    #[test]
322    fn complement() {
323        let np = NVector::from_lat_long_degrees(90.0, 0.0);
324        let sp = NVector::from_lat_long_degrees(-90.0, 0.0);
325        let northern = Cap::from_centre_and_radius(np, Angle::QUARTER_CIRCLE);
326        let southern = Cap::from_centre_and_radius(sp, Angle::QUARTER_CIRCLE);
327
328        let northern_complement = northern.complement();
329        assert_eq!(southern.centre, northern_complement.centre);
330        assert!((southern.radius.length2() - northern_complement.radius.length2()).abs() < 1e15);
331
332        let southern_complement = southern.complement();
333        assert_eq!(northern.centre, southern_complement.centre);
334        assert!((northern.radius.length2() - southern_complement.radius.length2()).abs() < 1e15);
335    }
336
337    #[test]
338    fn contains_position() {
339        let cap = Cap::from_centre_and_boundary_position(
340            NVector::from_lat_long_degrees(90.0, 0.0),
341            NVector::from_lat_long_degrees(0.0, 0.0),
342        );
343        assert!(cap.contains_position(NVector::from_lat_long_degrees(0.0, 0.0)));
344        assert!(cap.contains_position(NVector::from_lat_long_degrees(45.0, 45.0)));
345    }
346
347    #[test]
348    fn interior_contains_position() {
349        let cap = Cap::from_centre_and_boundary_position(
350            NVector::from_lat_long_degrees(90.0, 0.0),
351            NVector::from_lat_long_degrees(0.0, 0.0),
352        );
353        assert!(!cap.interior_contains_position(NVector::from_lat_long_degrees(0.0, 0.0)));
354        assert!(cap.interior_contains_position(NVector::from_lat_long_degrees(45.0, 45.0)));
355    }
356
357    #[test]
358    fn contains_cap() {
359        let c = Cap::from_centre_and_radius(
360            NVector::from_lat_long_degrees(30.0, 30.0),
361            Angle::from_degrees(10.0),
362        );
363        assert!(Cap::FULL.contains_cap(c));
364        assert!(c.contains_cap(Cap::EMPTY));
365
366        let o = Cap::from_centre_and_radius(
367            NVector::from_lat_long_degrees(30.0, 30.0),
368            Angle::from_degrees(20.0),
369        );
370        assert!(!c.contains_cap(o));
371        assert!(o.contains_cap(c));
372    }
373
374    #[test]
375    fn radius() {
376        assert_eq!(
377            Angle::QUARTER_CIRCLE,
378            Cap::from_centre_and_boundary_position(
379                NVector::from_lat_long_degrees(90.0, 0.0),
380                NVector::from_lat_long_degrees(0.0, 0.0)
381            )
382            .radius()
383            .round_d7()
384        );
385        assert_eq!(
386            Angle::from_radians(PI / 4.0),
387            Cap::from_centre_and_boundary_position(
388                NVector::from_lat_long_degrees(90.0, 0.0),
389                NVector::from_lat_long_degrees(45.0, 45.0)
390            )
391            .radius()
392            .round_d7()
393        );
394    }
395
396    #[test]
397    fn union() {
398        assert!(Cap::FULL.union(Cap::EMPTY).is_full());
399        assert!(Cap::EMPTY.union(Cap::FULL).is_full());
400
401        // a and b have same centre, but radius of a  < radius of b.
402        let a = Cap::from_centre_and_radius(
403            NVector::from_lat_long_degrees(50.0, 10.0),
404            Angle::from_degrees(0.2),
405        );
406        let b = Cap::from_centre_and_radius(
407            NVector::from_lat_long_degrees(50.0, 10.0),
408            Angle::from_degrees(0.3),
409        );
410
411        assert!(b.contains_cap(a));
412        assert_eq!(b, a.union(b));
413
414        assert_eq!(Cap::FULL, a.union(Cap::FULL));
415        assert_eq!(a, a.union(Cap::EMPTY));
416
417        // a and c have different centers, one entirely encompasses the other.
418        let c = Cap::from_centre_and_radius(
419            NVector::from_lat_long_degrees(51.0, 11.0),
420            Angle::from_degrees(1.5),
421        );
422        assert!(c.contains_cap(a));
423        assert_eq!(a.union(c).centre(), c.centre());
424        assert_eq!(a.union(c).radius(), c.radius());
425
426        // e is partially overlapping a.
427        let e = Cap::from_centre_and_radius(
428            NVector::from_lat_long_degrees(50.3, 10.3),
429            Angle::from_degrees(0.2),
430        );
431        assert!(!e.contains_cap(a));
432        let u = a.union(e);
433        let c = LatLong::from_nvector(u.centre());
434        assert_eq!(Angle::from_degrees(50.1501), c.latitude().round_d5());
435        assert_eq!(Angle::from_degrees(10.14953), c.longitude().round_d5());
436        assert_eq!(Angle::from_degrees(0.37815), u.radius().round_d5());
437    }
438
439    #[test]
440    fn boundary() {
441        assert!(Cap::EMPTY.boundary(1).is_empty());
442        assert!(Cap::FULL.boundary(1).is_empty());
443        let northern = Cap::from_centre_and_radius(
444            NVector::from_lat_long_degrees(90.0, 0.0),
445            Angle::QUARTER_CIRCLE,
446        );
447        assert_eq!(
448            vec![
449                LatLong::from_degrees(0.0, 180.0),
450                LatLong::from_degrees(0.0, 90.0),
451                LatLong::from_degrees(0.0, 0.0),
452                LatLong::from_degrees(0.0, -90.0)
453            ],
454            northern
455                .boundary(4)
456                .iter()
457                .map(|v| LatLong::from_nvector(*v).round_d7())
458                .collect::<Vec<_>>()
459        );
460        assert_eq!(
461            vec![
462                LatLong::from_degrees(0.0, 180.0),
463                LatLong::from_degrees(0.0, 60.0),
464                LatLong::from_degrees(0.0, -60.0)
465            ],
466            northern
467                .boundary(2)
468                .iter()
469                .map(|v| LatLong::from_nvector(*v).round_d7())
470                .collect::<Vec<_>>()
471        );
472    }
473}
jord/spherical/cap.rs

jord/spherical/
cap.rs
