axiolid_measure/
mesh.rs

1//! Deterministic raw triangle-mesh measures.
2use axiolid_core::{Point3, Tolerance, Vec3};
3use axiolid_mesh::{audit_mesh, MeshHealth, TriangleMeshView};
4use core::fmt;
5
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct SurfaceProperties {
8    pub area: f64,
9    pub centroid: Point3,
10}
11#[derive(Debug, Clone, Copy, PartialEq)]
12pub struct VolumeProperties {
13    pub signed_volume: f64,
14    pub centroid: Point3,
15}
16#[derive(Debug, Clone, PartialEq, Eq)]
17pub enum MeshMeasureError {
18    MeshNotSurfaceUsable(MeshHealth),
19    MeshNotVolumeUsable(MeshHealth),
20    ZeroSurfaceArea,
21    ZeroSignedVolume,
22}
23impl fmt::Display for MeshMeasureError {
24    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
25        f.write_str("mesh is not suitable for requested raw measure")
26    }
27}
28impl std::error::Error for MeshMeasureError {}
29
30pub fn surface_properties<M: TriangleMeshView + ?Sized>(
31    mesh: &M,
32    tolerance: Tolerance,
33) -> Result<SurfaceProperties, MeshMeasureError> {
34    let health = audit_mesh(mesh, tolerance);
35    if !health.is_surface_usable() || health.degenerate_triangles != 0 {
36        return Err(MeshMeasureError::MeshNotSurfaceUsable(health));
37    }
38    // Area uses edge differences, which are already local and well
39    // conditioned. The centroid is not: it accumulates absolute positions,
40    // so at large coordinates it loses the same precision `volume_properties`
41    // did. Re-base it the same way.
42    let base = mesh.position(0);
43    let (mut area, mut weighted) = (0.0f64, Vec3::ZERO);
44    for i in 0..mesh.triangle_count() {
45        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize) - base);
46        let weight = (b - a).cross(c - a).length() * 0.5;
47        area += weight;
48        weighted += (a + b + c) * (weight / 3.0);
49    }
50    if !area.is_finite() || area == 0.0 {
51        return Err(MeshMeasureError::ZeroSurfaceArea);
52    }
53    Ok(SurfaceProperties {
54        area,
55        // Shift the local-frame centroid back to world coordinates.
56        centroid: base + (weighted / area),
57    })
58}
59
60pub fn volume_properties<M: TriangleMeshView + ?Sized>(
61    mesh: &M,
62    tolerance: Tolerance,
63) -> Result<VolumeProperties, MeshMeasureError> {
64    let health = audit_mesh(mesh, tolerance);
65    if !health.is_closed_two_manifold() {
66        return Err(MeshMeasureError::MeshNotVolumeUsable(health));
67    }
68    // Sum about a local origin rather than the world origin.
69    //
70    // Volume and centroid are translation-invariant, so re-basing the operands
71    // changes nothing mathematically -- but it changes the conditioning
72    // completely. Evaluating `a . (b x c)` about the world origin scales every
73    // term with the CUBE of the distance to it, so a 0.1 m box on a national
74    // grid at 1e7 forms terms of order 1e21 whose cancellation has to produce
75    // 1e-3. That measured 25% volume error and a centroid 8160 km off the
76    // solid: a rule check reading either got a confidently wrong number.
77    //
78    // The first vertex is the local origin. Any point near the geometry works;
79    // this one is free and always on the mesh.
80    let base = mesh.position(0);
81    let (mut volume, mut weighted) = (0.0f64, Vec3::ZERO);
82    for i in 0..mesh.triangle_count() {
83        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize) - base);
84        let tetra = a.dot(b.cross(c)) / 6.0;
85        volume += tetra;
86        weighted += (a + b + c) * (tetra / 4.0);
87    }
88    if !volume.is_finite() || volume == 0.0 {
89        return Err(MeshMeasureError::ZeroSignedVolume);
90    }
91    Ok(VolumeProperties {
92        signed_volume: volume,
93        // The sum ran in local coordinates, so the centroid comes back
94        // relative to `base`. Shift it home.
95        centroid: base + (weighted / volume),
96    })
97}
98
99/// Second moments of the enclosed solid about the origin.
100///
101/// The `x`, `y`, `z` components are the integrals of `x^2`, `y^2` and
102/// `z^2` over the enclosed volume, i.e. the diagonal of the second-moment
103/// tensor for unit density. The classical inertia tensor diagonal is
104/// `(Iyy + Izz, Ixx + Izz, Ixx + Iyy)` from these, so callers can derive
105/// either convention without the provider guessing which one they meant.
106///
107/// # Method
108///
109/// The divergence theorem again, one order higher than the volume sum:
110/// for a closed oriented triangulation the integral of `x^2` over the
111/// enclosed region reduces to a sum over triangles of terms in the
112/// vertices' coordinates. Each triangle contributes
113/// `(nx / 60) * sum over the 10 symmetric monomials`, the standard
114/// closed-form tetrahedral moment about the origin.
115///
116/// Requires the same closed two-manifold input as [`volume_properties`]:
117/// an open shell has no enclosed region and the integral is meaningless,
118/// so it is refused rather than summed into a plausible number.
119pub fn second_moments<M: TriangleMeshView + ?Sized>(
120    mesh: &M,
121    tolerance: Tolerance,
122) -> Result<Vec3, MeshMeasureError> {
123    let health = audit_mesh(mesh, tolerance);
124    if !health.is_closed_two_manifold() {
125        return Err(MeshMeasureError::MeshNotVolumeUsable(health));
126    }
127    let mut moments = Vec3::ZERO;
128    for i in 0..mesh.triangle_count() {
129        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize));
130        // Signed volume of the tetrahedron on the origin, times 6.
131        let six_v = a.dot(b.cross(c));
132        for axis in 0..3 {
133            let (pa, pb, pc) = (a[axis], b[axis], c[axis]);
134            // Integral of t^2 over the tetrahedron, in barycentric closed
135            // form: (a^2 + b^2 + c^2 + ab + ac + bc) / 60 per unit 6V.
136            let quadratic = pa * pa + pb * pb + pc * pc + pa * pb + pa * pc + pb * pc;
137            moments[axis] += six_v * quadratic / 60.0;
138        }
139    }
140    if !moments.x.is_finite() || !moments.y.is_finite() || !moments.z.is_finite() {
141        return Err(MeshMeasureError::ZeroSignedVolume);
142    }
143    Ok(moments)
144}
axiolid_measure/mesh.rs

axiolid_measure/
mesh.rs
