Skip to main content

sphereql_core/
regions.rs

1use std::f64::consts::{PI, TAU};
2
3use crate::distance::angular_distance;
4use crate::error::SphereQlError;
5use crate::types::SphericalPoint;
6
7/// Spatial containment test for a [`SphericalPoint`].
8pub trait Contains {
9    /// Returns `true` if `point` lies inside (or on the boundary of) this region.
10    fn contains(&self, point: &SphericalPoint) -> bool;
11}
12
13/// A directional cone anchored at the sphere origin.
14///
15/// `apex` is retained for serialization and downstream tooling but is
16/// not consulted by [`Cone::contains`]; containment is decided purely by
17/// the angular distance from the test point to `axis`.
18#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
19pub struct Cone {
20    pub apex: SphericalPoint,
21    pub axis: SphericalPoint,
22    pub half_angle: f64,
23}
24
25impl Cone {
26    /// Creates a new `Cone` with validation.
27    ///
28    /// Returns an error if `half_angle` is not in (0, π].
29    pub fn new(
30        apex: SphericalPoint,
31        axis: SphericalPoint,
32        half_angle: f64,
33    ) -> Result<Self, SphereQlError> {
34        if half_angle <= 0.0 || half_angle > PI {
35            return Err(SphereQlError::InvalidConeAngle(half_angle));
36        }
37        Ok(Self {
38            apex,
39            axis,
40            half_angle,
41        })
42    }
43}
44
45impl Contains for Cone {
46    fn contains(&self, point: &SphericalPoint) -> bool {
47        let point_unit = SphericalPoint::new_unchecked(1.0, point.theta, point.phi);
48        let axis_unit = SphericalPoint::new_unchecked(1.0, self.axis.theta, self.axis.phi);
49        angular_distance(&point_unit, &axis_unit) <= self.half_angle
50    }
51}
52
53/// A spherical cap: all points within `half_angle` radians of `center` on S².
54///
55/// Unlike `Cone`, a cap has no apex — containment depends only on angular distance
56/// from `center`, regardless of radial distance.
57#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
58pub struct Cap {
59    pub center: SphericalPoint,
60    pub half_angle: f64,
61}
62
63impl Cap {
64    /// Creates a new `Cap` with validation.
65    ///
66    /// Returns an error if `half_angle` is not in (0, π].
67    pub fn new(center: SphericalPoint, half_angle: f64) -> Result<Self, SphereQlError> {
68        if half_angle <= 0.0 || half_angle > PI {
69            return Err(SphereQlError::InvalidCapAngle(half_angle));
70        }
71        Ok(Self { center, half_angle })
72    }
73}
74
75impl Contains for Cap {
76    fn contains(&self, point: &SphericalPoint) -> bool {
77        let point_unit = SphericalPoint::new_unchecked(1.0, point.theta, point.phi);
78        let center_unit = SphericalPoint::new_unchecked(1.0, self.center.theta, self.center.phi);
79        angular_distance(&point_unit, &center_unit) <= self.half_angle
80    }
81}
82
83/// A hollow spherical shell bounded by two radii.
84///
85/// Contains all points with `inner <= r <= outer`.
86#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
87pub struct Shell {
88    pub inner: f64,
89    pub outer: f64,
90}
91
92impl Shell {
93    /// Creates a new `Shell` with validation.
94    ///
95    /// Returns an error if `inner < 0` or `inner >= outer`.
96    pub fn new(inner: f64, outer: f64) -> Result<Self, SphereQlError> {
97        if inner < 0.0 || inner >= outer {
98            return Err(SphereQlError::InvalidShellBounds { inner, outer });
99        }
100        Ok(Self { inner, outer })
101    }
102}
103
104impl Contains for Shell {
105    fn contains(&self, point: &SphericalPoint) -> bool {
106        point.r >= self.inner && point.r <= self.outer
107    }
108}
109
110/// A latitude band: all points with polar angle `phi` in [`phi_min`, `phi_max`].
111///
112/// `phi = 0` is the north pole; `phi = π` is the south pole.
113#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
114pub struct Band {
115    pub phi_min: f64,
116    pub phi_max: f64,
117}
118
119impl Band {
120    /// Creates a new `Band` with validation.
121    ///
122    /// Returns an error if `phi_min < 0`, `phi_min >= phi_max`, or `phi_max > π`.
123    pub fn new(phi_min: f64, phi_max: f64) -> Result<Self, SphereQlError> {
124        if phi_min < 0.0 || phi_min >= phi_max || phi_max > PI {
125            return Err(SphereQlError::InvalidBandBounds { phi_min, phi_max });
126        }
127        Ok(Self { phi_min, phi_max })
128    }
129}
130
131impl Contains for Band {
132    fn contains(&self, point: &SphericalPoint) -> bool {
133        point.phi >= self.phi_min && point.phi <= self.phi_max
134    }
135}
136
137/// A longitude wedge: all points with azimuthal angle `theta` between `theta_min` and
138/// `theta_max`.
139///
140/// Supports wrap-around: if `theta_min > theta_max`, the wedge spans across the 0/2π
141/// boundary (e.g. 350° to 10°).
142#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
143pub struct Wedge {
144    pub theta_min: f64,
145    pub theta_max: f64,
146}
147
148impl Wedge {
149    /// Creates a new `Wedge` with validation.
150    ///
151    /// Returns an error if either bound is outside [0, 2π) or if `theta_min == theta_max`.
152    pub fn new(theta_min: f64, theta_max: f64) -> Result<Self, SphereQlError> {
153        if !(0.0..TAU).contains(&theta_min)
154            || !(0.0..TAU).contains(&theta_max)
155            || theta_min == theta_max
156        {
157            return Err(SphereQlError::InvalidWedgeBounds {
158                theta_min,
159                theta_max,
160            });
161        }
162        Ok(Self {
163            theta_min,
164            theta_max,
165        })
166    }
167
168    fn wraps(&self) -> bool {
169        self.theta_min > self.theta_max
170    }
171}
172
173impl Contains for Wedge {
174    fn contains(&self, point: &SphericalPoint) -> bool {
175        if self.wraps() {
176            point.theta >= self.theta_min || point.theta <= self.theta_max
177        } else {
178            point.theta >= self.theta_min && point.theta <= self.theta_max
179        }
180    }
181}
182
183/// A composable spatial region on S².
184///
185/// Regions can be primitive shapes (`Cone`, `Cap`, `Shell`, `Band`, `Wedge`) or
186/// boolean combinations (`Intersection`, `Union`). All implement [`Contains`].
187#[derive(Debug, Clone, PartialEq, serde::Serialize, serde::Deserialize)]
188#[non_exhaustive]
189pub enum Region {
190    Cone(Cone),
191    Cap(Cap),
192    Shell(Shell),
193    Band(Band),
194    Wedge(Wedge),
195    Intersection(Vec<Region>),
196    Union(Vec<Region>),
197}
198
199impl Region {
200    /// Creates a region that contains points inside **all** of the given regions.
201    pub fn intersection(regions: Vec<Region>) -> Self {
202        Region::Intersection(regions)
203    }
204
205    /// Creates a region that contains points inside **any** of the given regions.
206    pub fn union(regions: Vec<Region>) -> Self {
207        Region::Union(regions)
208    }
209}
210
211impl Contains for Region {
212    fn contains(&self, point: &SphericalPoint) -> bool {
213        match self {
214            Region::Cone(c) => c.contains(point),
215            Region::Cap(c) => c.contains(point),
216            Region::Shell(s) => s.contains(point),
217            Region::Band(b) => b.contains(point),
218            Region::Wedge(w) => w.contains(point),
219            Region::Intersection(regions) => regions.iter().all(|r| r.contains(point)),
220            Region::Union(regions) => regions.iter().any(|r| r.contains(point)),
221        }
222    }
223}
224
225#[cfg(test)]
226mod tests {
227    use super::*;
228    use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, PI};
229
230    fn point(r: f64, theta: f64, phi: f64) -> SphericalPoint {
231        SphericalPoint::new_unchecked(r, theta, phi)
232    }
233
234    // --- Cone tests ---
235
236    #[test]
237    fn cone_contains_point_inside() {
238        let cone = Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, FRAC_PI_4), FRAC_PI_2).unwrap();
239        let p = point(2.0, 0.0, FRAC_PI_4);
240        assert!(cone.contains(&p));
241    }
242
243    #[test]
244    fn cone_excludes_point_outside() {
245        let cone = Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.1), 0.05).unwrap();
246        let p = point(1.0, PI, FRAC_PI_2);
247        assert!(!cone.contains(&p));
248    }
249
250    #[test]
251    fn cone_contains_point_on_boundary() {
252        let half = FRAC_PI_4;
253        let cone = Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.0), half).unwrap();
254        // axis is at phi=0 (north pole), point at phi=half_angle should be on boundary
255        let p = point(1.0, 0.0, half);
256        assert!(cone.contains(&p));
257    }
258
259    #[test]
260    fn cone_various_half_angles() {
261        // Narrow cone
262        let narrow = Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, FRAC_PI_2), 0.01).unwrap();
263        let near = point(1.0, 0.0, FRAC_PI_2);
264        let far = point(1.0, 0.0, FRAC_PI_2 + 0.1);
265        assert!(narrow.contains(&near));
266        assert!(!narrow.contains(&far));
267
268        // Full hemisphere
269        let wide = Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, FRAC_PI_2), FRAC_PI_2).unwrap();
270        assert!(wide.contains(&point(1.0, 0.5, FRAC_PI_2 + 0.3)));
271    }
272
273    #[test]
274    fn cone_invalid_half_angle() {
275        assert!(Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.0), 0.0).is_err());
276        assert!(Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.0), -0.1).is_err());
277        assert!(Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.0), PI + 0.1).is_err());
278        // PI is valid
279        assert!(Cone::new(point(0.0, 0.0, 0.0), point(1.0, 0.0, 0.0), PI).is_ok());
280    }
281
282    // --- Cap tests ---
283
284    #[test]
285    fn cap_contains_point_inside() {
286        let cap = Cap::new(point(1.0, 0.0, FRAC_PI_2), FRAC_PI_4).unwrap();
287        let p = point(5.0, 0.0, FRAC_PI_2);
288        assert!(cap.contains(&p));
289    }
290
291    #[test]
292    fn cap_excludes_point_outside() {
293        let cap = Cap::new(point(1.0, 0.0, 0.1), 0.05).unwrap();
294        let p = point(1.0, PI, PI - 0.1);
295        assert!(!cap.contains(&p));
296    }
297
298    #[test]
299    fn cap_ignores_radius() {
300        let cap = Cap::new(point(1.0, 0.0, FRAC_PI_2), FRAC_PI_4).unwrap();
301        let near = point(0.1, 0.0, FRAC_PI_2);
302        let far = point(1000.0, 0.0, FRAC_PI_2);
303        assert!(cap.contains(&near));
304        assert!(cap.contains(&far));
305    }
306
307    #[test]
308    fn cap_invalid_half_angle() {
309        assert!(Cap::new(point(1.0, 0.0, 0.0), 0.0).is_err());
310        assert!(Cap::new(point(1.0, 0.0, 0.0), -1.0).is_err());
311    }
312
313    #[test]
314    fn cap_error_is_cap_specific() {
315        let err = Cap::new(point(1.0, 0.0, 0.0), 0.0).unwrap_err();
316        assert!(
317            matches!(err, SphereQlError::InvalidCapAngle(_)),
318            "expected InvalidCapAngle, got {err:?}"
319        );
320    }
321
322    // --- Shell tests ---
323
324    #[test]
325    fn shell_contains_point_inside() {
326        let shell = Shell::new(1.0, 5.0).unwrap();
327        assert!(shell.contains(&point(3.0, 0.0, 0.0)));
328    }
329
330    #[test]
331    fn shell_excludes_point_outside() {
332        let shell = Shell::new(1.0, 5.0).unwrap();
333        assert!(!shell.contains(&point(0.5, 0.0, 0.0)));
334        assert!(!shell.contains(&point(6.0, 0.0, 0.0)));
335    }
336
337    #[test]
338    fn shell_boundary_inclusive() {
339        let shell = Shell::new(1.0, 5.0).unwrap();
340        assert!(shell.contains(&point(1.0, 0.0, 0.0)));
341        assert!(shell.contains(&point(5.0, 0.0, 0.0)));
342    }
343
344    #[test]
345    fn shell_invalid_bounds() {
346        assert!(Shell::new(5.0, 1.0).is_err());
347        assert!(Shell::new(3.0, 3.0).is_err());
348        assert!(Shell::new(-1.0, 5.0).is_err());
349    }
350
351    // --- Band tests ---
352
353    #[test]
354    fn band_contains_point_inside() {
355        let band = Band::new(FRAC_PI_4, 3.0 * FRAC_PI_4).unwrap();
356        assert!(band.contains(&point(1.0, 0.0, FRAC_PI_2)));
357    }
358
359    #[test]
360    fn band_excludes_point_outside() {
361        let band = Band::new(FRAC_PI_4, FRAC_PI_2).unwrap();
362        assert!(!band.contains(&point(1.0, 0.0, 0.1)));
363        assert!(!band.contains(&point(1.0, 0.0, PI - 0.1)));
364    }
365
366    #[test]
367    fn band_boundary_inclusive() {
368        let band = Band::new(FRAC_PI_4, FRAC_PI_2).unwrap();
369        assert!(band.contains(&point(1.0, 0.0, FRAC_PI_4)));
370        assert!(band.contains(&point(1.0, 0.0, FRAC_PI_2)));
371    }
372
373    #[test]
374    fn band_poles() {
375        // Band covering north pole area
376        let north = Band::new(0.0 + 0.001, FRAC_PI_4).unwrap();
377        assert!(north.contains(&point(1.0, 0.0, 0.01)));
378        assert!(!north.contains(&point(1.0, 0.0, FRAC_PI_2)));
379
380        // Band covering south pole area
381        let south = Band::new(3.0 * FRAC_PI_4, PI).unwrap();
382        assert!(south.contains(&point(1.0, 0.0, PI - 0.1)));
383        assert!(!south.contains(&point(1.0, 0.0, FRAC_PI_4)));
384    }
385
386    #[test]
387    fn band_invalid_bounds() {
388        assert!(Band::new(FRAC_PI_2, FRAC_PI_4).is_err());
389        assert!(Band::new(FRAC_PI_4, FRAC_PI_4).is_err());
390        assert!(Band::new(-0.1, FRAC_PI_2).is_err());
391        assert!(Band::new(0.0, PI + 0.1).is_err());
392    }
393
394    // --- Wedge tests ---
395
396    #[test]
397    fn wedge_contains_normal_range() {
398        let wedge = Wedge::new(0.5, 2.0).unwrap();
399        assert!(wedge.contains(&point(1.0, 1.0, FRAC_PI_2)));
400        assert!(!wedge.contains(&point(1.0, 3.0, FRAC_PI_2)));
401    }
402
403    #[test]
404    fn wedge_wraparound() {
405        // 350° to 10° in radians: ~6.1087 to ~0.1745
406        let theta_min = 350.0_f64.to_radians();
407        let theta_max = 10.0_f64.to_radians();
408        let wedge = Wedge::new(theta_min, theta_max).unwrap();
409
410        let inside_high = point(1.0, 355.0_f64.to_radians(), FRAC_PI_2);
411        let inside_low = point(1.0, 5.0_f64.to_radians(), FRAC_PI_2);
412        let outside = point(1.0, 180.0_f64.to_radians(), FRAC_PI_2);
413
414        assert!(wedge.contains(&inside_high));
415        assert!(wedge.contains(&inside_low));
416        assert!(!wedge.contains(&outside));
417    }
418
419    #[test]
420    fn wedge_boundary_inclusive() {
421        let wedge = Wedge::new(1.0, 2.0).unwrap();
422        assert!(wedge.contains(&point(1.0, 1.0, FRAC_PI_2)));
423        assert!(wedge.contains(&point(1.0, 2.0, FRAC_PI_2)));
424    }
425
426    #[test]
427    fn wedge_rejects_invalid_theta() {
428        assert!(Wedge::new(-0.1, 1.0).is_err());
429        assert!(Wedge::new(0.0, 7.0).is_err());
430    }
431
432    // --- Compound region tests ---
433
434    #[test]
435    fn intersection_shell_and_band() {
436        let shell = Region::Shell(Shell::new(1.0, 5.0).unwrap());
437        let band = Region::Band(Band::new(FRAC_PI_4, 3.0 * FRAC_PI_4).unwrap());
438        let region = Region::intersection(vec![shell, band]);
439
440        // Inside both
441        assert!(region.contains(&point(3.0, 0.0, FRAC_PI_2)));
442        // Inside shell but outside band
443        assert!(!region.contains(&point(3.0, 0.0, 0.1)));
444        // Inside band but outside shell
445        assert!(!region.contains(&point(10.0, 0.0, FRAC_PI_2)));
446    }
447
448    #[test]
449    fn union_two_caps() {
450        let cap_a = Region::Cap(Cap::new(point(1.0, 0.0, 0.1), 0.2).unwrap());
451        let cap_b = Region::Cap(Cap::new(point(1.0, 0.0, PI - 0.1), 0.2).unwrap());
452        let region = Region::union(vec![cap_a, cap_b]);
453
454        // Near north pole (cap_a)
455        assert!(region.contains(&point(1.0, 0.0, 0.05)));
456        // Near south pole (cap_b)
457        assert!(region.contains(&point(1.0, 0.0, PI - 0.05)));
458        // Equator (neither)
459        assert!(!region.contains(&point(1.0, 0.0, FRAC_PI_2)));
460    }
461
462    #[test]
463    fn empty_intersection_contains_everything() {
464        let region = Region::intersection(vec![]);
465        assert!(region.contains(&point(1.0, 0.0, FRAC_PI_2)));
466    }
467
468    #[test]
469    fn empty_union_contains_nothing() {
470        let region = Region::union(vec![]);
471        assert!(!region.contains(&point(1.0, 0.0, FRAC_PI_2)));
472    }
473
474    // --- Region enum dispatch ---
475
476    #[test]
477    fn region_dispatches_to_inner_types() {
478        let shell_region = Region::Shell(Shell::new(1.0, 5.0).unwrap());
479        assert!(shell_region.contains(&point(3.0, 0.0, 0.0)));
480        assert!(!shell_region.contains(&point(10.0, 0.0, 0.0)));
481
482        let wedge_region = Region::Wedge(Wedge::new(0.5, 2.0).unwrap());
483        assert!(wedge_region.contains(&point(1.0, 1.0, FRAC_PI_2)));
484    }
485}