Skip to main content

fast_mesh_ops/
measure.rs

1//! Scalar measurements of a mesh: surface area and enclosed volume.
2//!
3//! [`measure`] reports the total triangle area of `mesh` (always defined)
4//! and, when the mesh is watertight, its enclosed volume. Watertight here
5//! reuses [`analyze`](crate::analyze): no boundary loop or open chain and
6//! no non-manifold edge. An open or non-manifold mesh still gets its area;
7//! only `volume` is `None`, with [`VolumeReason`] saying why.
8//!
9//! Volume follows the tolerant convention of this crate: it is the absolute
10//! value of the signed divergence-theorem sum, `|sum dot(p, cross(q,
11//! r)) / 6|`, so the sign never depends on whether the mesh is wound
12//! inward or outward. A mesh with inconsistent winding
13//! (`orientation_conflicts > 0`, for example from flipped triangles) still
14//! reports a volume, but the opposing triangles partly cancel, so the
15//! number is unreliable; the conflict count flags that. This favours
16//! returning a number with a warning over failing on slightly broken
17//! inputs; callers that need strictness should require
18//! `orientation_conflicts == 0` themselves. Duplicate triangles (the same
19//! vertex set twice) count twice, like any other triangle.
20//!
21//! Positions are read as the stored `f32` local coordinates and all
22//! geometry is computed in `f64`, like everywhere in this crate. Area and
23//! volume are translation invariant, so the origin only matters to
24//! [`measure_placed`], which maps each mesh point (origin plus local
25//! position) through a caller closure, typically a placement, before
26//! accumulating: the report is then in world units and exact under any
27//! positive scale, including non-uniform ones.
28
29use fast_mesh_store::{Mesh, VertexId};
30
31use crate::analyze::{Analysis, AnalyzeOptions, analyze};
32use crate::{Result, allocation_allowed};
33
34/// Why [`MeasureReport::volume`] is `None`.
35#[derive(Debug, Clone, Copy, PartialEq, Eq)]
36pub enum VolumeReason {
37    /// The mesh has a boundary: at least one loop or open chain.
38    Open,
39    /// The mesh has at least one non-manifold edge (reported first when
40    /// the mesh is both open and non-manifold).
41    NonManifold,
42}
43
44/// Scalar measurements of a mesh, see [`measure`].
45#[derive(Debug, Clone, PartialEq)]
46pub struct MeasureReport {
47    /// The total triangle area.
48    pub area: f64,
49    /// The signed divergence-theorem sum, `sum dot(p, cross(q, r)) / 6`.
50    /// Its sign follows the winding (negative for inward-wound closed
51    /// meshes); meaningless when the mesh is not watertight, since the sum
52    /// then depends on where the missing caps would be.
53    pub signed_volume: f64,
54    /// The enclosed volume, `|signed_volume|`, when the mesh is watertight
55    /// ([`closed`](MeasureReport::closed)); `None` otherwise, with
56    /// [`volume_reason`](MeasureReport::volume_reason) saying why.
57    pub volume: Option<f64>,
58    /// No boundary loop or open chain and no non-manifold edge.
59    pub closed: bool,
60    /// Why `volume` is `None`; `None` exactly when `closed`.
61    pub volume_reason: Option<VolumeReason>,
62    /// The boundary loops and open chains, as [`analyze`] reports them.
63    pub boundary_loops: usize,
64    /// How many edges have three or more triangles.
65    pub non_manifold_edges: usize,
66    /// How many interior edges are used in the same direction by both
67    /// their triangles (flipped winding nearby makes `volume` unreliable).
68    pub orientation_conflicts: usize,
69}
70
71/// The watertight verdict and defect counts behind a [`MeasureReport`].
72struct Verdict {
73    closed: bool,
74    volume_reason: Option<VolumeReason>,
75    boundary_loops: usize,
76    non_manifold_edges: usize,
77    orientation_conflicts: usize,
78}
79
80/// The verdict from an [`Analysis`]: watertight exactly when there is no
81/// boundary and no non-manifold edge.
82fn verdict(a: &Analysis) -> Verdict {
83    let closed = a.loops.is_empty() && a.non_manifold_edges == 0;
84    let volume_reason = if closed {
85        None
86    } else if a.non_manifold_edges > 0 {
87        Some(VolumeReason::NonManifold)
88    } else {
89        Some(VolumeReason::Open)
90    };
91    Verdict {
92        closed,
93        volume_reason,
94        boundary_loops: a.loops.len(),
95        non_manifold_edges: a.non_manifold_edges,
96        orientation_conflicts: a.orientation_conflicts,
97    }
98}
99
100/// The area and signed volume of the alive triangles, each mapped through
101/// `place` first. Triangles with a repeated vertex id contribute nothing.
102fn accumulate(
103    mesh: &Mesh,
104    origin: [f64; 3],
105    mut place: impl FnMut([f64; 3]) -> [f64; 3],
106) -> Result<(f64, f64)> {
107    allocation_allowed()?;
108    let positions = mesh.positions();
109    let mut point = |v: VertexId| {
110        let p = positions[v as usize];
111        place([
112            origin[0] + f64::from(p[0]),
113            origin[1] + f64::from(p[1]),
114            origin[2] + f64::from(p[2]),
115        ])
116    };
117    let (mut area, mut signed) = (0.0, 0.0);
118    for t in 0..mesh.triangle_slots() as u32 {
119        let Some([a, b, c]) = mesh.triangle(t) else {
120            continue;
121        };
122        if a == b || b == c || a == c {
123            continue;
124        }
125        let (p, q, r) = (point(a), point(b), point(c));
126        let u = [q[0] - p[0], q[1] - p[1], q[2] - p[2]];
127        let v = [r[0] - p[0], r[1] - p[1], r[2] - p[2]];
128        let n = [
129            u[1] * v[2] - u[2] * v[1],
130            u[2] * v[0] - u[0] * v[2],
131            u[0] * v[1] - u[1] * v[0],
132        ];
133        area += 0.5 * (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]).sqrt();
134        // Divergence theorem against the origin: translation invariant
135        // for closed meshes, like the area.
136        let cross = [
137            q[1] * r[2] - q[2] * r[1],
138            q[2] * r[0] - q[0] * r[2],
139            q[0] * r[1] - q[1] * r[0],
140        ];
141        signed += (p[0] * cross[0] + p[1] * cross[1] + p[2] * cross[2]) / 6.0;
142    }
143    Ok((area, signed))
144}
145
146/// Builds a [`MeasureReport`] from accumulated sums and a verdict.
147fn report(area: f64, signed_volume: f64, v: Verdict) -> MeasureReport {
148    MeasureReport {
149        area,
150        signed_volume,
151        volume: v.closed.then(|| signed_volume.abs()),
152        closed: v.closed,
153        volume_reason: v.volume_reason,
154        boundary_loops: v.boundary_loops,
155        non_manifold_edges: v.non_manifold_edges,
156        orientation_conflicts: v.orientation_conflicts,
157    }
158}
159
160/// Measures `mesh` in its local units: total area always, enclosed volume
161/// when watertight. See the [module docs](self) for the tolerant volume
162/// convention.
163///
164/// # Errors
165/// [`OpsError::Store`] with [`StoreError::OutOfMemory`] when a temporary
166/// cannot grow.
167///
168/// [`OpsError::Store`]: crate::OpsError::Store
169/// [`StoreError::OutOfMemory`]: fast_mesh_store::StoreError::OutOfMemory
170pub fn measure(mesh: &Mesh) -> Result<MeasureReport> {
171    let analysis = analyze(mesh, &AnalyzeOptions::default())?;
172    let (area, signed_volume) = accumulate(mesh, [0.0; 3], |p| p)?;
173    Ok(report(area, signed_volume, verdict(&analysis)))
174}
175
176/// Measures `mesh` with each mesh point (the origin plus a local position)
177/// mapped through `place` first, so with a placement the report is in
178/// world units, exact under any positive scale. Flags and the watertight
179/// verdict are placement independent and match [`measure`].
180///
181/// # Errors
182/// As [`measure`].
183pub fn measure_placed(
184    mesh: &Mesh,
185    origin: [f64; 3],
186    place: impl FnMut([f64; 3]) -> [f64; 3],
187) -> Result<MeasureReport> {
188    let analysis = analyze(mesh, &AnalyzeOptions::default())?;
189    let (area, signed_volume) = accumulate(mesh, origin, place)?;
190    Ok(report(area, signed_volume, verdict(&analysis)))
191}
192
193#[cfg(test)]
194mod tests {
195    use super::*;
196    use crate::test_util::*;
197
198    #[test]
199    fn unit_cube_volume_and_area() {
200        let mesh = box_mesh([0.0; 3], [1.0; 3]);
201        let m = measure(&mesh).unwrap();
202        assert!((m.area - 6.0).abs() < 1e-12, "area {}", m.area);
203        assert!((m.volume.unwrap() - 1.0).abs() < 1e-12, "{m:?}");
204        assert!(m.signed_volume.abs() - 1.0 < 1e-12);
205        assert!(m.closed && m.volume_reason.is_none());
206        assert_eq!(
207            (
208                m.boundary_loops,
209                m.non_manifold_edges,
210                m.orientation_conflicts
211            ),
212            (0, 0, 0)
213        );
214    }
215
216    #[test]
217    fn tetrahedron_volume_and_area() {
218        // (0,0,0), (1,0,0), (0,1,0), (0,0,1), outward wound: volume 1/6,
219        // area three halves plus the sqrt(2)-sided face.
220        let mesh = mesh_of(
221            &[0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
222            &[[0, 2, 1], [0, 3, 2], [0, 1, 3], [1, 2, 3]],
223        );
224        let m = measure(&mesh).unwrap();
225        assert!(
226            (m.area - (1.5 + 3f64.sqrt() / 2.0)).abs() < 1e-12,
227            "area {}",
228            m.area
229        );
230        assert!((m.volume.unwrap() - 1.0 / 6.0).abs() < 1e-12, "{m:?}");
231        assert!(m.closed && m.volume_reason.is_none());
232        assert_eq!(m.orientation_conflicts, 0);
233    }
234
235    #[test]
236    fn open_grid_has_area_but_no_volume() {
237        let mesh = grid(4);
238        let m = measure(&mesh).unwrap();
239        assert!((m.area - 16.0).abs() < 1e-9, "area {}", m.area);
240        assert_eq!(m.volume, None);
241        assert_eq!(m.volume_reason, Some(VolumeReason::Open));
242        assert!(!m.closed);
243        assert_eq!(m.boundary_loops, 1);
244    }
245
246    #[test]
247    fn non_manifold_fin_reports_that_reason() {
248        let m = measure(&fin()).unwrap();
249        assert_eq!(m.volume, None);
250        assert_eq!(m.volume_reason, Some(VolumeReason::NonManifold));
251        assert!(!m.closed);
252        assert!(m.non_manifold_edges > 0);
253    }
254
255    #[test]
256    fn flipped_triangles_keep_a_flagged_volume() {
257        // Tolerant: a volume is still reported, with the conflict count
258        // warning that opposing triangles cancelled.
259        let sphere = uv_sphere(8, 12, 1.0);
260        let whole = measure(&sphere).unwrap().volume.unwrap();
261        let flipped = flip(&sphere, &[0, 1, 2]);
262        let m = measure(&flipped).unwrap();
263        assert!(m.closed && m.volume.is_some());
264        assert!(m.orientation_conflicts > 0);
265        assert!(
266            (m.volume.unwrap() - whole).abs() > 1e-6,
267            "cancellation should move the volume"
268        );
269    }
270
271    #[test]
272    fn placed_cube_scales_area_and_volume() {
273        // A uniform scale of 2 about the mesh point: area x4, volume x8.
274        let mesh = box_mesh([0.0; 3], [1.0; 3]);
275        let m = measure_placed(&mesh, [0.0; 3], |p| [2.0 * p[0], 2.0 * p[1], 2.0 * p[2]]).unwrap();
276        assert!((m.area - 24.0).abs() < 1e-9, "area {}", m.area);
277        assert!((m.volume.unwrap() - 8.0).abs() < 1e-9, "{m:?}");
278        // Non-uniform scales are exact too: 2 x 1 x 1 doubles both.
279        let m = measure_placed(&mesh, [0.0; 3], |p| [2.0 * p[0], p[1], p[2]]).unwrap();
280        assert!((m.area - 10.0).abs() < 1e-9, "area {}", m.area);
281        assert!((m.volume.unwrap() - 2.0).abs() < 1e-9, "{m:?}");
282        // A translation alone changes nothing.
283        let m = measure_placed(&mesh, [0.0; 3], |p| [p[0] + 5.0, p[1] - 3.0, p[2] + 1.0]).unwrap();
284        assert!((m.area - 6.0).abs() < 1e-12);
285        assert!((m.volume.unwrap() - 1.0).abs() < 1e-12);
286    }
287
288    #[test]
289    fn placed_matches_local_flags() {
290        // Flags and the verdict do not depend on the mapping.
291        let mesh = damaged_sphere();
292        let a = measure(&mesh).unwrap();
293        let b = measure_placed(&mesh, [1.0, 2.0, 3.0], |p| [p[0] + 4.0, p[1], p[2]]).unwrap();
294        assert_eq!(
295            (
296                a.closed,
297                a.volume_reason,
298                a.boundary_loops,
299                a.non_manifold_edges,
300                a.orientation_conflicts
301            ),
302            (
303                b.closed,
304                b.volume_reason,
305                b.boundary_loops,
306                b.non_manifold_edges,
307                b.orientation_conflicts
308            )
309        );
310        assert!((a.area - b.area).abs() < 1e-9);
311    }
312}