axiolid_measure/
winding.rs1use core::fmt;
8
9use axiolid_core::{Point3, Tolerance};
10use axiolid_mesh::{audit_mesh, MeshHealth, TriangleMeshView};
11
12#[derive(Debug, Clone, Copy, PartialEq)]
14pub struct WindingNumber {
15 pub value: f64,
17 pub skipped_singular_triangles: usize,
20}
21
22#[derive(Debug, Clone, PartialEq, Eq)]
24pub enum WindingError {
25 MeshNotWindingUsable(MeshHealth),
28 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#[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 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 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}