axiolid_measure/

mesh.rs

//! Deterministic raw triangle-mesh measures.
use axiolid_core::{Point3, Tolerance, Vec3};
use axiolid_mesh::{audit_mesh, audit_mesh_cycle, MeshCycleHealth, MeshHealth, TriangleMeshView};
use core::fmt;

#[derive(Debug, Clone, Copy, PartialEq)]
pub struct SurfaceProperties {
    pub area: f64,
    pub centroid: Point3,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct VolumeProperties {
    pub signed_volume: f64,
    pub centroid: Point3,
}
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum MeshMeasureError {
    MeshNotSurfaceUsable(MeshHealth),
    MeshNotVolumeUsable(MeshHealth),
    ZeroSurfaceArea,
    ZeroSignedVolume,
}
impl fmt::Display for MeshMeasureError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        f.write_str("mesh is not suitable for requested raw measure")
    }
}
impl std::error::Error for MeshMeasureError {}

/// Area of one triangle given its corners about the mesh's local origin.
///
/// Shared by [`surface_properties`] and face selection so both sum the same
/// per-triangle values bit for bit.
pub(crate) fn triangle_area([a, b, c]: [Vec3; 3]) -> f64 {
    (b - a).cross(c - a).length() * 0.5
}

pub fn surface_properties<M: TriangleMeshView + ?Sized>(
    mesh: &M,
    tolerance: Tolerance,
) -> Result<SurfaceProperties, MeshMeasureError> {
    let health = audit_mesh(mesh, tolerance);
    if !health.is_surface_usable() || health.degenerate_triangles != 0 {
        return Err(MeshMeasureError::MeshNotSurfaceUsable(health));
    }
    // Area uses edge differences, which are already local and well
    // conditioned. The centroid is not: it accumulates absolute positions,
    // so at large coordinates it loses the same precision `volume_properties`
    // did. Re-base it the same way.
    let base = mesh.position(0);
    let (mut area, mut weighted) = (0.0f64, Vec3::ZERO);
    for i in 0..mesh.triangle_count() {
        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize) - base);
        let weight = triangle_area([a, b, c]);
        area += weight;
        weighted += (a + b + c) * (weight / 3.0);
    }
    if !area.is_finite() || area == 0.0 {
        return Err(MeshMeasureError::ZeroSurfaceArea);
    }
    Ok(SurfaceProperties {
        area,
        // Shift the local-frame centroid back to world coordinates.
        centroid: base + (weighted / area),
    })
}

pub fn volume_properties<M: TriangleMeshView + ?Sized>(
    mesh: &M,
    tolerance: Tolerance,
) -> Result<VolumeProperties, MeshMeasureError> {
    let health = audit_mesh(mesh, tolerance);
    if !health.is_closed_two_manifold() {
        return Err(MeshMeasureError::MeshNotVolumeUsable(health));
    }
    divergence_volume(mesh).ok_or(MeshMeasureError::ZeroSignedVolume)
}

/// Divergence-theorem sum shared by [`volume_properties`] and
/// [`cycle_volume_properties`]; `None` when the sum is zero or not finite.
fn divergence_volume<M: TriangleMeshView + ?Sized>(mesh: &M) -> Option<VolumeProperties> {
    // Sum about a local origin rather than the world origin.
    //
    // Volume and centroid are translation-invariant, so re-basing the operands
    // changes nothing mathematically -- but it changes the conditioning
    // completely. Evaluating `a . (b x c)` about the world origin scales every
    // term with the CUBE of the distance to it, so a 0.1 m box on a national
    // grid at 1e7 forms terms of order 1e21 whose cancellation has to produce
    // 1e-3. That measured 25% volume error and a centroid 8160 km off the
    // solid: a rule check reading either got a confidently wrong number.
    //
    // The first vertex is the local origin. Any point near the geometry works;
    // this one is free and always on the mesh.
    let base = mesh.position(0);
    let (mut volume, mut weighted) = (0.0f64, Vec3::ZERO);
    for i in 0..mesh.triangle_count() {
        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize) - base);
        let tetra = a.dot(b.cross(c)) / 6.0;
        volume += tetra;
        weighted += (a + b + c) * (tetra / 4.0);
    }
    if !volume.is_finite() || volume == 0.0 {
        return None;
    }
    Some(VolumeProperties {
        signed_volume: volume,
        // The sum ran in local coordinates, so the centroid comes back
        // relative to `base`. Shift it home.
        centroid: base + (weighted / volume),
    })
}

/// What [`cycle_volume_properties`] refuses.
#[non_exhaustive]
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum CycleVolumeError {
    /// The mesh is not a closed oriented cycle: an index or position is
    /// invalid, no triangle is usable, a triangle is degenerate, or an edge
    /// is used unequally often in each direction. The audit says which.
    NotClosedOrientedCycle(Box<MeshCycleHealth>),
    /// The sum is zero or not finite.
    ZeroSignedVolume,
}

impl fmt::Display for CycleVolumeError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::NotClosedOrientedCycle(_) => f.write_str("mesh is not a closed oriented cycle"),
            Self::ZeroSignedVolume => f.write_str("mesh has zero or non-finite signed volume"),
        }
    }
}

impl std::error::Error for CycleVolumeError {}

/// Signed volume and centroid of a mesh whose triangles form a closed
/// oriented cycle, the weaker input [`volume_properties`] does not accept.
///
/// A cycle is a mesh in which every undirected edge is used equally often in
/// each direction (see `axiolid_mesh::MeshCycleHealth`). It may hold shells
/// that touch along an edge, or an edge shared by four triangles. The
/// divergence sum of such a mesh does not depend on its base point and equals
/// the integral of the integer winding number `w` of the surface; on
/// two-manifold input the result is bit-identical to [`volume_properties`].
///
/// The result is the enclosed volume only when `w` is 0 or 1 everywhere,
/// which is the contract [`volume_properties`] has always had: overlapping or
/// duplicated shells count twice and an inverted shell subtracts. No edge
/// audit can see that; `axiolid_inspect::certify_unit_winding` certifies it
/// (ADR 0099 D6), bounding the volume where `w` is neither 0 nor 1 by the winding excess.
///
/// A T-junction leaves three unbalanced edges and is refused. Run
/// `axiolid_heal::split_t_junctions` (or `RepairAction::SplitTJunctions`)
/// first; it is opt-in because it only splits edges and creates no vertex.
/// A mesh whose vector area vanishes without being closed, such as a box
/// without top and bottom, is refused by the same audit.
pub fn cycle_volume_properties<M: TriangleMeshView + ?Sized>(
    mesh: &M,
    tolerance: Tolerance,
) -> Result<VolumeProperties, CycleVolumeError> {
    let audit = audit_mesh_cycle(mesh, tolerance);
    if !audit.is_closed_oriented_cycle() {
        return Err(CycleVolumeError::NotClosedOrientedCycle(Box::new(audit)));
    }
    divergence_volume(mesh).ok_or(CycleVolumeError::ZeroSignedVolume)
}

/// Second moments of the enclosed solid about the origin.
///
/// The `x`, `y`, `z` components are the integrals of `x^2`, `y^2` and
/// `z^2` over the enclosed volume, i.e. the diagonal of the second-moment
/// tensor for unit density. The classical inertia tensor diagonal is
/// `(Iyy + Izz, Ixx + Izz, Ixx + Iyy)` from these, so callers can derive
/// either convention without the provider guessing which one they meant.
///
/// # Method
///
/// The divergence theorem again, one order higher than the volume sum:
/// for a closed oriented triangulation the integral of `x^2` over the
/// enclosed region reduces to a sum over triangles of terms in the
/// vertices' coordinates. Each triangle contributes
/// `(nx / 60) * sum over the 10 symmetric monomials`, the standard
/// closed-form tetrahedral moment about the origin.
///
/// The sum runs about the mesh's first vertex, as [`volume_properties`]'s
/// does, and is shifted back to the origin by the parallel-axis terms: with
/// `x = x' + o`, `int x_i^2 = int x'_i^2 + 2 o_i int x'_i + o_i^2 V`
/// (#343). About the world origin each term grew with the fifth power of
/// the distance to it, and a box on a national grid lost its central
/// moments entirely to the cancellation.
///
/// Requires the same closed two-manifold input as [`volume_properties`]:
/// an open shell has no enclosed region and the integral is meaningless,
/// so it is refused rather than summed into a plausible number.
pub fn second_moments<M: TriangleMeshView + ?Sized>(
    mesh: &M,
    tolerance: Tolerance,
) -> Result<Vec3, MeshMeasureError> {
    let health = audit_mesh(mesh, tolerance);
    if !health.is_closed_two_manifold() {
        return Err(MeshMeasureError::MeshNotVolumeUsable(health));
    }
    let base = mesh.position(0);
    let (mut moments, mut first, mut volume) = (Vec3::ZERO, Vec3::ZERO, 0.0);
    for i in 0..mesh.triangle_count() {
        let [a, b, c] = mesh.triangle(i).map(|j| mesh.position(j as usize) - base);
        // Signed volume of the tetrahedron on `base`, times 6.
        let six_v = a.dot(b.cross(c));
        volume += six_v / 6.0;
        first += (a + b + c) * (six_v / 24.0);
        for axis in 0..3 {
            let (pa, pb, pc) = (a[axis], b[axis], c[axis]);
            // Integral of t^2 over the tetrahedron, in barycentric closed
            // form: (a^2 + b^2 + c^2 + ab + ac + bc) / 60 per unit 6V.
            let quadratic = pa * pa + pb * pb + pc * pc + pa * pb + pa * pc + pb * pc;
            moments[axis] += six_v * quadratic / 60.0;
        }
    }
    // Shift the moments about `base` home to the origin.
    let moments = moments + 2.0 * base * first + base * base * volume;
    if !moments.x.is_finite() || !moments.y.is_finite() || !moments.z.is_finite() {
        return Err(MeshMeasureError::ZeroSignedVolume);
    }
    Ok(moments)
}