axiolid_measure/
winding.rs

1//! Generalized winding-number accumulation over validated triangle views.
2//!
3//! This module reports a raw metric-like winding measure; it does not classify
4//! points as inside or outside. Callers retain shell-closure requirements,
5//! orientation policy, classification thresholds, and application semantics.
6
7use core::fmt;
8
9use axiolid_core::{Point3, Tolerance};
10use axiolid_mesh::{audit_mesh, MeshHealth, TriangleMeshView};
11
12/// Raw winding accumulation for one query point.
13#[derive(Debug, Clone, Copy, PartialEq)]
14pub struct WindingNumber {
15    /// Sum of oriented triangle solid angles divided by `4π`.
16    pub value: f64,
17    /// Faces ignored because the query point was within the supplied linear
18    /// tolerance of one of their vertices, where their solid angle is undefined.
19    pub skipped_singular_triangles: usize,
20}
21
22/// Failure to prepare or evaluate a winding query.
23#[derive(Debug, Clone, PartialEq, Eq)]
24pub enum WindingError {
25    /// The mesh contains invalid indices, non-finite coordinates, no usable
26    /// triangles, or degenerate triangles under the supplied tolerance.
27    MeshNotWindingUsable(MeshHealth),
28    /// Winding is undefined for a non-finite query point.
29    NonFinitePoint,
30}
31
32impl fmt::Display for WindingError {
33    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
34        match self {
35            Self::MeshNotWindingUsable(_) => {
36                formatter.write_str("mesh is not usable for winding accumulation")
37            }
38            Self::NonFinitePoint => formatter.write_str("winding query point must be finite"),
39        }
40    }
41}
42
43impl std::error::Error for WindingError {}
44
45/// A read-only, structurally prepared mesh for repeated winding queries.
46///
47/// Preparation performs the O(triangles) structural audit once. Individual
48/// queries only read caller-owned geometry and have no global state, mutation,
49/// threading assumption, or source identity. A batch/provider API can reuse
50/// this exact prepared-input contract later.
51#[derive(Clone, Copy)]
52pub struct WindingMesh<'a, M: TriangleMeshView + ?Sized> {
53    mesh: &'a M,
54    tolerance: Tolerance,
55}
56
57impl<M: TriangleMeshView + ?Sized> fmt::Debug for WindingMesh<'_, M> {
58    fn fmt(&self, formatter: &mut fmt::Formatter<'_>) -> fmt::Result {
59        formatter
60            .debug_struct("WindingMesh")
61            .field("positions", &self.mesh.position_count())
62            .field("triangles", &self.mesh.triangle_count())
63            .field("tolerance", &self.tolerance)
64            .finish()
65    }
66}
67
68impl<'a, M: TriangleMeshView + ?Sized> WindingMesh<'a, M> {
69    /// Validate and prepare a foreign or owned triangle mesh.
70    ///
71    /// Closedness is intentionally not required: generalized winding returns a
72    /// continuous raw value for open surfaces too. Degenerate faces are rejected
73    /// because their solid angles are not defined by this operation.
74    pub fn prepare(mesh: &'a M, tolerance: Tolerance) -> Result<Self, WindingError> {
75        let health = audit_mesh(mesh, tolerance);
76        if !health.is_surface_usable() || health.degenerate_triangles != 0 {
77            return Err(WindingError::MeshNotWindingUsable(health));
78        }
79        Ok(Self { mesh, tolerance })
80    }
81
82    /// Accumulate the generalized winding number at one finite query point.
83    pub fn winding_number(&self, point: Point3) -> Result<WindingNumber, WindingError> {
84        if !point.is_finite() {
85            return Err(WindingError::NonFinitePoint);
86        }
87
88        let mut solid_angle = 0.0;
89        let mut skipped_singular_triangles = 0;
90        let singular_distance_squared = self.tolerance.linear().powi(2);
91
92        for triangle_index in 0..self.mesh.triangle_count() {
93            let indices = self
94                .mesh
95                .triangle(triangle_index)
96                .map(|index| index as usize);
97            let [a, b, c] = indices.map(|index| self.mesh.position(index) - point);
98            let squared_lengths = [a.length_squared(), b.length_squared(), c.length_squared()];
99            let singular = squared_lengths.iter().any(|&squared_length| {
100                squared_length == 0.0
101                    || (self.tolerance.linear() > 0.0 && squared_length < singular_distance_squared)
102            });
103            if singular {
104                skipped_singular_triangles += 1;
105                continue;
106            }
107
108            let [length_a, length_b, length_c] = squared_lengths.map(f64::sqrt);
109            let numerator = a.dot(b.cross(c));
110            let denominator = length_a * length_b * length_c
111                + a.dot(b) * length_c
112                + b.dot(c) * length_a
113                + c.dot(a) * length_b;
114            solid_angle += 2.0 * numerator.atan2(denominator);
115        }
116
117        Ok(WindingNumber {
118            value: solid_angle / (4.0 * std::f64::consts::PI),
119            skipped_singular_triangles,
120        })
121    }
122}