ogeom_algo/distance/
mesh.rs

1//! Distances to a triangle mesh, through its hierarchy.
2
3use ogeom_core::{FastMap, OgeomResult, parallel};
4use ogeom_math::{Point, Vector};
5use ogeom_mesh::{Bvh, Hit};
6
7use super::{Closest, DistanceOptions, MeshFeature, signed, unit, unit_or_zero};
8
9/// The signed distance from each point to the mesh `index` was built over,
10/// with the closest point, the feature it lies on and the normal there; in
11/// the points' order, `None` where nothing is within
12/// [`DistanceOptions::max`].
13///
14/// The closest point is [`Bvh::nearest`]'s, and with a facing limit
15/// [`Bvh::nearest_where`]'s over the triangles facing within the limit of
16/// the point's normal: exactly what brute force over the same triangles
17/// finds, ties to the lowest triangle. A point that is not finite, or with
18/// a facing limit a normal that has no direction, answers `None`. See the
19/// [module](super) for how the sign is read.
20///
21/// # Errors
22///
23/// As [`DistanceOptions::check`];
24/// [`OgeomError::Dimension`](ogeom_core::OgeomError::Dimension) if there are
25/// normals but not one per point, and
26/// [`OgeomError::Domain`](ogeom_core::OgeomError::Domain) if there is a
27/// facing limit and no normals.
28pub fn to_mesh(
29    index: &Bvh,
30    points: &[Point],
31    normals: Option<&[Vector]>,
32    options: &DistanceOptions,
33) -> OgeomResult<Vec<Option<Closest<MeshFeature>>>> {
34    to_mesh_with(index, &Pseudonormals::of(index), points, normals, options)
35}
36
37/// As [`to_mesh`], with the mesh's normals already gathered from `index`,
38/// for a caller that measures to the same mesh many times.
39pub(crate) fn to_mesh_with(
40    index: &Bvh,
41    pseudo: &Pseudonormals,
42    points: &[Point],
43    normals: Option<&[Vector]>,
44    options: &DistanceOptions,
45) -> OgeomResult<Vec<Option<Closest<MeshFeature>>>> {
46    let facing = options.checked(points.len(), normals)?;
47    let hits = match (facing, normals) {
48        (Some(cos), Some(normals)) => {
49            let wanted: Vec<Option<Vector>> = normals.iter().map(|&n| unit(n)).collect();
50            index.nearest_many_where(points, options.max, |i, t| {
51                wanted[i].is_some_and(|m| pseudo.faces(t, m, cos))
52            })
53        }
54        _ => index.nearest_many(points, options.max),
55    };
56    Ok(parallel::map_ordered(&hits, |i, hit| {
57        hit.map(|hit| pseudo.closest(points[i], &hit, options.signed))
58    }))
59}
60
61/// A mesh's normals at its faces, edges and vertices, the corners welded
62/// where their positions are equal.
63pub(crate) struct Pseudonormals {
64    /// Each triangle's unit normal, zero for a degenerate one.
65    face: Vec<Vector>,
66    /// Each triangle's edges' pseudonormals, edge `k` from corner `k`: the
67    /// sum of the unit normals of every triangle on the welded edge.
68    edge: Vec<[Vector; 3]>,
69    /// Each triangle's corners' welded vertices.
70    corner: Vec<[u32; 3]>,
71    /// Each welded vertex's pseudonormal: its triangles' unit normals
72    /// weighted by their angles at it.
73    vertex: Vec<Vector>,
74}
75
76/// A triangle's corners, or its edges by the corner each starts at.
77const CORNERS: [u8; 3] = [0, 1, 2];
78
79/// A position's key: its coordinates' bits, with negative zero as zero so
80/// that equal positions weld.
81fn key(p: Point) -> [u64; 3] {
82    [p.x, p.y, p.z].map(|c| (c + 0.0).to_bits())
83}
84
85impl Pseudonormals {
86    pub(crate) fn of(index: &Bvh) -> Self {
87        let count = index.triangle_count();
88        let triangles: Vec<[Point; 3]> = (0..count)
89            .filter_map(|t| index.triangle(u32::try_from(t).ok()?))
90            .collect();
91        let mut welded: FastMap<[u64; 3], u32> = FastMap::default();
92        let mut vertex: Vec<Vector> = Vec::new();
93        let mut face = Vec::with_capacity(count);
94        let mut corner = Vec::with_capacity(count);
95        for t in &triangles {
96            let ids = t.map(|p| {
97                let next = u32::try_from(vertex.len()).unwrap_or(u32::MAX);
98                let id = *welded.entry(key(p)).or_insert(next);
99                if id == next {
100                    vertex.push(Vector::ZERO);
101                }
102                id
103            });
104            let normal = unit_or_zero((t[1] - t[0]).cross(t[2] - t[0]));
105            for k in 0..3 {
106                let (a, b, c) = (t[k], t[(k + 1) % 3], t[(k + 2) % 3]);
107                let (u, v) = (b - a, c - a);
108                let angle = u.cross(v).magnitude().atan2(u.dot(v));
109                vertex[ids[k] as usize] += normal * angle;
110            }
111            face.push(normal);
112            corner.push(ids);
113        }
114        let mut edges: FastMap<(u32, u32), Vector> = FastMap::default();
115        let edge_of = |ids: &[u32; 3], k: usize| {
116            let (a, b) = (ids[k], ids[(k + 1) % 3]);
117            (a.min(b), a.max(b))
118        };
119        for (ids, &normal) in corner.iter().zip(&face) {
120            for k in 0..3 {
121                let sum = edges.entry(edge_of(ids, k)).or_insert(Vector::ZERO);
122                *sum += normal;
123            }
124        }
125        let edge = corner
126            .iter()
127            .map(|ids| [0, 1, 2].map(|k| edges[&edge_of(ids, k)]))
128            .collect();
129        Self {
130            face,
131            edge,
132            corner,
133            vertex,
134        }
135    }
136
137    /// Whether triangle `t` faces within the limit whose cosine is `cos` of
138    /// the unit normal `m`. A degenerate triangle faces nowhere.
139    fn faces(&self, t: u32, m: Vector, cos: f64) -> bool {
140        let n = self.face[t as usize];
141        n != Vector::ZERO && n.dot(m) >= cos
142    }
143
144    /// The hit as the feature it lies on, by its barycentric weights: a
145    /// corner where two weights are zero, the edge between the other two
146    /// where one is, the face where none is.
147    fn closest(&self, p: Point, hit: &Hit, sign: bool) -> Closest<MeshFeature> {
148        let t = hit.triangle;
149        let at = t as usize;
150        let zero = hit.barycentric.map(|w| w == 0.0);
151        let (feature, normal) = match zero {
152            [false, false, false] => (MeshFeature::Face(t), self.face[at]),
153            [true, true, _] | [true, _, true] | [_, true, true] => {
154                let k = zero.iter().position(|z| !z).unwrap_or(0);
155                (
156                    MeshFeature::Vertex(t, CORNERS[k]),
157                    self.vertex[self.corner[at][k] as usize],
158                )
159            }
160            _ => {
161                let j = zero.iter().position(|&z| z).unwrap_or(0);
162                let k = (j + 1) % 3;
163                (MeshFeature::Edge(t, CORNERS[k]), self.edge[at][k])
164            }
165        };
166        Closest {
167            distance: signed(hit.distance, p, hit.point, normal, sign),
168            point: hit.point,
169            feature,
170            normal: unit_or_zero(normal),
171        }
172    }
173}
ogeom_algo/distance/mesh.rs

ogeom_algo/distance/
mesh.rs
