axiolid_decimate/
collapse.rs

1//! Edge collapse with validity checks and a measured deviation bound.
2
3use axiolid_core::{Point3, Tolerance};
4use axiolid_mesh::TriMesh;
5use std::collections::{BTreeMap, BTreeSet};
6use thiserror::Error;
7
8/// What the caller wants back.
9#[derive(Debug, Clone, Copy, PartialEq)]
10#[non_exhaustive]
11pub enum DecimateTarget {
12    /// Reduce until at most this many triangles remain.
13    ///
14    /// The deviation bound still applies: the budget is honoured only as far
15    /// as it can be without exceeding `max_deviation`.
16    TriangleBudget(usize),
17    /// Collapse every edge whose removal stays within this deviation.
18    MaxDeviation(f64),
19}
20
21/// Why a decimation could not run.
22#[derive(Debug, Clone, PartialEq, Error)]
23#[non_exhaustive]
24pub enum DecimateError {
25    /// The index buffer is not a whole number of triangles.
26    #[error("index buffer length {0} is not a multiple of 3")]
27    RaggedIndices(usize),
28    /// A triangle references a vertex that does not exist.
29    #[error("triangle {0} references vertex {1}, which is out of range")]
30    IndexOutOfRange(usize, u32),
31    /// A deviation bound must be a positive, finite length.
32    #[error("deviation bound {0} is not a positive finite length")]
33    InvalidBound(f64),
34}
35
36/// What a decimation actually did.
37///
38/// The deviation is measured, not estimated: it is the largest distance any
39/// collapsed vertex moved from its original position. A caller that asked
40/// for 1mm accuracy can check it got it.
41#[derive(Debug, Clone, PartialEq)]
42#[non_exhaustive]
43pub struct DecimateReport {
44    /// Triangles before.
45    pub input_triangles: usize,
46    /// Triangles after.
47    pub output_triangles: usize,
48    /// Edge collapses performed.
49    pub collapses: usize,
50    /// Collapses rejected because they would have created a defect.
51    ///
52    /// A non-zero count is not a failure: it is the reason the result is
53    /// still a valid manifold rather than a faster, broken one.
54    pub rejected_unsafe: usize,
55    /// Collapses rejected because they would have exceeded the bound.
56    pub rejected_deviation: usize,
57    /// Largest distance any vertex moved, in model units.
58    pub max_deviation: f64,
59}
60
61impl DecimateReport {
62    /// Whether the mesh was left untouched.
63    pub fn is_noop(&self) -> bool {
64        self.collapses == 0
65    }
66}
67
68/// Decimate a triangle mesh by edge collapse.
69///
70/// Deterministic: candidate edges are ordered by length then by vertex
71/// index, so the same input and target produce the same output on every
72/// run, matching the determinism discipline the plan contract established
73/// in v0.6.
74///
75/// # Errors
76///
77/// Refuses a ragged index buffer, out-of-range indices, and a non-positive
78/// deviation bound.
79pub fn decimate(
80    mesh: &TriMesh,
81    target: DecimateTarget,
82    tolerance: Tolerance,
83) -> Result<(TriMesh, DecimateReport), DecimateError> {
84    if mesh.indices.len() % 3 != 0 {
85        return Err(DecimateError::RaggedIndices(mesh.indices.len()));
86    }
87    let vertex_count = mesh.positions.len();
88    for (t, chunk) in mesh.indices.chunks_exact(3).enumerate() {
89        for &index in chunk {
90            if index as usize >= vertex_count {
91                return Err(DecimateError::IndexOutOfRange(t, index));
92            }
93        }
94    }
95
96    let bound = match target {
97        DecimateTarget::MaxDeviation(d) => {
98            if !d.is_finite() || d <= 0.0 {
99                return Err(DecimateError::InvalidBound(d));
100            }
101            d
102        }
103        // A budget still needs a ceiling, or "reduce to N triangles" would
104        // licence arbitrary damage. The caller's tolerance is that ceiling.
105        DecimateTarget::TriangleBudget(_) => tolerance.linear(),
106    };
107
108    let budget = match target {
109        DecimateTarget::TriangleBudget(n) => n,
110        DecimateTarget::MaxDeviation(_) => 0,
111    };
112
113    run(mesh, budget, bound)
114}
115
116/// The collapse loop.
117///
118/// Each candidate edge is collapsed to its midpoint. The move is accepted
119/// only when it neither exceeds the deviation bound nor creates a defect.
120fn run(
121    mesh: &TriMesh,
122    budget: usize,
123    bound: f64,
124) -> Result<(TriMesh, DecimateReport), DecimateError> {
125    let input_triangles = mesh.indices.len() / 3;
126    let mut positions = mesh.positions.clone();
127    let mut triangles: Vec<[u32; 3]> = mesh
128        .indices
129        .chunks_exact(3)
130        .map(|c| [c[0], c[1], c[2]])
131        .collect();
132
133    // Deterministic candidate order: shortest edges first, ties broken by
134    // vertex index. Sorting by a float alone would leave equal-length edges
135    // in hash order and make the output depend on iteration chance.
136    let mut candidates: Vec<(u32, u32)> = unique_edges(&triangles).into_iter().collect();
137    candidates.sort_by(|a, b| {
138        let la = (positions[a.0 as usize] - positions[a.1 as usize]).length();
139        let lb = (positions[b.0 as usize] - positions[b.1 as usize]).length();
140        la.partial_cmp(&lb)
141            .unwrap_or(std::cmp::Ordering::Equal)
142            .then(a.cmp(b))
143    });
144
145    let mut report = DecimateReport {
146        input_triangles,
147        output_triangles: input_triangles,
148        collapses: 0,
149        rejected_unsafe: 0,
150        rejected_deviation: 0,
151        max_deviation: 0.0,
152    };
153    let mut moved = vec![0.0_f64; positions.len()];
154    let mut alive: Vec<bool> = vec![true; positions.len()];
155
156    for (u, v) in candidates {
157        if triangles.len() <= budget.max(4) && budget > 0 {
158            break;
159        }
160        if !alive[u as usize] || !alive[v as usize] {
161            continue;
162        }
163        let midpoint = (positions[u as usize] + positions[v as usize]) / 2.0;
164
165        // Deviation is cumulative: a vertex that already moved carries that
166        // history, so repeated collapses cannot drift past the bound one
167        // small step at a time.
168        let deviation = (midpoint - positions[u as usize])
169            .length()
170            .max((midpoint - positions[v as usize]).length())
171            + moved[u as usize].max(moved[v as usize]);
172        if deviation > bound {
173            report.rejected_deviation += 1;
174            continue;
175        }
176
177        match try_collapse(&triangles, &positions, u, v, midpoint) {
178            Some(next) => {
179                triangles = next;
180                positions[u as usize] = midpoint;
181                moved[u as usize] = deviation;
182                alive[v as usize] = false;
183                report.collapses += 1;
184                report.max_deviation = report.max_deviation.max(deviation);
185            }
186            None => report.rejected_unsafe += 1,
187        }
188    }
189
190    report.output_triangles = triangles.len();
191    Ok((compact(&triangles, &positions), report))
192}
193
194/// Attempt one collapse, returning the new triangle list or nothing.
195///
196/// Rejects the two ways a collapse damages a mesh:
197///
198/// - **Inversion.** A triangle whose normal flips has turned inside out. A
199///   decimator that permits this produces exactly the defect
200///   `axiolid-heal`'s `OrientOutward` exists to repair.
201/// - **Non-manifold edges.** Collapsing an edge whose endpoints share
202///   neighbours other than the two triangles on it welds unrelated sheets
203///   together.
204///
205/// Boundary vertices are not pinned: a collapse may move a vertex on an open
206/// boundary to an edge midpoint like any other vertex, and the deviation
207/// bound is what limits how far the silhouette moves.
208fn try_collapse(
209    triangles: &[[u32; 3]],
210    positions: &[Point3],
211    u: u32,
212    v: u32,
213    midpoint: Point3,
214) -> Option<Vec<[u32; 3]>> {
215    // The link condition: the shared neighbourhood of u and v must be
216    // exactly the two triangles on the edge. More than that, and the
217    // collapse creates a non-manifold edge.
218    let nu = neighbours(triangles, u);
219    let nv = neighbours(triangles, v);
220    let shared = nu.intersection(&nv).count();
221    if shared != 2 {
222        return None;
223    }
224
225    let mut next = Vec::with_capacity(triangles.len());
226    for &t in triangles {
227        let touches_u = t.contains(&u);
228        let touches_v = t.contains(&v);
229        if touches_u && touches_v {
230            // Degenerates to a line: this is the triangle being removed.
231            continue;
232        }
233        let mapped = t.map(|c| if c == v { u } else { c });
234        if touches_u || touches_v {
235            let before = normal(positions, t);
236            let after = normal_with(positions, mapped, u, midpoint);
237            // A flipped normal means the triangle turned inside out. Zero
238            // area after the move is equally unacceptable: it contributes
239            // nothing and confuses every downstream audit.
240            if after.length_squared() == 0.0 || before.dot(after) <= 0.0 {
241                return None;
242            }
243        }
244        next.push(mapped);
245    }
246    Some(next)
247}
248
249/// Vertices sharing a triangle with `vertex`.
250fn neighbours(triangles: &[[u32; 3]], vertex: u32) -> BTreeSet<u32> {
251    let mut set = BTreeSet::new();
252    for t in triangles {
253        if t.contains(&vertex) {
254            for &c in t {
255                if c != vertex {
256                    set.insert(c);
257                }
258            }
259        }
260    }
261    set
262}
263
264/// Every undirected edge, each appearing once.
265fn unique_edges(triangles: &[[u32; 3]]) -> BTreeSet<(u32, u32)> {
266    let mut set = BTreeSet::new();
267    for t in triangles {
268        for (a, b) in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])] {
269            set.insert((a.min(b), a.max(b)));
270        }
271    }
272    set
273}
274
275/// Unnormalised triangle normal.
276fn normal(positions: &[Point3], t: [u32; 3]) -> axiolid_core::Vec3 {
277    let a = positions[t[0] as usize];
278    let b = positions[t[1] as usize];
279    let c = positions[t[2] as usize];
280    (b - a).cross(c - a)
281}
282
283/// Triangle normal with one vertex moved to a candidate position.
284fn normal_with(
285    positions: &[Point3],
286    t: [u32; 3],
287    moved_index: u32,
288    moved_to: Point3,
289) -> axiolid_core::Vec3 {
290    let at = |i: u32| {
291        if i == moved_index {
292            moved_to
293        } else {
294            positions[i as usize]
295        }
296    };
297    let (a, b, c) = (at(t[0]), at(t[1]), at(t[2]));
298    (b - a).cross(c - a)
299}
300
301/// Drop unreferenced vertices and renumber.
302///
303/// Collapsed vertices leave holes in the position array. Emitting them would
304/// leave unreferenced vertices that make the result look defective to the
305/// v0.7 diagnosis.
306fn compact(triangles: &[[u32; 3]], positions: &[Point3]) -> TriMesh {
307    let mut remap: BTreeMap<u32, u32> = BTreeMap::new();
308    let mut kept = Vec::new();
309    let mut indices = Vec::with_capacity(triangles.len() * 3);
310    for t in triangles {
311        for &corner in t {
312            let next = u32::try_from(kept.len()).unwrap_or(u32::MAX);
313            let slot = *remap.entry(corner).or_insert_with(|| {
314                kept.push(positions[corner as usize]);
315                next
316            });
317            indices.push(slot);
318        }
319    }
320    TriMesh::new(kept, indices)
321}