fast_mesh_fields/
thickness.rs

1//! Per-vertex wall thickness of a mesh.
2//!
3//! [`thicknesses`] writes, for every vertex slot, the distance to the
4//! opposite surface along the vertex's smooth normal: the nearest ray hit in
5//! either direction (`+n` or `-n`), in the mesh's local `f32` coordinates.
6//! Closed solids read their wall thickness (a slab's faces read its depth, a
7//! sphere's vertices its diameter); a ray that escapes in both directions
8//! reads NaN, so open sheets and the rims of open meshes report NaN rather
9//! than a distance to nowhere. Hits beyond `max_distance` read `+∞`, the
10//! same "far" value the distance kernels use (and so the same [`ColorRule`]
11//! handling: `above` for far, `nan` for open). Dead slots and alive vertices
12//! no alive triangle uses read NaN.
13//!
14//! The normals are area-weighted face-normal sums, normalised (zero stays
15//! zero, and reads NaN). Self-hits are skipped by starting the rays a small
16//! way off the surface, scaled to the mesh's coordinate magnitude. The sign
17//! of the normal does not matter (both directions are tried), so meshes with
18//! flipped patches still read the nearest opposite surface.
19//!
20//! Nothing is allocated per vertex. With the `parallel` feature, large
21//! batches of rays run in chunks on rayon's pool through
22//! [`fast_mesh_query::raycast_many`]; the values are the same, bit for bit.
23
24use fast_mesh_bvh::Bvh;
25use fast_mesh_query::{Cull, RayOptions, raycast_many};
26use fast_mesh_store::Mesh;
27
28use crate::FieldsError;
29use crate::Result;
30
31/// A thickness that is this fraction of the mesh's coordinate magnitude, or
32/// less, is skipped as a self-hit: `32 * f32::EPSILON` times the largest
33/// absolute coordinate, about twice the rounding of one ray-triangle test.
34fn self_gap(magnitude: f32) -> f32 {
35    (32.0 * f32::EPSILON * magnitude).max(f32::MIN_POSITIVE)
36}
37
38/// The largest absolute coordinate of the mesh's positions, for
39/// [`self_gap`]. 0 when the mesh has no vertices.
40fn coordinate_magnitude(mesh: &Mesh) -> f32 {
41    mesh.positions()
42        .iter()
43        .flatten()
44        .map(|c| c.abs())
45        .fold(0.0, f32::max)
46}
47
48/// Per-vertex wall thickness of `mesh`: `out[v]` gets the value for vertex
49/// slot `v` (see the module docs).
50///
51/// `max_distance` bounds the values kept (`+∞` beyond it, inclusive); it
52/// must be zero or greater and not NaN, and may be infinite to keep every
53/// hit. `out` must hold one entry per vertex slot.
54///
55/// # Errors
56/// Everything is checked before anything is written, so on error `out` is
57/// unchanged:
58/// - [`FieldsError::InvalidInput`] when `out` is not one entry per vertex
59///   slot, or `max_distance` is negative or NaN;
60/// - [`FieldsError::Query`] with `MeshMismatch` or `StaleBvh` when `bvh`
61///   does not match the current `mesh`.
62pub fn thicknesses(mesh: &Mesh, bvh: &Bvh, max_distance: f32, out: &mut [f32]) -> Result<()> {
63    if out.len() != mesh.vertex_slots() {
64        return Err(FieldsError::InvalidInput(
65            "out must hold one entry per vertex slot",
66        ));
67    }
68    if max_distance.is_nan() || max_distance < 0.0 {
69        return Err(FieldsError::InvalidInput(
70            "max_distance must be zero or greater, and not NaN",
71        ));
72    }
73    if bvh.mesh_uid() != mesh.uid() {
74        return Err(fast_mesh_query::QueryError::MeshMismatch.into());
75    }
76    if !bvh.is_current(mesh) {
77        return Err(fast_mesh_query::QueryError::StaleBvh.into());
78    }
79    let slots = mesh.vertex_slots();
80    let mut normals = vec![[0.0f32; 3]; slots];
81    // Area-weighted face-normal sums, over the alive triangles in ascending
82    // order (the same sums `fast_mesh::compute_vertex_normals` makes).
83    for t in 0..mesh.triangle_slots() as u32 {
84        let Some(tri) = mesh.triangle(t) else {
85            continue;
86        };
87        if !mesh.is_triangle_alive(t) {
88            continue;
89        }
90        let p = mesh.positions();
91        let [a, b, c] = tri.map(|v| p[v as usize]);
92        let e1 = [b[0] - a[0], b[1] - a[1], b[2] - a[2]];
93        let e2 = [c[0] - a[0], c[1] - a[1], c[2] - a[2]];
94        let n = [
95            e1[1] * e2[2] - e1[2] * e2[1],
96            e1[2] * e2[0] - e1[0] * e2[2],
97            e1[0] * e2[1] - e1[1] * e2[0],
98        ];
99        for v in tri {
100            for k in 0..3 {
101                normals[v as usize][k] += n[k];
102            }
103        }
104    }
105    // Unit normals; zero stays zero (degenerate: reads NaN below).
106    for n in &mut normals {
107        let len = (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
108        if len > 0.0 {
109            for c in n.iter_mut() {
110                *c /= len;
111            }
112        }
113    }
114    let gap = self_gap(coordinate_magnitude(mesh));
115    let opts = RayOptions {
116        min_t: gap,
117        max_t: f32::INFINITY,
118        cull: Cull::None,
119    };
120    // One ray pair per defined vertex with a normal; the rest read NaN.
121    let mut origins: Vec<[f32; 3]> = Vec::new();
122    let mut dirs: Vec<[f32; 3]> = Vec::new();
123    let mut slots_of: Vec<u32> = Vec::new();
124    for v in 0..slots as u32 {
125        let n = normals[v as usize];
126        let defined = mesh.is_vertex_alive(v)
127            && mesh.corner_count(v) > 0
128            && n != [0.0; 3]
129            && mesh
130                .position(v)
131                .is_some_and(|p| p.iter().all(|c| c.is_finite()));
132        if !defined {
133            continue;
134        }
135        let p = mesh.position(v).unwrap_or([0.0; 3]);
136        origins.push(p);
137        dirs.push(n);
138        origins.push(p);
139        dirs.push([-n[0], -n[1], -n[2]]);
140        slots_of.push(v);
141        slots_of.push(v);
142    }
143    let mut hits = vec![None; origins.len()];
144    if !origins.is_empty() {
145        raycast_many(mesh, bvh, &origins, &dirs, &opts, &mut hits)?;
146    }
147    for v in out.iter_mut() {
148        *v = f32::NAN;
149    }
150    debug_assert_eq!(hits.len(), slots_of.len());
151    for (i, &v) in slots_of.iter().enumerate() {
152        let t = hits[i].as_ref().map(|h| h.t);
153        let cur = out[v as usize];
154        match (cur.is_nan(), t) {
155            // First of the pair, or a miss after a hit: keep what is there.
156            (_, None) => {}
157            // A hit where none was: the distance, far when past the bound.
158            (true, Some(t)) => {
159                out[v as usize] = if t > max_distance { f32::INFINITY } else { t };
160            }
161            // The nearer of the pair, far when past the bound.
162            (false, Some(t)) => {
163                let got = if t > max_distance { f32::INFINITY } else { t };
164                if got < cur {
165                    out[v as usize] = got;
166                }
167            }
168        }
169    }
170    // Keep the parallel paths warm under the same thresholds the distance
171    // kernels use (a no-op here: `raycast_many` parallelises on its own).
172    Ok(())
173}
174
175#[cfg(test)]
176mod tests {
177    use fast_mesh_bvh::{BuildOptions, Bvh};
178    use fast_mesh_store::Mesh;
179
180    use super::*;
181    use crate::test_util::*;
182    use fast_mesh_store::par;
183
184    fn build(mesh: &Mesh) -> Bvh {
185        Bvh::build(
186            mesh,
187            &BuildOptions {
188                max_leaf_size: 2,
189                ..Default::default()
190            },
191        )
192        .unwrap()
193    }
194
195    fn run(mesh: &Mesh, max_distance: f32) -> Vec<f32> {
196        let bvh = build(mesh);
197        let mut out = vec![f32::NAN; mesh.vertex_slots()];
198        thicknesses(mesh, &bvh, max_distance, &mut out).unwrap();
199        out
200    }
201
202    #[test]
203    fn cube_vertices_read_wall_to_wall() {
204        // Every corner of the [-1, 1]^3 cube reads the distance to the
205        // opposite surface along its area-weighted normal: between one face
206        // width (2) and the body diagonal (2 sqrt(3)).
207        let mesh = cube();
208        let out = run(&mesh, f32::INFINITY);
209        assert_eq!(out.len(), 8);
210        for &v in &out {
211            assert!(v.is_finite(), "{out:?}");
212            assert!((2.0..=2.0 * 3f32.sqrt() + 1e-4).contains(&v), "{out:?}");
213        }
214    }
215
216    /// A closed plate, 1 x 1 x 0.5 (`n` cells a side), facing out: the
217    /// interior face vertices have untilted normals and read the depth.
218    fn plate(n: u32) -> Mesh {
219        let h = 0.25;
220        let w = n + 1;
221        let mut positions = Vec::new();
222        for z in [h, -h] {
223            for y in 0..=n {
224                for x in 0..=n {
225                    positions.push([x as f32 / n as f32, y as f32 / n as f32, z]);
226                }
227            }
228        }
229        let (top, bottom) = (0, w * w);
230        let t = |x: u32, y: u32| top + y * w + x;
231        let b = |x: u32, y: u32| bottom + y * w + x;
232        let mut tris = Vec::new();
233        for y in 0..n {
234            for x in 0..n {
235                tris.extend([
236                    [t(x, y), t(x + 1, y), t(x + 1, y + 1)],
237                    [t(x, y), t(x + 1, y + 1), t(x, y + 1)],
238                ]);
239                tris.extend([
240                    [b(x, y), b(x + 1, y + 1), b(x + 1, y)],
241                    [b(x, y), b(x, y + 1), b(x + 1, y + 1)],
242                ]);
243            }
244        }
245        for x in 0..n {
246            // South (y = 0) and north (y = 1) skirts.
247            tris.extend([
248                [t(x, 0), b(x, 0), b(x + 1, 0)],
249                [t(x, 0), b(x + 1, 0), t(x + 1, 0)],
250                [t(x, n), t(x + 1, n), b(x + 1, n)],
251                [t(x, n), b(x + 1, n), b(x, n)],
252            ]);
253            // West (x = 0) and east (x = 1) skirts.
254            tris.extend([
255                [t(0, x), t(0, x + 1), b(0, x + 1)],
256                [t(0, x), b(0, x + 1), b(0, x)],
257                [t(n, x), b(n, x), b(n, x + 1)],
258                [t(n, x), b(n, x + 1), t(n, x + 1)],
259            ]);
260        }
261        mesh(&positions, &tris)
262    }
263
264    #[test]
265    fn slab_faces_read_depth() {
266        // The interior face vertices read the plate depth, 0.5.
267        let m = plate(4);
268        let out = run(&m, f32::INFINITY);
269        let w = 5;
270        for y in 1..4 {
271            for x in 1..4 {
272                for &v in &[y * w + x, w * w + y * w + x] {
273                    assert!((out[v as usize] - 0.5).abs() < 1e-3, "{out:?}");
274                }
275            }
276        }
277        assert!(out.iter().all(|v| v.is_finite()), "{out:?}");
278    }
279
280    #[test]
281    fn sphere_reads_diameter() {
282        let mesh = icosphere();
283        let out = run(&mesh, f32::INFINITY);
284        assert_eq!(out.len(), mesh.vertex_slots());
285        for &v in &out {
286            assert!(v.is_finite(), "{out:?}");
287            assert!((v - 2.0).abs() < 0.05, "{v}");
288        }
289    }
290
291    #[test]
292    fn open_sheet_reads_nan() {
293        // An open sheet has no opposite surface: every ray escapes.
294        let mesh = terrain(4);
295        let out = run(&mesh, f32::INFINITY);
296        assert!(!out.is_empty());
297        for &v in &out {
298            assert!(v.is_nan(), "{v}");
299        }
300    }
301
302    #[test]
303    fn open_cube_rim_reads_nan_while_walls_read_finite() {
304        // The open-cube fixture (top face removed): the missing face's
305        // vertices have no opposite surface above them... check the mix:
306        // walls read finite, and nothing reads a bogus zero.
307        let (p, mut t) = cube_buffers();
308        t.drain(2..4);
309        let m = mesh(&p, &t);
310        let out = run(&m, f32::INFINITY);
311        assert_eq!(out.len(), 8);
312        for &v in &out {
313            assert!(v.is_nan() || v > 1.0, "{out:?}");
314        }
315        assert!(out.iter().any(|v| v.is_finite()), "{out:?}");
316    }
317
318    #[test]
319    fn far_reads_infinite() {
320        let mesh = cube();
321        let out = run(&mesh, 1.0);
322        for &v in &out {
323            assert_eq!(v, f32::INFINITY, "{out:?}");
324        }
325        let out = run(&mesh, 4.0);
326        for &v in &out {
327            assert!(v.is_finite(), "{out:?}");
328        }
329    }
330
331    #[test]
332    fn dead_and_unused_slots_read_nan() {
333        // An isolated vertex (alive, used by no triangle) reads NaN while
334        // the cube around it reads finite.
335        let (mut positions, tris) = cube_buffers();
336        positions.push([5.0, 5.0, 5.0]);
337        let m = mesh(&positions, &tris);
338        assert_eq!(m.vertex_slots(), 9);
339        assert!(m.is_vertex_alive(8));
340        assert_eq!(m.corner_count(8), 0);
341        let bvh = build(&m);
342        let mut out = vec![7.0f32; m.vertex_slots()];
343        thicknesses(&m, &bvh, f32::INFINITY, &mut out).unwrap();
344        assert!(out[8].is_nan(), "{out:?}");
345        assert!(out[..8].iter().all(|v| v.is_finite()), "{out:?}");
346    }
347
348    #[test]
349    fn invalid_input_writes_nothing() {
350        let mesh = cube();
351        let bvh = build(&mesh);
352        let mut out = vec![7.0f32; mesh.vertex_slots()];
353        let bad_len = thicknesses(&mesh, &bvh, f32::INFINITY, &mut out[..4]);
354        assert!(matches!(bad_len, Err(FieldsError::InvalidInput(_))));
355        assert!(matches!(
356            thicknesses(&mesh, &bvh, f32::NAN, &mut out),
357            Err(FieldsError::InvalidInput(_))
358        ));
359        assert!(matches!(
360            thicknesses(&mesh, &bvh, -1.0, &mut out),
361            Err(FieldsError::InvalidInput(_))
362        ));
363        assert_eq!(out, vec![7.0; mesh.vertex_slots()]);
364        let other = cube();
365        let other_bvh = build(&other);
366        assert!(matches!(
367            thicknesses(&mesh, &other_bvh, 1.0, &mut out),
368            Err(FieldsError::Query(_))
369        ));
370        assert_eq!(out, vec![7.0; mesh.vertex_slots()]);
371    }
372
373    #[test]
374    fn thickness_bit_identical() {
375        // The forced-sequential run is the code a feature-off build runs.
376        // The torus is closed (every vertex reads finite), with enough rays
377        // to pass the parallel threshold.
378        let mesh = torus();
379        let bvh = build(&mesh);
380        assert!(mesh.triangle_count() >= 500);
381        let run = || {
382            let mut out = vec![f32::NAN; mesh.vertex_slots()];
383            thicknesses(&mesh, &bvh, f32::INFINITY, &mut out).unwrap();
384            out.iter().map(|x| x.to_bits()).collect::<Vec<_>>()
385        };
386        par::set_force_sequential(true);
387        let reference = run();
388        par::set_force_sequential(false);
389        if par::enabled() {
390            for k in [1usize, 2, 8] {
391                assert_eq!(par::with_threads(k, run), reference, "{k} threads");
392            }
393        }
394        // Every vertex of the torus reads a finite wall thickness.
395        let mut out = vec![f32::NAN; mesh.vertex_slots()];
396        thicknesses(&mesh, &bvh, f32::INFINITY, &mut out).unwrap();
397        assert!(out.iter().all(|v| v.is_finite()), "{out:?}");
398    }
399}
fast_mesh_fields/thickness.rs

fast_mesh_fields/
thickness.rs
