geop_core_geometry/
aabb.rs

1//! Cheap axis-aligned bounding box helpers shared by [`crate::nurb_curve`]
2//! and [`crate::nurb_surface`] — see `NurbCurve`/`NurbSurface`'s own `aabb`
3//! field doc comment for why one padded-to-3 `[Scalar; 3]` per curve/surface is
4//! enough to serve as a full bounding box under interval arithmetic.
5
6use geop_core_math::{scalars::Scalar, vector::Vector};
7
8use crate::nurb_curve::NurbCurve;
9
10/// Axis-aligned bounding box of `control_points`'s dehomogenized Cartesian
11/// coordinates (`D - 1` axes), padded to 3 with [`Scalar::ENTIRE`] so both
12/// 2-D pcurves (`D = 3`) and 3-D curves/surfaces (`D = 4`) share one result
13/// type — callers only ever read the first `D - 1` (see
14/// [`aabb_could_overlap`]'s own `c` parameter).
15///
16/// Each axis's bound is a *single* `Scalar`, not a separate `(min, max)` pair:
17/// `Scalar` is itself an interval type here, so folding every control point's
18/// coordinate on that axis into one `Scalar` via [`Scalar::union`] already
19/// produces exactly that axis's `[min, max]` extent — "3 scalar values" is
20/// the whole bounding box.
21///
22/// Every weight is definitely positive by construction
23/// (`NurbCurve`/`NurbSurface::try_new`), so dehomogenizing cannot fail. Each
24/// coordinate is divided by its weight, one rounding, rather than multiplied
25/// by a reciprocal, two: the box is the tightest the control points give.
26pub(crate) fn compute_aabb<const D: usize>(control_points: &[Vector<D>]) -> [Scalar; 3] {
27    let mut acc: [Option<Scalar>; 3] = [None; 3];
28    for p in control_points {
29        for (axis, slot) in acc.iter_mut().enumerate().take(D - 1) {
30            let coord = p[axis]
31                .div(p[D - 1])
32                .expect("weights are definitely positive by construction");
33            *slot = Some(match *slot {
34                None => coord,
35                Some(a) => a.union(coord),
36            });
37        }
38    }
39    acc.map(|a| a.unwrap_or(Scalar::ENTIRE))
40}
41
42/// The box around all of `boxes` (each a curve's or a surface's, see
43/// [`crate::nurb_curve::NurbCurve::aabb`]): per axis the union of theirs.
44/// None for no boxes.
45pub fn union_aabb<'a>(boxes: impl IntoIterator<Item = &'a [Scalar; 3]>) -> Option<[Scalar; 3]> {
46    boxes
47        .into_iter()
48        .copied()
49        .reduce(|a, b| std::array::from_fn(|k| a[k].union(b[k])))
50}
51
52/// True if two cached [`compute_aabb`] boxes could overlap, comparing only
53/// the first `c` (the objects' actual Cartesian dimension) of the 3 padded
54/// axes. `false` is a hard proof of separation, and cheap: `O(c)` scalar
55/// `could_be_equal`s against an already-cached value.
56pub(crate) fn aabb_could_overlap(a: &[Scalar; 3], b: &[Scalar; 3], c: usize) -> bool {
57    (0..c).all(|i| a[i].could_be_equal(b[i]))
58}
59
60/// True if a cached [`compute_aabb`] box could contain `point` — the same
61/// check as [`aabb_could_overlap`], with `point` standing in for a
62/// zero-extent "box" of its own (each of its `C` coordinates is itself an
63/// interval already, so no separate point-vs-box case is needed).
64pub(crate) fn aabb_could_contain<const C: usize>(aabb: &[Scalar; 3], point: &Vector<C>) -> bool {
65    (0..C).all(|i| aabb[i].could_be_equal(point[i]))
66}
67
68/// True if `curve` could meet the cached box `aabb` in its first `c` axes.
69///
70/// For a straight segment — degree one between two control points, which
71/// traces the segment between them whatever their weights — whether the
72/// segment itself could: far tighter than its own box when it runs
73/// diagonally, and a ray cast across a solid to classify a point is a long
74/// diagonal segment, whose box overlaps nearly every face's. Each axis
75/// confines the segment's parameter to where it lies within the box's
76/// extent on that axis; the segment meets the box only if those ranges
77/// and `[0, 1]` share a point. For any other curve, whether the two boxes
78/// could overlap. `false` is a hard proof of separation either way.
79pub(crate) fn curve_could_meet_aabb<const D: usize>(
80    curve: &NurbCurve<D>,
81    aabb: &[Scalar; 3],
82    c: usize,
83) -> bool {
84    if !aabb_could_overlap(&curve.aabb, aabb, c) {
85        return false;
86    }
87    let [p, q] = curve.control_points[..] else {
88        return true;
89    };
90    if curve.degree != 1 {
91        return true;
92    }
93    let point = |cp: Vector<D>, i: usize| cp[i].div(cp[D - 1]);
94    let mut range = Scalar::ZERO.union(Scalar::ONE);
95    for (i, extent) in aabb.iter().enumerate().take(c) {
96        let (Ok(a), Ok(b)) = (point(p, i), point(q, i)) else {
97            return true;
98        };
99        let d = b.sub(a);
100        if d.could_be_equal(Scalar::ZERO) {
101            // Hardly moving along this axis, the segment's parameter is not
102            // confined by it; the boxes overlap on it, checked above.
103            continue;
104        }
105        let (Ok(t0), Ok(t1)) = (extent.lower().sub(a).div(d), extent.upper().sub(a).div(d)) else {
106            continue;
107        };
108        let along = t0.union(t1);
109        if !along.could_be_equal(range) {
110            return false;
111        }
112        range = along.intersect(range);
113    }
114    true
115}
geop_core_geometry/aabb.rs

geop_core_geometry/
aabb.rs
