axiolid_inspect/winding.rs
1//! Certificate that a closed oriented cycle's winding number is 0 or 1
2//! (#420, ADR 0099 D6).
3//!
4//! `axiolid_measure::cycle_volume_properties` returns the integral of the
5//! surface's integer winding number `w`, which is the enclosed volume only
6//! where `w` is 0 or 1: duplicated or overlapping shells count twice and an
7//! inverted shell subtracts. This module bounds how much volume breaks that.
8//!
9//! # Moments of the winding number
10//!
11//! The prism decomposition of [`crate::overlap`] needs no solid, only a
12//! closed oriented 2-cycle. Above a plane `z = z0` below every vertex, the
13//! winding number at a point is the signed count of faces above it,
14//!
15//! ```text
16//! w = sum over faces f of s_f * 1{ z0 < z < h_f(x, y) }
17//! ```
18//!
19//! with `s_f` the sign of the face's shadow `f'` (`+1` facing up). Then
20//!
21//! ```text
22//! integral of w = sum over f of s_f * integral over f' of (h_f - z0)
23//! integral of w^2 = sum over f, g of s_f s_g * integral over f' and g' of (min(h_f, h_g) - z0)
24//! ```
25//!
26//! The first is the signed volume and the second the unclamped pair sum
27//! behind `intersection_volume`, here of the mesh against itself over the
28//! faces as given (no flip of an inward-wound mesh, which would change the
29//! sign of `w`).
30//!
31//! # The excess
32//!
33//! `w` is an integer, so `w^2 - w = w (w - 1)` is 0 where `w` is 0 or 1 and
34//! at least 2 elsewhere. Hence
35//!
36//! ```text
37//! vol{ w not in {0, 1} } <= E := (integral of w^2 - integral of w) / 2
38//! = integral of w (w - 1) / 2
39//! ```
40//!
41//! `E` weighs a point by `w (w - 1) / 2`, which is 0 at `w` = 0 or 1, 1 at
42//! `w` = -1 or 2 and 3 at `w` = 3: it
43//! equals the volume only where `w` is in {-1..2}, and a triple-welded unit
44//! cube has `E` = 3 over a bad volume of 1. `E` is 0 exactly when `w` is 0
45//! or 1 almost everywhere. Both
46//! integrals are outward-rounded intervals, so the excess is an interval
47//! clamped below at 0; its upper end bounds the volume where the winding
48//! number leaves {0, 1}. The returned [`VolumeInterval`] brackets `E`, not
49//! that volume, so its "true volume" wording reads as `E` here. Edge-touching shells and a cavity (an inverted
50//! inner shell, `w` = 0 inside it) keep the excess at 0; a duplicated,
51//! overlapping or inverted shell makes it positive.
52
53use crate::overlap::{clamp, faces_as_given, lowest, pair_sum, signed_volume, Iv, VolumeInterval};
54use axiolid_core::Tolerance;
55use axiolid_mesh::{audit_mesh_cycle, MeshCycleHealth, TriangleMeshView};
56
57/// Why the winding number was not certified.
58#[derive(Debug, Clone, PartialEq)]
59#[non_exhaustive]
60pub enum WindingError {
61 /// A coordinate is not finite.
62 NonFinite,
63 /// The mesh is not a closed oriented cycle; the audit says why.
64 NotClosedOrientedCycle(Box<MeshCycleHealth>),
65 /// The winding excess `E = integral of w(w-1)/2` (an upper bound on the
66 /// volume with `w` outside {0, 1}, equal to it where `w` is in
67 /// {-1..2}) is certified greater than the budget.
68 Exceeds {
69 /// Bracket of `E`, its lower end above the budget.
70 excess: VolumeInterval,
71 },
72 /// The bracket straddles the budget: neither certified within nor over.
73 Undecided {
74 /// Bracket of `E`.
75 excess: VolumeInterval,
76 },
77 /// The budget is not a finite positive number.
78 InvalidBudget,
79}
80
81impl core::fmt::Display for WindingError {
82 fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
83 match self {
84 Self::NonFinite => f.write_str("a coordinate is not finite"),
85 Self::NotClosedOrientedCycle(_) => {
86 f.write_str("the mesh is not a closed oriented cycle")
87 }
88 Self::Exceeds { excess } => write!(
89 f,
90 "the winding excess (an upper bound on the volume with winding number outside 0 and 1) is at least {}, over the budget",
91 excess.lower
92 ),
93 Self::Undecided { excess } => write!(
94 f,
95 "the winding excess (an upper bound on the volume with winding number outside 0 and 1) lies in [{}, {}], straddling the budget",
96 excess.lower, excess.upper
97 ),
98 Self::InvalidBudget => f.write_str("the budget is not a finite positive number"),
99 }
100 }
101}
102
103impl core::error::Error for WindingError {}
104
105/// Certifies that the winding excess `E = integral of w(w-1)/2` of a
106/// closed oriented cycle, an upper bound on the volume with `w` not in
107/// {0, 1} that equals it where `w` is in {-1..2}, is at most `max_excess`,
108/// and returns the bracket of `E` (`upper <= max_excess`; the interval's
109/// "volume" is `E`, not the bad volume).
110///
111/// `max_excess` must be a small positive volume scaled to the model: the
112/// bracket is outward-rounded, so a budget of 0 can never be met.
113///
114/// # Errors
115///
116/// [`WindingError::InvalidBudget`] for a zero, negative or non-finite
117/// `max_excess`, [`WindingError::NonFinite`] for a non-finite position,
118/// [`WindingError::NotClosedOrientedCycle`] with the audit when the mesh is
119/// not a closed oriented cycle, [`WindingError::Exceeds`] when the excess is
120/// certified above the budget and [`WindingError::Undecided`] when its
121/// bracket straddles it.
122pub fn certify_unit_winding<M: TriangleMeshView + ?Sized>(
123 mesh: &M,
124 tolerance: Tolerance,
125 max_excess: f64,
126) -> Result<VolumeInterval, WindingError> {
127 if !max_excess.is_finite() || max_excess <= 0.0 {
128 return Err(WindingError::InvalidBudget);
129 }
130 if (0..mesh.position_count()).any(|i| !mesh.position(i).is_finite()) {
131 return Err(WindingError::NonFinite);
132 }
133 let audit = audit_mesh_cycle(mesh, tolerance);
134 if !audit.is_closed_oriented_cycle() {
135 return Err(WindingError::NotClosedOrientedCycle(Box::new(audit)));
136 }
137 let base = lowest(mesh);
138 let faces = faces_as_given(mesh);
139 let second = pair_sum(&faces, &faces, base);
140 let first = signed_volume(&faces, base);
141 let excess = clamp(second.sub(first).mul(Iv::point(0.5)), 0.0, f64::INFINITY);
142 if excess.upper <= max_excess {
143 Ok(excess)
144 } else if excess.lower > max_excess {
145 Err(WindingError::Exceeds { excess })
146 } else {
147 Err(WindingError::Undecided { excess })
148 }
149}