fast_mesh_fields/
distance.rs

1//! Signed distances from points to a mesh.
2
3use fast_mesh_bvh::Bvh;
4use fast_mesh_query::{ClosestHit, ClosestOptions, QueryError, closest_point, closest_point_many};
5use fast_mesh_store::{Mesh, par};
6
7use crate::{FieldsError, Result, SignData};
8
9/// Points from which [`signed_distances`] runs in parallel chunks with the
10/// `parallel` feature. A point costs about a microsecond (more in winding
11/// mode). Measured on the i5-13600K against a 1M-triangle mesh: 256 points
12/// break even (0.18 ms either way), 512 run 2.3 times faster in parallel.
13pub(crate) const PAR_MIN_POINTS: usize = 512;
14
15/// Points per parallel chunk of [`signed_distances`].
16pub(crate) const POINT_CHUNK: usize = 128;
17
18/// Signed distances from `points_in_b` (in mesh B's local coordinates) to
19/// mesh B: `out[i]` gets the value for `points_in_b[i]`.
20///
21/// - The magnitude is the distance to the closest point of B, from
22///   [`fast_mesh_query::closest_point`].
23/// - The sign comes from `sign` (see [`SignData`]): negative behind the
24///   surface (inside), positive in front of it (outside). A point exactly on
25///   the surface gets `+0.0`; a point whose side is undecided (exactly in the
26///   tangent plane of its closest feature) counts as outside.
27/// - A point farther than `max_distance` from the surface (inclusive bound)
28///   gets `f32::INFINITY`, meaning far, with no sign. That includes points
29///   deeper inside than `max_distance`. An infinite `max_distance` computes
30///   every point.
31///
32/// Nothing is allocated per point. With the `parallel` feature, large
33/// batches run in chunks on rayon's pool, each chunk into its own buffer;
34/// the values are the same, bit for bit.
35///
36/// # Errors
37/// Everything is checked before anything is written, so on error `out` is
38/// unchanged:
39/// - [`FieldsError::InvalidInput`] when the slices differ in length, a point
40///   is not finite (NaN or infinite), or `max_distance` is negative or NaN;
41/// - [`FieldsError::Query`] with `MeshMismatch` or `StaleBvh` when `bvh_b`
42///   does not match the current `mesh_b`;
43/// - [`FieldsError::InvalidInput`] when `sign` is not current for `mesh_b`
44///   and `bvh_b` (see [`SignData::is_current`]).
45pub fn signed_distances(
46    mesh_b: &Mesh,
47    bvh_b: &Bvh,
48    sign: &SignData,
49    points_in_b: &[[f32; 3]],
50    max_distance: f32,
51    out: &mut [f32],
52) -> Result<()> {
53    if points_in_b.len() != out.len() {
54        return Err(FieldsError::InvalidInput(
55            "points and out must have the same length",
56        ));
57    }
58    if max_distance.is_nan() || max_distance < 0.0 {
59        return Err(FieldsError::InvalidInput(
60            "max_distance must be zero or greater, and not NaN",
61        ));
62    }
63    if !points_in_b.iter().flatten().all(|c| c.is_finite()) {
64        return Err(FieldsError::InvalidInput("points must be finite"));
65    }
66    if bvh_b.mesh_uid() != mesh_b.uid() {
67        return Err(QueryError::MeshMismatch.into());
68    }
69    if !bvh_b.is_current(mesh_b) {
70        return Err(QueryError::StaleBvh.into());
71    }
72    if !sign.is_current(mesh_b, bvh_b) {
73        return Err(FieldsError::InvalidInput(
74            "sign data is not current for this mesh and BVH",
75        ));
76    }
77    // The closest points, validated once per batch: `closest_point_many`
78    // checks its inputs a single time, where per-point `closest_point` calls
79    // revalidate the BVH on every point. The hits are identical either way.
80    let opts = ClosestOptions { max_distance };
81    let sequential =
82        points_in_b.len() < PAR_MIN_POINTS || !par::enabled() || par::force_sequential();
83    let mut hits: Vec<Option<ClosestHit>> = Vec::new();
84    hits.try_reserve_exact(points_in_b.len())
85        .map_err(|_| FieldsError::OutOfMemory)?;
86    hits.resize(points_in_b.len(), None);
87    if sequential {
88        par::set_force_sequential(true);
89    }
90    let found = closest_point_many(mesh_b, bvh_b, points_in_b, &opts, &mut hits);
91    if sequential {
92        par::set_force_sequential(false);
93    }
94    found?;
95    if sequential {
96        for ((&p, hit), v) in points_in_b.iter().zip(&hits).zip(out.iter_mut()) {
97            *v = hit_value(mesh_b, bvh_b, sign, p, hit.as_ref());
98        }
99        return Ok(());
100    }
101    // Every point is computed alone, so chunking changes no value.
102    let oom = || FieldsError::OutOfMemory;
103    let parts = par::map_range(points_in_b.len(), POINT_CHUNK, PAR_MIN_POINTS, oom, |r| {
104        let mut part = Vec::new();
105        part.try_reserve_exact(r.len()).map_err(|_| oom())?;
106        for i in r {
107            part.push(hit_value(
108                mesh_b,
109                bvh_b,
110                sign,
111                points_in_b[i],
112                hits[i].as_ref(),
113            ));
114        }
115        Ok(part)
116    })?;
117    for (chunk, part) in out.chunks_mut(POINT_CHUNK).zip(parts) {
118        chunk.copy_from_slice(&part);
119    }
120    Ok(())
121}
122
123/// The signed distance of one point, as [`signed_distances`] computes it
124/// (the inputs already checked).
125pub(crate) fn point_value(
126    mesh: &Mesh,
127    bvh: &Bvh,
128    sign: &SignData,
129    p: [f32; 3],
130    max_distance: f32,
131) -> Result<f32> {
132    let opts = ClosestOptions { max_distance };
133    let hit = closest_point(mesh, bvh, p, &opts)?;
134    Ok(hit_value(mesh, bvh, sign, p, hit.as_ref()))
135}
136
137/// The signed distance of `p` given its closest point `hit` (`None`: far).
138pub(crate) fn hit_value(
139    mesh: &Mesh,
140    bvh: &Bvh,
141    sign: &SignData,
142    p: [f32; 3],
143    hit: Option<&ClosestHit>,
144) -> f32 {
145    match hit {
146        None => f32::INFINITY,
147        Some(hit) if hit.distance == 0.0 => 0.0,
148        Some(hit) if sign.is_inside(mesh, bvh, p, hit) => -hit.distance,
149        Some(hit) => hit.distance,
150    }
151}
152
153#[cfg(test)]
154mod tests {
155    use fast_mesh_query::{ClosestOptions, Feature, QueryError, closest_point};
156    use fast_mesh_store::{Mesh, par};
157
158    use super::*;
159    use crate::test_util::*;
160    use crate::{FieldsError, SignMode};
161
162    fn distances(mesh: &mut Mesh, points: &[[f32; 3]], max_distance: f32) -> (SignData, Vec<f32>) {
163        let bvh = build_bvh(mesh);
164        let sign = SignData::build(mesh, &bvh, SignMode::Auto).unwrap();
165        let mut out = vec![f32::NAN; points.len()];
166        signed_distances(mesh, &bvh, &sign, points, max_distance, &mut out).unwrap();
167        (sign, out)
168    }
169
170    #[test]
171    fn closed_meshes_match_exact_sign() {
172        let mut rng = Lcg(7);
173        let meshes = [
174            ("cube", cube()),
175            ("icosphere", icosphere()),
176            ("torus", torus()),
177            ("tetrahedron", tetrahedron()),
178            ("star prism", star_prism()),
179        ];
180        for (name, mut mesh) in meshes {
181            let near = diagonal(&mesh) as f32 / 20.0;
182            let points = sample_points(&mesh, &mut rng, 2000, near);
183            let (sign, out) = distances(&mut mesh, &points, f32::INFINITY);
184            assert_eq!(sign.status().mode, SignMode::Pseudonormal, "{name}");
185            assert!(sign.status().closed, "{name}");
186            let checked = check_against_exact(&mesh, &points, &out);
187            assert!(checked > 1900, "{name}: only {checked} signs checked");
188            assert!(out.iter().any(|&v| v < 0.0), "{name}: no point inside");
189            assert!(out.iter().any(|&v| v > 0.0), "{name}: no point outside");
190        }
191    }
192
193    #[test]
194    fn closed_meshes_match_exact_sign_at_vertices_and_edges() {
195        // Points pushed off vertices and edge midpoints, where the closest
196        // feature is a vertex or an edge rather than a face.
197        let mut rng = Lcg(11);
198        let meshes = [
199            ("cube", cube()),
200            ("icosphere", icosphere()),
201            ("torus", torus()),
202            ("tetrahedron", tetrahedron()),
203            ("star prism", star_prism()),
204        ];
205        for (name, mut mesh) in meshes {
206            let d = diagonal(&mesh) as f32;
207            let mut points = Vec::new();
208            for t in 0..mesh.triangle_slots() as u32 {
209                let tri = mesh
210                    .triangle(t)
211                    .unwrap()
212                    .map(|v| mesh.positions()[v as usize]);
213                for k in 0..3 {
214                    let (a, b) = (tri[k], tri[(k + 1) % 3]);
215                    let mid: [f32; 3] = std::array::from_fn(|i| (a[i] + b[i]) / 2.0);
216                    for base in [a, mid] {
217                        let off = rng.point([-0.05 * d; 3], [0.05 * d; 3]);
218                        points.push(std::array::from_fn(|i| base[i] + off[i]));
219                    }
220                }
221            }
222            let (_, out) = distances(&mut mesh, &points, f32::INFINITY);
223            let checked = check_against_exact(&mesh, &points, &out);
224            assert!(
225                checked * 10 > points.len() * 9,
226                "{name}: only {checked} of {} signs checked",
227                points.len()
228            );
229        }
230    }
231
232    #[test]
233    fn open_sheet_sign_past_the_rim_uses_boundary_edges() {
234        // Points just past the rim y = 0 of a terrain sheet, above and below
235        // it: the closest feature is a boundary edge, whose normal is that of
236        // its single face.
237        let n = 16;
238        let mut mesh = terrain(n);
239        let mut rng = Lcg(13);
240        let mut points = Vec::new();
241        let mut below = Vec::new();
242        for i in 0..1000 {
243            let x = 1.0 + (n - 2) as f32 * rng.next();
244            let y = -(0.02 + 0.3 * rng.next());
245            let h = 0.02 + 0.3 * rng.next();
246            let z = terrain_surface(&mesh, n, f64::from(x), 0.0) as f32;
247            let under = i % 2 == 0;
248            points.push([x, y, if under { z - h } else { z + h }]);
249            below.push(under);
250        }
251        let (_, out) = distances(&mut mesh, &points, f32::INFINITY);
252        let bvh = build_bvh(&mesh);
253        let mut on_edge = 0;
254        for ((&p, &value), &under) in points.iter().zip(&out).zip(&below) {
255            let hit = closest_point(&mesh, &bvh, p, &ClosestOptions::default());
256            if matches!(hit.unwrap().unwrap().feature, Feature::Edge(_)) {
257                on_edge += 1;
258            }
259            assert_eq!(value < 0.0, under, "sign at {p:?}: value {value}");
260        }
261        assert!(on_edge > 900, "only {on_edge} closest points on edges");
262    }
263
264    #[test]
265    fn boundary_signal_marks_open_rim() {
266        // Points past the rim report a boundary feature; points over the
267        // interior do not.
268        let n = 16;
269        let mut mesh = terrain(n);
270        let bvh = build_bvh(&mesh);
271        let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
272        assert!(!sign.status().closed);
273        let mut rng = Lcg(13);
274        let mut rim_boundary = 0;
275        for _ in 0..200 {
276            let x = 1.0 + (n - 2) as f32 * rng.next();
277            let y = -(0.02 + 0.3 * rng.next());
278            let z = terrain_surface(&mesh, n, f64::from(x), 0.0) as f32;
279            let p = [x, y, z + 0.1];
280            let hit = closest_point(&mesh, &bvh, p, &ClosestOptions::default())
281                .unwrap()
282                .unwrap();
283            assert_eq!(
284                SignData::is_boundary_hit(&mesh, &hit),
285                Some(true),
286                "at {p:?}"
287            );
288            rim_boundary += 1;
289        }
290        assert_eq!(rim_boundary, 200);
291        let mut interior_clean = 0;
292        for _ in 0..200 {
293            let (x, y) = (4.0 + 8.0 * rng.next(), 4.0 + 8.0 * rng.next());
294            let z = terrain_surface(&mesh, n, f64::from(x), f64::from(y)) as f32;
295            let p = [x, y, z + 0.1];
296            let hit = closest_point(&mesh, &bvh, p, &ClosestOptions::default())
297                .unwrap()
298                .unwrap();
299            if hit.feature == Feature::Face {
300                assert_eq!(SignData::is_boundary_hit(&mesh, &hit), Some(false));
301                interior_clean += 1;
302            }
303        }
304        assert!(interior_clean > 150, "only {interior_clean} face hits");
305    }
306
307    #[test]
308    fn open_sheet_penetration_is_negative() {
309        let n = 16;
310        let mut mesh = terrain(n);
311        let mut rng = Lcg(3);
312        let mut points = Vec::new();
313        let mut below = Vec::new();
314        for i in 0..2000 {
315            let (x, y) = (4.0 + 8.0 * rng.next(), 4.0 + 8.0 * rng.next());
316            let h = 0.01 + 1.5 * rng.next();
317            let z = terrain_surface(&mesh, n, f64::from(x), f64::from(y)) as f32;
318            let under = i % 2 == 0;
319            points.push([x, y, if under { z - h } else { z + h }]);
320            below.push(under);
321        }
322        let (sign, out) = distances(&mut mesh, &points, f32::INFINITY);
323        let status = sign.status();
324        assert_eq!(status.mode, SignMode::Pseudonormal);
325        assert!(!status.closed);
326        assert!(status.consistent_orientation);
327        let tol = 1e-5 * diagonal(&mesh);
328        for ((&p, &value), &under) in points.iter().zip(&out).zip(&below) {
329            let brute = brute_distance(&mesh, p);
330            assert!(
331                (f64::from(value.abs()) - brute).abs() <= tol,
332                "magnitude at {p:?}"
333            );
334            if brute > tol {
335                assert_eq!(value < 0.0, under, "sign at {p:?}: value {value}");
336            }
337        }
338    }
339
340    /// Splits the edge of triangle `t` that starts at corner `k` at its
341    /// midpoint `m`: `(a, b, c)` becomes `(a, m, c)` and `(m, b, c)`, and the
342    /// zero-area triangle `(a, b, m)` fills the gap against the neighbor
343    /// across `a`–`b`, so the mesh stays closed and consistently oriented.
344    /// Returns the zero-area triangle and the midpoint.
345    fn split_with_sliver(
346        positions: &mut Vec<[f32; 3]>,
347        tris: &mut Vec<[u32; 3]>,
348        t: usize,
349        k: usize,
350    ) -> ([u32; 3], [f32; 3]) {
351        let tri = tris[t];
352        let (a, b, c) = (tri[k], tri[(k + 1) % 3], tri[(k + 2) % 3]);
353        let (pa, pb) = (positions[a as usize], positions[b as usize]);
354        let pm = std::array::from_fn(|i| (pa[i] + pb[i]) / 2.0);
355        positions.push(pm);
356        let m = (positions.len() - 1) as u32;
357        tris[t] = [a, m, c];
358        tris.push([m, b, c]);
359        ([a, b, m], pm)
360    }
361
362    #[test]
363    fn degenerate_faces_do_not_break_signs() {
364        let cases = [
365            // A face diagonal, then two cube edges (convex, 90 degrees).
366            (cube_buffers(), vec![(0, 2), (2, 0), (6, 1)]),
367            // Two opposite edges of a tetrahedron, sharp enough that one
368            // face normal alone gives wrong signs next to them.
369            (tetrahedron_buffers(), vec![(0, 0), (3, 1)]),
370        ];
371        for ((mut positions, mut tris), splits) in cases {
372            let mut slivers = Vec::new();
373            let mut mids = Vec::new();
374            for (t, k) in splits {
375                let (s, m) = split_with_sliver(&mut positions, &mut tris, t, k);
376                slivers.push(s);
377                mids.push(m);
378            }
379            let count = slivers.len();
380            // The slivers get the smallest ids, so they win distance ties.
381            slivers.extend(tris);
382            let mut mesh = mesh(&positions, &slivers);
383            let bvh = build_bvh(&mesh);
384            let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
385            let status = sign.status();
386            assert_eq!(status.degenerate_triangles, count);
387            assert!(status.closed);
388            assert!(status.consistent_orientation);
389            assert_eq!(status.mode, SignMode::Pseudonormal);
390
391            let mut rng = Lcg(5);
392            let mut points = sample_points(&mesh, &mut rng, 1000, 0.2);
393            for m in &mids {
394                for _ in 0..300 {
395                    let off = rng.point([-0.4; 3], [0.4; 3]);
396                    points.push(std::array::from_fn(|i| m[i] + off[i]));
397                }
398                // Straight out from the midpoint (the shapes are centred on
399                // the origin): the closest point is the midpoint itself,
400                // shared by a sliver.
401                points.push(m.map(|c| c * 1.5));
402                points.push(m.map(|c| c * 0.5));
403            }
404            let mut out = vec![f32::NAN; points.len()];
405            signed_distances(&mesh, &bvh, &sign, &points, f32::INFINITY, &mut out).unwrap();
406            let checked = check_against_exact(&mesh, &points, &out);
407            assert!(
408                checked * 10 > points.len() * 9,
409                "only {checked} of {} signs checked",
410                points.len()
411            );
412
413            // The fallback for a reported zero-area triangle really ran.
414            let on_sliver = points
415                .iter()
416                .filter(|&&p| {
417                    let hit = closest_point(&mesh, &bvh, p, &ClosestOptions::default());
418                    (hit.unwrap().unwrap().triangle as usize) < count
419                })
420                .count();
421            assert!(on_sliver > 0);
422        }
423    }
424
425    #[test]
426    fn pinch_vertex_signs_match_exact_away_from_pinch() {
427        // Two cubes meeting at one vertex: closed, one non-manifold vertex.
428        // Signs match the exact winding number away from the pinch.
429        let (mut p, mut t) = cube_buffers();
430        let (p2, t2) = cube_buffers();
431        p.extend(p2[1..].iter().map(|q| q.map(|c| c + 2.0)));
432        t.extend(
433            t2.iter()
434                .map(|tri| tri.map(|v| if v == 0 { 7 } else { v + 7 })),
435        );
436        let mut mesh = mesh(&p, &t);
437        let bvh = build_bvh(&mesh);
438        let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
439        let status = sign.status();
440        assert!(status.closed && status.consistent_orientation);
441        assert_eq!(status.non_manifold_vertices, 1);
442        assert_eq!(status.mode, SignMode::Pseudonormal);
443        let pinch = [1.0f32, 1.0, 1.0];
444        let mut rng = Lcg(41);
445        let points = sample_points(&mesh, &mut rng, 2000, 0.2);
446        let mut out = vec![f32::NAN; points.len()];
447        signed_distances(&mesh, &bvh, &sign, &points, f32::INFINITY, &mut out).unwrap();
448        let tol = 1e-5 * diagonal(&mesh);
449        let mut checked = 0;
450        for (&pt, &value) in points.iter().zip(&out) {
451            let to_pinch = ((pt[0] - pinch[0]).powi(2)
452                + (pt[1] - pinch[1]).powi(2)
453                + (pt[2] - pinch[2]).powi(2))
454            .sqrt() as f64;
455            if to_pinch < 0.2 || brute_distance(&mesh, pt) <= tol {
456                continue;
457            }
458            let inside = exact_winding(&mesh, pt) >= 0.5;
459            assert_eq!(value < 0.0, inside, "sign at {pt:?}: value {value}");
460            checked += 1;
461        }
462        assert!(checked > 1500, "only {checked} signs checked");
463    }
464
465    #[test]
466    fn t_junction_vertex_signs_are_consistent() {
467        // A T-junction: the right square's mid-edge vertex lies on the left
468        // square's edge without splitting it. The mesh has a slit (open) but
469        // consistent orientation; signs match the side of the sheet away from
470        // the slit.
471        let pos = vec![
472            [0.0, 0.0, 0.0],
473            [1.0, 0.0, 0.0],
474            [1.0, 1.0, 0.0],
475            [0.0, 1.0, 0.0],
476            [2.0, 0.0, 0.0],
477            [2.0, 1.0, 0.0],
478            [1.0, 0.5, 0.0],
479            [2.0, 0.5, 0.0],
480        ];
481        let tris = vec![
482            [0, 1, 2],
483            [0, 2, 3],
484            [1, 4, 7],
485            [1, 7, 6],
486            [6, 7, 5],
487            [6, 5, 2],
488        ];
489        let mut mesh = mesh(&pos, &tris);
490        let bvh = build_bvh(&mesh);
491        let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
492        let status = sign.status();
493        assert!(!status.closed);
494        assert!(status.consistent_orientation);
495        assert_eq!(status.mode, SignMode::Pseudonormal);
496        let mut rng = Lcg(43);
497        for _ in 0..500 {
498            let x = rng.next() * 2.0;
499            let y = rng.next();
500            if (x - 1.0).abs() < 0.1 {
501                continue;
502            }
503            let h = 0.05 + rng.next() * 0.2;
504            for under in [false, true] {
505                let p = [x, y, if under { -h } else { h }];
506                let mut out = [f32::NAN; 1];
507                signed_distances(&mesh, &bvh, &sign, &[p], f32::INFINITY, &mut out).unwrap();
508                assert_eq!(out[0] < 0.0, under, "sign at {p:?}: {}", out[0]);
509            }
510        }
511    }
512
513    #[test]
514    fn open_mesh_with_flipped_patch_reports_inconsistent() {
515        // A terrain with one interior triangle flipped: open, inconsistent
516        // orientation, pseudonormal mode, and wrong signs on that patch.
517        let base = terrain(8);
518        let n = 8u32;
519        let flip = (4 * n + 4) * 2;
520        let mut tris: Vec<[u32; 3]> = (0..base.triangle_slots() as u32)
521            .map(|t| base.triangle(t).unwrap())
522            .collect();
523        tris[flip as usize].swap(1, 2);
524        let positions: Vec<[f32; 3]> = base.positions().to_vec();
525        let mut mesh = mesh(&positions, &tris);
526        let bvh = build_bvh(&mesh);
527        let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
528        let status = sign.status();
529        assert!(!status.closed);
530        assert!(!status.consistent_orientation);
531        assert_eq!(status.mode, SignMode::Pseudonormal);
532        let tri = tris[flip as usize].map(|v| positions[v as usize]);
533        let centroid: [f32; 3] = std::array::from_fn(|k| (tri[0][k] + tri[1][k] + tri[2][k]) / 3.0);
534        for h in [0.05, 0.15] {
535            for under in [false, true] {
536                let p = [
537                    centroid[0],
538                    centroid[1],
539                    centroid[2] + if under { -h } else { h },
540                ];
541                let mut out = [f32::NAN; 1];
542                signed_distances(&mesh, &bvh, &sign, &[p], f32::INFINITY, &mut out).unwrap();
543                assert_ne!(out[0] < 0.0, under, "flipped sign at {p:?}: {}", out[0]);
544            }
545        }
546    }
547
548    #[test]
549    fn pseudonormal_signs_on_open_cube() {
550        // The open-cube fixture (top face removed): center reads inside,
551        // far outside reads outside, under pseudonormals.
552        let (p, mut t) = cube_buffers();
553        t.drain(2..4);
554        let mut mesh = mesh(&p, &t);
555        let bvh = build_bvh(&mesh);
556        let sign = SignData::build(&mut mesh, &bvh, SignMode::Pseudonormal).unwrap();
557        let status = sign.status();
558        assert!(!status.closed);
559        assert!(status.consistent_orientation);
560        assert_eq!(status.mode, SignMode::Pseudonormal);
561        let points = [
562            ([0.0, 0.0, 0.0], true),
563            ([0.5, -0.3, 0.2], true),
564            ([3.0, 0.0, 0.0], false),
565            ([0.0, 3.0, 0.0], false),
566            ([0.0, 0.0, -3.0], false),
567        ];
568        let pts: Vec<[f32; 3]> = points.iter().map(|(p, _)| *p).collect();
569        let mut out = vec![f32::NAN; pts.len()];
570        signed_distances(&mesh, &bvh, &sign, &pts, f32::INFINITY, &mut out).unwrap();
571        for ((p, inside), value) in points.iter().zip(out) {
572            assert_eq!(value < 0.0, *inside, "at {p:?}: {value}");
573        }
574    }
575
576    #[test]
577    fn far_points_are_infinite() {
578        let mut mesh = cube();
579        let bvh = build_bvh(&mesh);
580        let sign = SignData::build(&mut mesh, &bvh, SignMode::Pseudonormal).unwrap();
581        let points = [
582            [3.0, 0.0, 0.0],   // 2 outside: far
583            [1.25, 0.0, 0.0],  // 0.25 outside
584            [0.0, 0.0, 0.0],   // 1 inside: deeper than max_distance, far
585            [0.0, 0.0, -0.75], // 0.25 inside
586            [1.0, 0.5, 0.25],  // on the surface
587            [1.5, 0.0, 0.0],   // exactly at max_distance: kept
588        ];
589        let mut out = [f32::NAN; 6];
590        signed_distances(&mesh, &bvh, &sign, &points, 0.5, &mut out).unwrap();
591        assert_eq!(out[0], f32::INFINITY);
592        assert_eq!(out[1], 0.25);
593        assert_eq!(out[2], f32::INFINITY);
594        assert_eq!(out[3], -0.25);
595        assert_eq!(out[4].to_bits(), 0.0f32.to_bits());
596        assert_eq!(out[5], 0.5);
597
598        // An infinite max_distance computes every point.
599        signed_distances(&mesh, &bvh, &sign, &points, f32::INFINITY, &mut out).unwrap();
600        assert_eq!(out[0], 2.0);
601        assert_eq!(out[2], -1.0);
602
603        // Invalid input is rejected before anything is written.
604        let mut out = [7.0f32; 6];
605        let mut bad = points;
606        bad[3][1] = f32::NAN;
607        let invalid = |r: Result<()>| matches!(r, Err(FieldsError::InvalidInput(_)));
608        assert!(invalid(signed_distances(
609            &mesh, &bvh, &sign, &bad, 0.5, &mut out
610        )));
611        bad[3][1] = f32::INFINITY;
612        assert!(invalid(signed_distances(
613            &mesh, &bvh, &sign, &bad, 0.5, &mut out
614        )));
615        assert!(invalid(signed_distances(
616            &mesh,
617            &bvh,
618            &sign,
619            &points,
620            f32::NAN,
621            &mut out
622        )));
623        assert!(invalid(signed_distances(
624            &mesh, &bvh, &sign, &points, -1.0, &mut out
625        )));
626        assert!(invalid(signed_distances(
627            &mesh,
628            &bvh,
629            &sign,
630            &points[..5],
631            0.5,
632            &mut out
633        )));
634        assert_eq!(out, [7.0; 6]);
635
636        // A mesh with no triangles is far from everything.
637        let mut empty = Mesh::new();
638        let ebvh = build_bvh(&empty);
639        let esign = SignData::build(&mut empty, &ebvh, SignMode::Auto).unwrap();
640        let mut out = [0.0f32; 2];
641        signed_distances(&empty, &ebvh, &esign, &points[..2], f32::INFINITY, &mut out).unwrap();
642        assert_eq!(out, [f32::INFINITY; 2]);
643    }
644
645    /// Thread counts the parallel paths are compared under.
646    const POOLS: [usize; 3] = [1, 2, 8];
647
648    fn bits(values: &[f32]) -> Vec<u32> {
649        values.iter().map(|x| x.to_bits()).collect()
650    }
651
652    #[test]
653    fn signed_distances_bit_identical() {
654        // The forced-sequential run is the code a feature-off build runs.
655        let mut mesh = terrain(317);
656        assert!(mesh.triangle_count() >= 200_000);
657        let bvh = build_bvh(&mesh);
658        // Still dozens of chunks each (`POINT_CHUNK`), at a fraction of the
659        // debug-build time of larger batches.
660        let points = sample_points(&mesh, &mut Lcg(31), 30_000, 2.0);
661        for (mode, count) in [(SignMode::Pseudonormal, 30_000), (SignMode::Winding, 8_000)] {
662            let sign = SignData::build(&mut mesh, &bvh, mode).unwrap();
663            let points = &points[..count];
664            for max_distance in [f32::INFINITY, 0.5] {
665                let run = || {
666                    let mut out = vec![f32::NAN; points.len()];
667                    signed_distances(&mesh, &bvh, &sign, points, max_distance, &mut out).unwrap();
668                    bits(&out)
669                };
670                par::set_force_sequential(true);
671                let reference = run();
672                par::set_force_sequential(false);
673                if !par::enabled() {
674                    // Every run below would be this one again.
675                    continue;
676                }
677                for k in POOLS {
678                    let got = par::with_threads(k, run);
679                    assert!(got == reference, "{mode:?}, {max_distance}: {k} threads");
680                }
681                let before = par::parallel_entries();
682                assert!(run() == reference, "{mode:?}, {max_distance}: global pool");
683                assert!(par::parallel_entries() > before, "ran sequentially");
684            }
685        }
686    }
687
688    #[test]
689    fn small_inputs_stay_sequential() {
690        let mut mesh = icosphere();
691        let before = par::parallel_entries();
692        let points = sample_points(&mesh, &mut Lcg(2), 400, 0.2);
693        for mode in [SignMode::Pseudonormal, SignMode::Winding] {
694            let bvh = build_bvh(&mesh);
695            let sign = SignData::build(&mut mesh, &bvh, mode).unwrap();
696            let mut out = vec![f32::NAN; points.len()];
697            signed_distances(&mesh, &bvh, &sign, &points, f32::INFINITY, &mut out).unwrap();
698        }
699        assert_eq!(par::parallel_entries(), before);
700    }
701
702    #[test]
703    fn stale_or_foreign_inputs_are_rejected() {
704        let mut mesh = cube();
705        let mut bvh = build_bvh(&mesh);
706        let sign = SignData::build(&mut mesh, &bvh, SignMode::Auto).unwrap();
707        let points = [[2.0, 0.0, 0.0]];
708        let mut out = [7.0f32];
709
710        let other = cube();
711        let other_bvh = build_bvh(&other);
712        assert_eq!(
713            signed_distances(&mesh, &other_bvh, &sign, &points, 1.0, &mut out),
714            Err(FieldsError::Query(QueryError::MeshMismatch))
715        );
716        assert_eq!(
717            signed_distances(&other, &other_bvh, &sign, &points, 1.0, &mut out),
718            Err(FieldsError::InvalidInput(
719                "sign data is not current for this mesh and BVH"
720            ))
721        );
722
723        let mut tx = mesh.begin_edit();
724        tx.set_positions(&[0], &[[-1.0, -1.0, -1.5]]).unwrap();
725        let patch = tx.commit().unwrap();
726        assert_eq!(
727            signed_distances(&mesh, &bvh, &sign, &points, 1.0, &mut out),
728            Err(FieldsError::Query(QueryError::StaleBvh))
729        );
730        assert_eq!(
731            SignData::build(&mut mesh, &bvh, SignMode::Auto).err(),
732            Some(FieldsError::Query(QueryError::StaleBvh))
733        );
734        bvh.update(&mut mesh, patch.changes()).unwrap();
735        assert_eq!(
736            signed_distances(&mesh, &bvh, &sign, &points, 1.0, &mut out),
737            Err(FieldsError::InvalidInput(
738                "sign data is not current for this mesh and BVH"
739            ))
740        );
741        assert_eq!(out, [7.0]);
742    }
743}
fast_mesh_fields/distance.rs

fast_mesh_fields/
distance.rs
