loam_shape/
isovolume.rs

1//! For a 1-Lipschitz field, the half-diagonal margin encloses its nonpositive set within the sampled domain.
2//! Hart, Sphere Tracing, 1996, §2.
3
4use glam::{Vec3, Vec4};
5
6use crate::Shape;
7
8#[derive(Clone, Copy, Debug, PartialEq, Eq)]
9struct CellBox<const D: usize> {
10    lo: [usize; D],
11    hi: [usize; D],
12}
13
14/// Encloses `{p : f(p) <= 0}` within the sampled domain; [`Self::clipped`] detects contact with its boundary.
15#[derive(Clone, Debug)]
16pub struct Isovolume<const D: usize> {
17    origin: [f32; D],
18    cell: f32,
19    boxes: Vec<CellBox<D>>,
20    occupied_cells: usize,
21    clipped: bool,
22}
23
24impl<const D: usize> Isovolume<D> {
25    /// Requires a 1-Lipschitz signed distance; `resolution` counts cells along the longest axis.
26    pub fn extract(
27        min: [f32; D],
28        max: [f32; D],
29        resolution: usize,
30        sdf: impl Fn([f32; D]) -> f32,
31    ) -> Self {
32        assert!(resolution > 0, "resolution must be positive");
33        let mut extent = [0.0f32; D];
34        for (e, (lo, hi)) in extent.iter_mut().zip(min.iter().zip(max.iter())) {
35            *e = hi - lo;
36            assert!(*e > 0.0, "domain must be non-degenerate on every axis");
37        }
38
39        let cell = extent.iter().copied().fold(0.0f32, f32::max) / resolution as f32;
40        let mut counts = [0usize; D];
41        for (n, e) in counts.iter_mut().zip(extent.iter()) {
42            *n = ((e / cell).ceil() as usize).max(1);
43        }
44
45        let margin = 0.5 * cell * (D as f32).sqrt();
46        let total: usize = counts.iter().product();
47        let mut occupied = vec![false; total];
48        let mut occupied_cells = 0usize;
49        let mut clipped = false;
50        for (index, marked) in occupied.iter_mut().enumerate() {
51            let coords = unflatten(index, &counts);
52            if sdf(cell_centre(&min, cell, &coords)) > margin {
53                continue;
54            }
55            *marked = true;
56            occupied_cells += 1;
57            clipped |= coords
58                .iter()
59                .zip(counts.iter())
60                .any(|(&c, &n)| c == 0 || c + 1 == n);
61        }
62
63        let mut covered = vec![false; total];
64        let mut boxes = Vec::new();
65        let mut seed = 0usize;
66        loop {
67            while seed < total && (covered[seed] || !occupied[seed]) {
68                seed += 1;
69            }
70            if seed == total {
71                break;
72            }
73            let grown = grow(unflatten(seed, &counts), &occupied, &counts);
74            visit_cells(&grown, |coords| {
75                covered[flatten(coords, &counts)] = true;
76                true
77            });
78            boxes.push(grown);
79        }
80
81        Self {
82            origin: min,
83            cell,
84            boxes,
85            occupied_cells,
86            clipped,
87        }
88    }
89
90    pub fn piece_count(&self) -> usize {
91        self.boxes.len()
92    }
93
94    pub fn occupied_cells(&self) -> usize {
95        self.occupied_cells
96    }
97
98    pub fn cell_size(&self) -> f32 {
99        self.cell
100    }
101
102    /// Upper bound on how far the cover extends past the true surface.
103    pub fn enclosure_margin(&self) -> f32 {
104        self.cell * (D as f32).sqrt()
105    }
106
107    /// Exact volume of the union of the pieces: that union is the marked cell set.
108    pub fn volume(&self) -> f32 {
109        self.occupied_cells as f32 * self.cell.powi(D as i32)
110    }
111
112    /// True when a marked cell touched the sampled domain's boundary.
113    pub fn clipped(&self) -> bool {
114        self.clipped
115    }
116
117    /// Corner-to-corner, in the grid frame.
118    pub fn piece_bounds(&self, index: usize) -> ([f32; D], [f32; D]) {
119        let b = self.boxes[index];
120        let mut lo = [0.0f32; D];
121        let mut hi = [0.0f32; D];
122        for k in 0..D {
123            lo[k] = self.origin[k] + b.lo[k] as f32 * self.cell;
124            hi[k] = self.origin[k] + (b.hi[k] + 1) as f32 * self.cell;
125        }
126        (lo, hi)
127    }
128
129    /// Linear in piece count: a verification query, not a broadphase.
130    pub fn contains(&self, p: [f32; D]) -> bool {
131        (0..self.boxes.len()).any(|i| {
132            let (lo, hi) = self.piece_bounds(i);
133            (0..D).all(|k| p[k] >= lo[k] && p[k] <= hi[k])
134        })
135    }
136}
137
138impl Isovolume<3> {
139    /// Pose is extrinsic per the [`Shape`] contract, so the hull is origin-centred.
140    pub fn colliders(&self) -> Vec<(Vec3, Shape)> {
141        (0..self.boxes.len())
142            .map(|i| {
143                let (lo, hi) = self.piece_bounds(i);
144                let centre = Vec3::from_array(lo).lerp(Vec3::from_array(hi), 0.5);
145                let h = 0.5 * (Vec3::from_array(hi) - Vec3::from_array(lo));
146                let mut vertices = Vec::with_capacity(8);
147                for sz in [-1.0f32, 1.0] {
148                    for sy in [-1.0f32, 1.0] {
149                        for sx in [-1.0f32, 1.0] {
150                            vertices.push(Vec3::new(sx * h.x, sy * h.y, sz * h.z));
151                        }
152                    }
153                }
154                (centre, Shape::ConvexPolytope3D { vertices })
155            })
156            .collect()
157    }
158}
159
160impl Isovolume<4> {
161    pub fn colliders(&self) -> Vec<(Vec4, Shape)> {
162        (0..self.boxes.len())
163            .map(|i| {
164                let (lo, hi) = self.piece_bounds(i);
165                let centre = Vec4::from_array(lo).lerp(Vec4::from_array(hi), 0.5);
166                let h = 0.5 * (Vec4::from_array(hi) - Vec4::from_array(lo));
167                let mut vertices = Vec::with_capacity(16);
168                for sw in [-1.0f32, 1.0] {
169                    for sz in [-1.0f32, 1.0] {
170                        for sy in [-1.0f32, 1.0] {
171                            for sx in [-1.0f32, 1.0] {
172                                vertices.push(Vec4::new(sx * h.x, sy * h.y, sz * h.z, sw * h.w));
173                            }
174                        }
175                    }
176                }
177                (centre, Shape::ConvexPolytope4D { vertices })
178            })
179            .collect()
180    }
181}
182
183fn cell_centre<const D: usize>(origin: &[f32; D], cell: f32, coords: &[usize; D]) -> [f32; D] {
184    let mut p = [0.0f32; D];
185    for k in 0..D {
186        p[k] = origin[k] + (coords[k] as f32 + 0.5) * cell;
187    }
188    p
189}
190
191fn flatten<const D: usize>(coords: &[usize; D], counts: &[usize; D]) -> usize {
192    let mut index = 0;
193    for k in (0..D).rev() {
194        index = index * counts[k] + coords[k];
195    }
196    index
197}
198
199fn unflatten<const D: usize>(mut index: usize, counts: &[usize; D]) -> [usize; D] {
200    let mut coords = [0usize; D];
201    for k in 0..D {
202        coords[k] = index % counts[k];
203        index /= counts[k];
204    }
205    coords
206}
207
208fn visit_cells<const D: usize>(b: &CellBox<D>, mut f: impl FnMut(&[usize; D]) -> bool) {
209    let mut coords = b.lo;
210    loop {
211        if !f(&coords) {
212            return;
213        }
214        let mut axis = 0;
215        loop {
216            if axis == D {
217                return;
218            }
219            if coords[axis] < b.hi[axis] {
220                coords[axis] += 1;
221                break;
222            }
223            coords[axis] = b.lo[axis];
224            axis += 1;
225        }
226    }
227}
228
229fn all_occupied<const D: usize>(b: &CellBox<D>, occupied: &[bool], counts: &[usize; D]) -> bool {
230    let mut ok = true;
231    visit_cells(b, |coords| {
232        ok = occupied[flatten(coords, counts)];
233        ok
234    });
235    ok
236}
237
238fn grow<const D: usize>(seed: [usize; D], occupied: &[bool], counts: &[usize; D]) -> CellBox<D> {
239    let mut b = CellBox { lo: seed, hi: seed };
240    loop {
241        let mut grew = false;
242        for axis in 0..D {
243            if b.lo[axis] > 0 {
244                let mut layer = b;
245                layer.lo[axis] = b.lo[axis] - 1;
246                layer.hi[axis] = b.lo[axis] - 1;
247                if all_occupied(&layer, occupied, counts) {
248                    b.lo[axis] -= 1;
249                    grew = true;
250                }
251            }
252            if b.hi[axis] + 1 < counts[axis] {
253                let mut layer = b;
254                layer.lo[axis] = b.hi[axis] + 1;
255                layer.hi[axis] = b.hi[axis] + 1;
256                if all_occupied(&layer, occupied, counts) {
257                    b.hi[axis] += 1;
258                    grew = true;
259                }
260            }
261        }
262        if !grew {
263            return b;
264        }
265    }
266}
267
268#[cfg(test)]
269mod tests {
270    use super::*;
271
272    fn torus_3d(major: f32, minor: f32) -> impl Fn([f32; 3]) -> f32 {
273        move |p| {
274            let radial = (p[0] * p[0] + p[2] * p[2]).sqrt() - major;
275            (radial * radial + p[1] * p[1]).sqrt() - minor
276        }
277    }
278
279    fn torus_4d(major: f32, minor: f32) -> impl Fn([f32; 4]) -> f32 {
280        move |p| {
281            let radial = (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]).sqrt() - major;
282            (radial * radial + p[3] * p[3]).sqrt() - minor
283        }
284    }
285
286    const TORUS_MAJOR: f32 = 1.0;
287    const TORUS_MINOR: f32 = 0.3;
288    const TORUS_BOUND: f32 = 1.6;
289
290    fn torus_volume_3d() -> f32 {
291        2.0 * std::f32::consts::PI.powi(2) * TORUS_MAJOR * TORUS_MINOR * TORUS_MINOR
292    }
293
294    fn extract_torus_3d(resolution: usize) -> Isovolume<3> {
295        Isovolume::extract(
296            [-TORUS_BOUND; 3],
297            [TORUS_BOUND; 3],
298            resolution,
299            torus_3d(TORUS_MAJOR, TORUS_MINOR),
300        )
301    }
302
303    fn direction_3d(i: usize) -> [f32; 3] {
304        const ALPHA_1: f32 = 0.754_877_7;
305        const ALPHA_2: f32 = 0.569_840_3;
306        let u = (0.5 + ALPHA_1 * i as f32).fract();
307        let v = (0.5 + ALPHA_2 * i as f32).fract();
308        let z = 2.0 * u - 1.0;
309        let r = (1.0 - z * z).max(0.0).sqrt();
310        let phi = 2.0 * std::f32::consts::PI * v;
311        [r * phi.cos(), r * phi.sin(), z]
312    }
313
314    fn surface_point_3d(
315        sdf: &impl Fn([f32; 3]) -> f32,
316        from: [f32; 3],
317        dir: [f32; 3],
318    ) -> Option<[f32; 3]> {
319        const STEP: f32 = 0.01;
320        const STEPS: usize = 800;
321        let at = |t: f32| {
322            [
323                from[0] + dir[0] * t,
324                from[1] + dir[1] * t,
325                from[2] + dir[2] * t,
326            ]
327        };
328        if sdf(at(0.0)) > 0.0 {
329            return None;
330        }
331        let mut lo = 0.0f32;
332        let hi = (1..=STEPS).map(|i| i as f32 * STEP).find(|&t| {
333            if sdf(at(t)) <= 0.0 {
334                lo = t;
335                false
336            } else {
337                true
338            }
339        })?;
340        let mut hi = hi;
341        for _ in 0..40 {
342            let mid = 0.5 * (lo + hi);
343            if sdf(at(mid)) <= 0.0 {
344                lo = mid;
345            } else {
346                hi = mid;
347            }
348        }
349        Some(at(0.5 * (lo + hi)))
350    }
351
352    #[test]
353    fn rectangular_occupied_region_merges_into_one_piece() {
354        fn check<const D: usize>() {
355            let volume = Isovolume::extract([-2.0; D], [2.0; D], 16, |p| {
356                p.into_iter().map(f32::abs).fold(0.0, f32::max) - 0.5
357            });
358            assert_eq!(volume.piece_count(), 1);
359            assert_eq!(volume.piece_bounds(0), ([-0.75; D], [0.75; D]));
360            assert_eq!(volume.occupied_cells(), 6usize.pow(D as u32));
361        }
362        check::<3>();
363        check::<4>();
364    }
365
366    #[test]
367    fn every_isosurface_sample_lies_inside_a_piece() {
368        let sdf = torus_3d(TORUS_MAJOR, TORUS_MINOR);
369        for resolution in [24, 32, 48] {
370            let volume = extract_torus_3d(resolution);
371            assert!(!volume.clipped());
372            let mut checked = 0;
373            for i in 0..512 {
374                let angle = i as f32 * std::f32::consts::TAU / 512.0;
375                let core = [TORUS_MAJOR * angle.cos(), 0.0, TORUS_MAJOR * angle.sin()];
376                let dir = direction_3d(i);
377                let Some(p) = surface_point_3d(&sdf, core, dir) else {
378                    continue;
379                };
380                assert!(
381                    sdf(p).abs() < 1.0e-4,
382                    "bisection did not land on the surface: f = {}",
383                    sdf(p)
384                );
385                assert!(
386                    volume.contains(p),
387                    "resolution {resolution}: surface point {p:?} is outside every piece"
388                );
389                checked += 1;
390            }
391            assert!(checked > 400, "only {checked} usable surface samples");
392        }
393    }
394
395    #[test]
396    fn every_interior_probe_lies_inside_a_piece() {
397        let sdf = torus_3d(TORUS_MAJOR, TORUS_MINOR);
398        let volume = extract_torus_3d(32);
399        const PROBES: usize = 61;
400        let mut interior = 0;
401        for iz in 0..PROBES {
402            for iy in 0..PROBES {
403                for ix in 0..PROBES {
404                    let coord = |i: usize| {
405                        -TORUS_BOUND + 2.0 * TORUS_BOUND * (i as f32 + 0.31) / PROBES as f32
406                    };
407                    let p = [coord(ix), coord(iy), coord(iz)];
408                    if sdf(p) > 0.0 {
409                        continue;
410                    }
411                    interior += 1;
412                    assert!(volume.contains(p), "interior probe {p:?} is not covered");
413                }
414            }
415        }
416        assert!(interior > 1000, "only {interior} interior probes");
417    }
418
419    #[test]
420    fn no_piece_corner_exceeds_the_enclosure_margin() {
421        let sdf = torus_3d(TORUS_MAJOR, TORUS_MINOR);
422        let volume = extract_torus_3d(32);
423        let margin = volume.enclosure_margin();
424        for i in 0..volume.piece_count() {
425            let (lo, hi) = volume.piece_bounds(i);
426            for cz in [lo[2], hi[2]] {
427                for cy in [lo[1], hi[1]] {
428                    for cx in [lo[0], hi[0]] {
429                        let d = sdf([cx, cy, cz]);
430                        assert!(d <= margin + 1.0e-5, "corner distance {d} exceeds {margin}");
431                    }
432                }
433            }
434        }
435    }
436
437    #[test]
438    fn cover_volume_tightens_towards_the_true_volume() {
439        let truth = torus_volume_3d();
440        let mut previous = f32::INFINITY;
441        for resolution in [16, 32, 64] {
442            let ratio = extract_torus_3d(resolution).volume() / truth;
443            assert!(
444                ratio > 1.0,
445                "resolution {resolution}: cover {ratio} misses volume"
446            );
447            assert!(
448                ratio < previous,
449                "resolution {resolution}: ratio {ratio} did not improve on {previous}"
450            );
451            previous = ratio;
452        }
453        assert!(
454            previous < 1.6,
455            "finest cover is still {previous}x the solid"
456        );
457    }
458
459    #[test]
460    fn bounds_that_cut_the_solid_are_reported_as_clipped() {
461        let tight = Isovolume::extract([-0.8; 3], [0.8; 3], 24, torus_3d(TORUS_MAJOR, TORUS_MINOR));
462        assert!(tight.clipped());
463        assert!(!extract_torus_3d(24).clipped());
464    }
465
466    #[test]
467    fn a_field_with_no_solid_yields_no_pieces() {
468        let empty = Isovolume::extract([-1.0; 3], [1.0; 3], 8, |_| 10.0);
469        assert_eq!(empty.piece_count(), 0);
470        assert_eq!(empty.occupied_cells(), 0);
471        assert_eq!(empty.volume(), 0.0);
472        assert!(empty.colliders().is_empty());
473        assert!(!empty.contains([0.0; 3]));
474    }
475
476    #[test]
477    fn colliders_are_origin_centred_hulls_matching_their_piece_bounds() {
478        let volume = extract_torus_3d(16);
479        let colliders = volume.colliders();
480        assert_eq!(colliders.len(), volume.piece_count());
481        for (i, (centre, shape)) in colliders.iter().enumerate() {
482            let Shape::ConvexPolytope3D { vertices } = shape else {
483                panic!("expected ConvexPolytope3D, got {:?}", shape.kind());
484            };
485            assert_eq!(vertices.len(), 8);
486            let (lo, hi) = volume.piece_bounds(i);
487            let expect = 0.5 * (Vec3::from_array(lo) + Vec3::from_array(hi));
488            assert!(centre.distance(expect) < 1.0e-5);
489            let sum: Vec3 = vertices.iter().copied().sum();
490            assert!(sum.length() < 1.0e-5, "hull is not origin-centred");
491            let half = 0.5 * (Vec3::from_array(hi) - Vec3::from_array(lo));
492            for v in vertices {
493                assert!((v.abs() - half).abs().max_element() < 1.0e-5);
494            }
495        }
496    }
497
498    #[test]
499    fn colliders_4d_are_origin_centred_sixteen_vertex_hulls() {
500        let volume = Isovolume::extract(
501            [-TORUS_BOUND; 4],
502            [TORUS_BOUND; 4],
503            10,
504            torus_4d(TORUS_MAJOR, TORUS_MINOR),
505        );
506        let colliders = volume.colliders();
507        assert_eq!(colliders.len(), volume.piece_count());
508        assert!(volume.piece_count() > 0);
509        for (i, (centre, shape)) in colliders.iter().enumerate() {
510            let Shape::ConvexPolytope4D { vertices } = shape else {
511                panic!("expected ConvexPolytope4D, got {:?}", shape.kind());
512            };
513            assert_eq!(vertices.len(), 16);
514            let (lo, hi) = volume.piece_bounds(i);
515            let expect = 0.5 * (Vec4::from_array(lo) + Vec4::from_array(hi));
516            assert!(centre.distance(expect) < 1.0e-5);
517            let half = 0.5 * (Vec4::from_array(hi) - Vec4::from_array(lo));
518            for v in vertices {
519                assert!((v.abs() - half).abs().max_element() < 1.0e-5);
520            }
521        }
522    }
523
524    #[test]
525    fn every_interior_probe_lies_inside_a_piece_in_4d() {
526        let sdf = torus_4d(TORUS_MAJOR, TORUS_MINOR);
527        let volume = Isovolume::extract([-TORUS_BOUND; 4], [TORUS_BOUND; 4], 16, &sdf);
528        assert!(!volume.clipped());
529        const PROBES: usize = 21;
530        let coord = |i: usize| -TORUS_BOUND + 2.0 * TORUS_BOUND * (i as f32 + 0.37) / PROBES as f32;
531        let mut interior = 0;
532        for iw in 0..PROBES {
533            for iz in 0..PROBES {
534                for iy in 0..PROBES {
535                    for ix in 0..PROBES {
536                        let p = [coord(ix), coord(iy), coord(iz), coord(iw)];
537                        if sdf(p) > 0.0 {
538                            continue;
539                        }
540                        interior += 1;
541                        assert!(volume.contains(p), "interior probe {p:?} is not covered");
542                    }
543                }
544            }
545        }
546        assert!(interior > 1000, "only {interior} interior probes");
547    }
548}
loam_shape/isovolume.rs

loam_shape/
isovolume.rs
