phys_geom/

volume.rs

// Copyright (C) 2020-2025 phys-geom authors. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

//! Volume computation and mesh processing utilities.
//!
//! This module provides functionality for computing volumes of geometric shapes and
//! processing triangle meshes. It includes both individual shape volume calculations
//! and mesh-based operations like signed volume computation and bounding box calculation.
//!
//! # Volume Computation
//!
//! The module provides a [`ComputeVolume`] trait that is implemented for all
//! supported geometric shapes. This allows for uniform volume calculation across
//! different shape types.
//!
//! # Mesh Processing
//!
//! Additional utilities are provided for working with triangle meshes:
//! - Signed volume computation for watertight meshes
//! - Axis-aligned bounding box computation
//!
//! # Examples
//!
//! ```rust
//! use phys_geom::shape::{Sphere, Cuboid};
//! use phys_geom::volume::ComputeVolume;
//!
//! // Compute volume of individual shapes
//! let sphere = Sphere::new(1.0);
//! let sphere_volume = sphere.compute_volume(); // 4π/3 ≈ 4.1888
//!
//! let cuboid = Cuboid::new([2.0, 3.0, 1.0]);
//! let cuboid_volume = cuboid.compute_volume(); // 6.0
//! ```
//!
//! # Mathematical Accuracy
//!
//! All volume calculations use mathematically accurate formulas:
//! - Sphere: V = (4/3)πr³
//! - Capsule: V = πr²(h + 4r/3) where h is the cylindrical height
//! - Cuboid: V = l × w × h
//! - Cylinder: V = πr²h
//!
//! The calculations are optimized for performance while maintaining numerical precision.

use crate::math::{Real, PI};
use crate::shape::{Capsule, Cuboid, Cylinder, InfinitePlane, Sphere, Tetrahedron};
use crate::Aabb3;

/// Trait for computing the volume of geometric shapes.
///
/// This trait provides a uniform interface for volume calculation across different
/// geometric shape types. All shapes implement this trait, allowing for generic
/// volume computation algorithms.
///
/// # Examples
///
/// ```rust
/// use phys_geom::shape::{Sphere, Cuboid};
/// use phys_geom::volume::ComputeVolume;
/// use phys_geom::Real;
///
/// fn calculate_total_volume<T: ComputeVolume>(shapes: &[T]) -> Real {
///     shapes.iter().map(|s| s.compute_volume()).sum()
/// }
///
/// // Works with homogeneous collections
/// let spheres = vec![Sphere::new(1.0), Sphere::new(2.0)];
/// let sphere_total = calculate_total_volume(&spheres);
///
/// let cuboids = vec![Cuboid::new([1.0, 1.0, 1.0]), Cuboid::new([2.0, 2.0, 2.0])];
/// let cuboid_total = calculate_total_volume(&cuboids);
/// ```
///
/// # Implementation Notes
///
/// Implementations should use mathematically accurate formulas and be optimized
/// for performance. Volume calculations should be deterministic and return
/// positive values for valid shapes.
pub trait ComputeVolume {
    /// Computes the volume of this shape.
    ///
    /// # Returns
    ///
    /// The volume as a `Real` (f32 or f64 depending on features).
    /// Always returns a non-negative value for valid shapes.
    ///
    /// # Examples
    ///
    /// ```rust
    /// use phys_geom::shape::Sphere;
    /// use phys_geom::volume::ComputeVolume;
    ///
    /// let sphere = Sphere::new(2.0);
    /// let volume = sphere.compute_volume();
    /// assert!(volume > 0.0);
    /// ```
    fn compute_volume(&self) -> Real;
}

impl ComputeVolume for Sphere {
    #[inline]
    fn compute_volume(&self) -> Real {
        const A: Real = 4.0 * PI / 3.0;
        let r = self.radius();
        A * r * r * r
    }
}

impl ComputeVolume for Capsule {
    #[inline]
    fn compute_volume(&self) -> Real {
        const A: Real = 4.0 / 3.0;
        let r = self.radius();
        let h = self.height();
        PI * r * r * (h + A * r)
    }
}

impl ComputeVolume for Cuboid {
    #[inline]
    fn compute_volume(&self) -> Real {
        let half_length = self.half_length();
        half_length[0] * half_length[1] * half_length[2] * 8.0
    }
}

impl ComputeVolume for Cylinder {
    #[inline]
    fn compute_volume(&self) -> Real {
        let r = self.radius();
        let half_v = PI * r * r * self.half_height();
        half_v + half_v
    }
}

impl ComputeVolume for InfinitePlane {
    #[inline]
    fn compute_volume(&self) -> Real {
        0.0
    }
}

/// Computes the signed volume of a triangle mesh using the divergence theorem.
///
/// This function calculates the signed volume of a closed, watertight triangle mesh
/// by summing the signed volumes of tetrahedra formed by each triangle face and the origin.
/// The sign of the volume depends on the winding order of the triangle vertices.
///
/// # Mathematical Background
///
/// The signed volume is computed using the formula:
/// V = Σ(V_tetrahedron) for each triangle, where each tetrahedron is formed
/// by the origin and the triangle vertices.
///
/// # Vertex Ordering
///
/// - **Counter-clockwise** (when viewed from outside): Positive volume
/// - **Clockwise** (when viewed from outside): Negative volume
///
/// The vertex ordering should be consistent across all triangles for accurate results.
///
/// # Use Cases
///
/// - **Mesh validation**: Checking if a mesh is watertight (volume should be consistent)
/// - **Physics simulations**: Computing mass properties of complex objects
/// - **3D modeling**: Calculating volumes of imported meshes
/// - **Computational geometry**: Determining inside/outside relationships
///
/// # Arguments
///
/// * `vertices` - An iterator over triangle vertices, yielding 3 vertices per triangle
///
/// # Returns
///
/// The signed volume as a `Real` value.
///
/// # Panics
///
/// Panics if the number of vertices is not a multiple of 3, indicating malformed
/// triangle data.
///
/// # Examples
///
/// ```rust
/// use phys_geom::volume::triangle_mesh_signed;
///
/// // A simple tetrahedron represented as 4 triangles
/// let vertices = vec![
///     // Triangle 1
///     [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0],
///     // Triangle 2
///     [0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0],
///     // ... more triangles
/// ];
///
/// let volume = triangle_mesh_signed(vertices.into_iter());
/// ```
///
/// # Performance
///
/// The algorithm processes each triangle exactly once, making it O(n) where n is
/// the number of triangles. The computation is efficient and suitable for large meshes.
#[inline]
pub fn triangle_mesh_signed(mut vertices: impl Iterator<Item = [Real; 3]>) -> Real {
    let mut volume = 0.0;
    while let Some(v0) = vertices.next() {
        let v1 = vertices
            .next()
            .expect("The vertices should be a multiple of 3.");
        let v2 = vertices
            .next()
            .expect("The vertices should be a multiple of 3.");
        let dv = Tetrahedron {
            a: [0.0, 0.0, 0.0],
            b: v0,
            c: v1,
            d: v2,
        }
        .signed_volume();
        volume += dv;
    }
    volume
}

/// Computes the axis-aligned bounding box (AABB) of a mesh or point cloud.
///
/// This function calculates the smallest axis-aligned bounding box that contains
/// all the provided vertices. It's useful for:
/// - Quick collision detection broad-phase
/// - Spatial partitioning and culling
/// - Mesh visualization and processing
/// - Bounding volume hierarchy construction
///
/// # Algorithm
///
/// The function iterates through all vertices, tracking the minimum and maximum
/// coordinate values along each axis. The resulting AABB is defined by these
/// extreme points.
///
/// # Arguments
///
/// * `vertices` - An iterator over 3D points [x, y, z]
///
/// # Returns
///
/// An [`Aabb3`] that contains all input vertices.
///
/// # Special Cases
///
/// - **Empty iterator**: Returns an AABB with MAX/MIN values (anti-infinite)
/// - **Single point**: Returns a zero-volume AABB at that point
/// - **Coplanar points**: Returns a flat AABB with zero volume in one or more dimensions
///
/// # Examples
///
/// ```rust
/// use phys_geom::volume::compute_mesh_bound;
///
/// let vertices = vec![
///     [0.0, 0.0, 0.0],
///     [1.0, 0.0, 0.0],
///     [0.5, 1.0, 0.0],
///     [0.5, 0.5, 1.0],
/// ];
///
/// let aabb = compute_mesh_bound(vertices.into_iter());
///
/// // The AABB should contain all points
/// assert_eq!(aabb.min(), [0.0, 0.0, 0.0].into());
/// assert_eq!(aabb.max(), [1.0, 1.0, 1.0].into());
/// ```
///
/// # Performance
///
/// The algorithm runs in O(n) time where n is the number of vertices, with O(1)
/// additional memory usage. It processes each vertex exactly once.
///
/// # Numerical Considerations
///
/// The function uses the floating-point min/max operations directly. For very
/// large or very small coordinate values, consider the precision limitations
/// of the underlying `Real` type.
pub fn compute_mesh_bound(vertices: impl Iterator<Item = [Real; 3]>) -> Aabb3 {
    let mut min = [Real::MAX; 3];
    let mut max = [Real::MIN; 3];
    for point in vertices {
        min = [
            min[0].min(point[0]),
            min[1].min(point[1]),
            min[2].min(point[2]),
        ];
        max = [
            max[0].max(point[0]),
            max[1].max(point[1]),
            max[2].max(point[2]),
        ];
    }
    Aabb3::new_unchecked(min, max)
}

#[cfg(test)]
mod tests {
    use crate::shape::{InfinitePlane, Tetrahedron};
    use crate::volume::ComputeVolume;

    #[test]
    fn test_infinite_plane() {
        let infinite_plane = InfinitePlane::default();
        assert_eq!(infinite_plane.compute_volume(), 0.0);
    }

    #[test]
    fn test_triangle_mesh() {
        let tet = [
            [1.0, 0.0, 0.0],
            [2.0, 0.0, 0.0],
            [1.0, 1.0, 0.0],
            [1.0, 0.0, 1.0],
        ];
        let target_volume = Tetrahedron {
            a: tet[0],
            b: tet[1],
            c: tet[2],
            d: tet[3],
        }
        .signed_volume();
        let indices = [0, 2, 1, 0, 1, 3, 0, 3, 2, 1, 2, 3];

        let triangle_iter = indices.iter().map(|&i| tet[i]);

        let mesh_volume = crate::volume::triangle_mesh_signed(triangle_iter);
        assert_eq!(target_volume, mesh_volume);
        assert_eq!(mesh_volume, 1.0 / 6.0);
    }
}