axiolid_heal/
tjunction.rs

1//! Splitting triangle edges at vertices that lie on them (T-junctions).
2//!
3//! A cut opening leaves a long edge on one face matched by two short edges
4//! on the other: the mesh looks closed, but no edge is shared, so the
5//! oriented-cycle audit counts three unbalanced edges and no volume is
6//! measured. Splitting the long side at the vertex restores the pairing.
7//!
8//! Only edges the audit finds unbalanced are split, and only at vertices
9//! that are endpoints of such an edge, so a closed mesh is never touched.
10//! A vertex lies on an edge when it is strictly inside it (`0 < t < 1`),
11//! at most `tolerance.linear()` from its line and further than
12//! `tolerance.linear()` from both endpoints. No vertex is created or moved,
13//! so per-vertex data is untouched; corner data follows the split corners.
14//!
15//! The split is declined, and the triangle kept as it was, when it would
16//! give a wrong answer: the triangle's own third corner lies on its edge, a
17//! piece would be degenerate or face against the parent, or the two ends of
18//! a split side hold different values in a corner-indexed buffer, so no
19//! single value could be given to the new corner, or a per-vertex value at
20//! the inserted vertex is not what interpolating its side's endpoints gives.
21//! Declines are counted.
22//!
23//! Recommended pipeline before `axiolid_measure::cycle_volume_properties`:
24//! weld vertices, drop degenerate elements, split T-junctions, with the same
25//! tolerance.
26
27use std::collections::{BTreeMap, BTreeSet};
28use std::ops::ControlFlow;
29
30use axiolid_core::{Aabb, Point3, Tolerance};
31use axiolid_mesh::{split_triangle_at_side_corners, TriMesh};
32use axiolid_spatial::{Bvh, SpatialIndex, SpatialItem};
33
34use crate::carry;
35use crate::mesh::MeshHealError;
36
37/// What [`split_t_junctions`] did.
38#[non_exhaustive]
39#[derive(Debug, Clone, PartialEq)]
40pub struct TJunctionSplit {
41    /// Vertex-on-edge incidences found on unbalanced edges, counted once per
42    /// undirected edge however many triangles use it.
43    pub t_junctions: usize,
44    /// Triangles that were split.
45    pub split_triangles: usize,
46    /// Triangles added: each split triangle's pieces minus one.
47    pub triangles_added: usize,
48    /// Triangles with a vertex on an edge that were kept unsplit.
49    pub declined_triangles: usize,
50    /// Largest distance from a used vertex to the edge it was placed on.
51    /// The split triangles stay within this (at most the tolerance) of the
52    /// originals; it is reported, never refused.
53    pub max_deviation: f64,
54}
55
56/// Split triangle edges at vertices lying on them within `tolerance`.
57///
58/// The first piece of a split triangle keeps its index; the other pieces are
59/// appended in ascending parent order, so triangle indices stay stable.
60/// Returns the mesh unchanged when nothing is split.
61///
62/// # Errors
63/// [`MeshHealError::RaggedIndices`] when the index buffer is not a whole
64/// number of triangles.
65pub fn split_t_junctions(
66    mesh: &TriMesh,
67    tolerance: Tolerance,
68) -> Result<(TriMesh, TJunctionSplit), MeshHealError> {
69    if !mesh.indices.len().is_multiple_of(3) {
70        return Err(MeshHealError::RaggedIndices {
71            indices: mesh.indices.len(),
72        });
73    }
74    let linear = tolerance.linear();
75    let limit = linear.powi(4);
76    let count = mesh.indices.len() / 3;
77    let mut report = TJunctionSplit {
78        t_junctions: 0,
79        split_triangles: 0,
80        triangles_added: 0,
81        declined_triangles: 0,
82        max_deviation: 0.0,
83    };
84
85    // The audit's usability: valid indices, finite positions, and a doubled
86    // area above the tolerance (compared squared).
87    let usable: Vec<Option<[u32; 3]>> = (0..count)
88        .map(|t| {
89            let ids = [
90                mesh.indices[t * 3],
91                mesh.indices[t * 3 + 1],
92                mesh.indices[t * 3 + 2],
93            ];
94            let [a, b, c] = ids.map(|i| mesh.positions.get(i as usize).copied());
95            let (a, b, c) = (a?, b?, c?);
96            if !(a.is_finite() && b.is_finite() && c.is_finite()) {
97                return None;
98            }
99            ((b - a).cross(c - a).length_squared() > limit).then_some(ids)
100        })
101        .collect();
102
103    // Net direction per undirected edge over usable triangles.
104    let mut net: BTreeMap<(u32, u32), i64> = BTreeMap::new();
105    for ids in usable.iter().flatten() {
106        for s in 0..3 {
107            let (u, v) = (ids[s], ids[(s + 1) % 3]);
108            *net.entry((u.min(v), u.max(v))).or_default() += if u < v { 1 } else { -1 };
109        }
110    }
111    let unbalanced: Vec<(u32, u32)> = net
112        .iter()
113        .filter(|&(_, &winding)| winding != 0)
114        .map(|(&edge, _)| edge)
115        .collect();
116    let candidates: BTreeSet<u32> = unbalanced.iter().flat_map(|&(a, b)| [a, b]).collect();
117    if candidates.is_empty() {
118        return Ok((mesh.clone(), report));
119    }
120
121    let index = Bvh::build(
122        candidates
123            .iter()
124            .map(|&v| SpatialItem::new(v, Aabb::from_point(mesh.positions[v as usize]))),
125    );
126
127    // Vertices on each unbalanced edge, with their distance from it, in
128    // ascending (parameter, vertex index) order from the low endpoint.
129    let mut on_edge: BTreeMap<(u32, u32), Vec<(u32, f64)>> = BTreeMap::new();
130    // Parameter of each vertex along its edge, from the low endpoint.
131    let mut parameter: BTreeMap<(u32, u32, u32), f64> = BTreeMap::new();
132    for &(a, b) in &unbalanced {
133        let (pa, pb) = (mesh.positions[a as usize], mesh.positions[b as usize]);
134        let mut probe = Aabb::from_point(pa);
135        probe.extend(pb);
136        probe.min = Point3::new(
137            probe.min.x - linear,
138            probe.min.y - linear,
139            probe.min.z - linear,
140        );
141        probe.max = Point3::new(
142            probe.max.x + linear,
143            probe.max.y + linear,
144            probe.max.z + linear,
145        );
146        let ab = pb - pa;
147        let length_squared = ab.length_squared();
148        let mut found: Vec<(f64, u32, f64)> = Vec::new();
149        index.visit_aabb(&probe, &mut |&m: &u32| {
150            if m != a && m != b {
151                let pm = mesh.positions[m as usize];
152                let t = (pm - pa).dot(ab) / length_squared;
153                let distance = ab.cross(pm - pa).length() / length_squared.sqrt();
154                if t > 0.0
155                    && t < 1.0
156                    && distance <= linear
157                    && (pm - pa).length() > linear
158                    && (pm - pb).length() > linear
159                {
160                    found.push((t, m, distance));
161                }
162            }
163            ControlFlow::Continue(())
164        });
165        if found.is_empty() {
166            continue;
167        }
168        found.sort_by(|x, y| x.0.total_cmp(&y.0).then(x.1.cmp(&y.1)));
169        report.t_junctions += found.len();
170        for &(t, m, _) in &found {
171            parameter.insert((a, b, m), t);
172        }
173        on_edge.insert((a, b), found.into_iter().map(|(_, m, d)| (m, d)).collect());
174    }
175    if on_edge.is_empty() {
176        return Ok((mesh.clone(), report));
177    }
178
179    let aligned = carry::corner_buffers_aligned(mesh);
180    let mut indices = mesh.indices.clone();
181    let mut source: Vec<usize> = (0..indices.len()).collect();
182    for (t, ids) in usable.iter().enumerate() {
183        let Some(ids) = *ids else {
184            continue;
185        };
186        // Points on side `s`, from `ids[s]` to `ids[s + 1]`.
187        let sides: [Vec<(u32, f64)>; 3] = std::array::from_fn(|s| {
188            let (u, v) = (ids[s], ids[(s + 1) % 3]);
189            let mut points = on_edge
190                .get(&(u.min(v), u.max(v)))
191                .cloned()
192                .unwrap_or_default();
193            if u > v {
194                points.reverse();
195            }
196            points
197        });
198        if sides.iter().all(Vec::is_empty) {
199            continue;
200        }
201        let own_corner_on_edge =
202            (0..3).any(|s| sides[s].iter().any(|&(m, _)| m == ids[(s + 2) % 3]));
203        let corners_split_cleanly = aligned
204            && (0..3).all(|s| {
205                sides[s].is_empty() || carry::corners_agree(mesh, t * 3 + s, t * 3 + (s + 1) % 3)
206            });
207        // New corners reuse a vertex that already carries its own per-vertex
208        // data, so that data must be what interpolation would have given.
209        let vertex_data_consistent = (0..3).all(|s| {
210            let (u, v) = (ids[s], ids[(s + 1) % 3]);
211            let (low, high) = (u.min(v), u.max(v));
212            sides[s].iter().all(|&(m, _)| {
213                parameter
214                    .get(&(low, high, m))
215                    .is_some_and(|&at| carry::vertex_data_agrees(mesh, low, high, m, at))
216            })
217        });
218        let pieces = (!own_corner_on_edge && corners_split_cleanly && vertex_data_consistent)
219            .then(|| triangulate(t, ids, &sides, mesh, limit))
220            .flatten();
221        let Some(Pieces {
222            vertices,
223            corners,
224            tokens,
225        }) = pieces
226        else {
227            report.declined_triangles += 1;
228            continue;
229        };
230        let pieces = tokens.len() / 3;
231        for (slot, &token) in tokens.iter().enumerate() {
232            let (vertex, corner) = (vertices[token], corners[token]);
233            if slot < 3 {
234                indices[t * 3 + slot] = vertex;
235                source[t * 3 + slot] = corner;
236            } else {
237                indices.push(vertex);
238                source.push(corner);
239            }
240        }
241        report.split_triangles += 1;
242        report.triangles_added += pieces - 1;
243        for &(_, distance) in sides.iter().flatten() {
244            report.max_deviation = report.max_deviation.max(distance);
245        }
246    }
247    if report.split_triangles == 0 {
248        return Ok((mesh.clone(), report));
249    }
250    let mut out = mesh.clone();
251    out.indices = indices;
252    carry::remap_corners(&mut out, &source);
253    Ok((out, report))
254}
255
256/// A split triangle as tokens: `0..3` are its corners, the rest its side
257/// points in side order.
258struct Pieces {
259    /// Vertex of every token.
260    vertices: Vec<u32>,
261    /// Source corner of every token; a side point takes the corner its side
262    /// starts at.
263    corners: Vec<usize>,
264    /// The pieces as tokens, three per piece.
265    tokens: Vec<usize>,
266}
267
268/// Triangulate triangle `t` with its side points; `None` when a piece would
269/// be degenerate or face against the parent.
270fn triangulate(
271    t: usize,
272    ids: [u32; 3],
273    sides: &[Vec<(u32, f64)>; 3],
274    mesh: &TriMesh,
275    limit: f64,
276) -> Option<Pieces> {
277    let mut vertices = ids.to_vec();
278    let mut corners: Vec<usize> = (0..3).map(|i| t * 3 + i).collect();
279    let mut side_tokens: [Vec<usize>; 3] = Default::default();
280    for (s, points) in sides.iter().enumerate() {
281        for &(m, _) in points {
282            side_tokens[s].push(vertices.len());
283            vertices.push(m);
284            corners.push(t * 3 + s);
285        }
286    }
287    let mut tokens = Vec::new();
288    if !split_triangle_at_side_corners([0, 1, 2], &side_tokens, &mut tokens) || tokens.len() <= 3 {
289        return None;
290    }
291    let point = |token: usize| mesh.positions[vertices[token] as usize];
292    let parent = (point(1) - point(0)).cross(point(2) - point(0));
293    let sound = tokens.chunks(3).all(|piece| {
294        let normal = (point(piece[1]) - point(piece[0])).cross(point(piece[2]) - point(piece[0]));
295        normal.length_squared() > limit && normal.dot(parent) > 0.0
296    });
297    sound.then_some(Pieces {
298        vertices,
299        corners,
300        tokens,
301    })
302}
axiolid_heal/tjunction.rs

axiolid_heal/
tjunction.rs
