1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
// SPDX-License-Identifier: MIT OR Apache-2.0
// Copyright (c) 2025 lacklustr@protonmail.com https://github.com/eadf
use crate::common::macros::{integrity_assert_eq, integrity_assert_unique, integrity_println};
use crate::corner_table::CornerIndex;
use crate::isotropic_remesh::IsotropicRemeshAlgo;
use rayon::iter::{IntoParallelIterator, ParallelIterator};
use std::fmt::Debug;
use vector_traits::num_traits::{AsPrimitive, Float};
use vector_traits::prelude::{GenericVector3, SimdUpgradable};
impl<S, V, const ENABLE_UNSAFE: bool> IsotropicRemeshAlgo<S, V, ENABLE_UNSAFE>
where
S: crate::common::sealed::ScalarType,
f64: AsPrimitive<S>,
V: Debug + Copy + From<[S; 3]> + Into<[S; 3]> + Sync + 'static,
{
/// Flip edges to improve triangle quality
/// An edge is flipped if:
/// 1. The two triangles are nearly coplanar (configurable threshold)
/// 2. The new edge would be shorter than the current edge
/// 3. The flip would not create degenerate triangles
pub(crate) fn flip_edges_aspect_ratio(&mut self, aspect_quality_weight: S) -> bool {
let num_corners = self.corner_table.vertex_of_corners().len() as u32;
let mut total_changes = false;
// Usually converges in 2-3 passes
let mut _pass = 1;
while _pass < 6 {
let mut changes_this_pass = false;
// Iterate through all edges
let mut candidates: Vec<_> = (0..num_corners)
.into_par_iter()
.map(CornerIndex)
.filter_map(|corner| {
if self.corner_table.is_corner_deleted(corner) {
return None;
}
if corner > self.corner_table.twin(corner) {
return None; // Will process from twin instead
}
self.should_flip_edge_aspect_ratio(corner, aspect_quality_weight)
})
.collect();
candidates.sort_unstable_by(|a, b| {
(a.2, a.3)
.partial_cmp(&(b.2, b.3))
.unwrap_or(std::cmp::Ordering::Equal)
});
while let Some((c0, c1, _quality, _edge)) = candidates.pop() {
if c1 != self.corner_table.twin(c0) {
// c1 is no longer twin of c0. Could this ever happen?
continue;
}
// recheck if we should flip this edge, but ignore the aspect ratio quality this time
if let Some((c0, _c1, _q, _)) = self.should_flip_edge_dihedral(c0) {
self.flip_manifold_edge(c0);
changes_this_pass = true;
#[cfg(feature = "integrity_check")]
self.check_mesh_integrity(
format!(
"after flip edge {}-{} quality:{_q}",
self.corner_table.data.dbg_corner(c0),
self.corner_table.data.dbg_corner(_c1)
)
.as_str(),
)
.unwrap();
}
}
if !changes_this_pass {
break;
}
total_changes |= changes_this_pass;
_pass += 1;
}
integrity_println!("flip_edges converged after {} passes", _pass);
total_changes
}
/// Determines if an edge should be flipped
/// ```text
///
/// c₀o
/// V₃ ─────────────── V₂ V₃ ─────────────── V₂
/// │c₀p c₀n / │ │ \ c₁n c₁ │
/// │ / │ │ \ │
/// │ t₀ / c₁ │ │ c₀p\ t₁ │
/// │ / │ => │ \ │
/// c₀no │ / │ c₁no │ \ │
/// │ / │ │ \ │
/// │c₀ / t₁ │ │ t₀ \ c₁p │
/// │ / │ │ \ │
/// │/ c₁n c₁p │ │ c₀ c₀n \ │
/// V₀ ────────────── V₁ V₀ ────────────── V₁
/// c₁o
///
/// Input: c₀: The corner index at V₀ in the V₀V₂V₃ triangle.
/// An edge should be flipped if:
/// 1. The two triangles are nearly coplanar (configurable threshold)
/// 2. The new edge would be shorter than the current edge (TODO: should this be here?)
/// 3. The flip would not create degenerate triangles
/// 4. The valence quality is improved
/// 5. The triangle quality is improved (better angles/aspect ratios)
/// ```
fn should_flip_edge_aspect_ratio(
&self,
c0: CornerIndex,
quality_threshold: S,
) -> Option<(CornerIndex, CornerIndex, S, u64)> {
// Get the four vertices of the quad formed by the two triangles
// c0 is in triangle ABC, c1 is in triangle BAD
let (c0n, c0p) = self.corner_table.next_prev(c0); // corner V₃, V₂ in t₀
let c1p = self.corner_table.opposite(c0p); // corner at V₁ in t₁
let c1 = self.corner_table.next(c1p);
integrity_assert_eq!(c1, self.corner_table.twin(c0));
let v0 = self.corner_table.vertex(c0);
let v2 = self.corner_table.vertex(c0n);
let v3 = self.corner_table.vertex(c0p);
let v1 = self.corner_table.vertex(c1p);
integrity_assert_eq!(v0, self.corner_table.vertex(self.corner_table.prev(c1p)));
// Check if all vertices are valid and distinct
integrity_assert_unique!(v0, v1, v2, v3);
// Get actual vertex positions
let pos_v0 = self.vertex(v0).to_simd();
let pos_v1 = self.vertex(v1).to_simd();
let pos_v2 = self.vertex(v2).to_simd();
let pos_v3 = self.vertex(v3).to_simd();
integrity_assert_unique!(pos_v0, pos_v1, pos_v2, pos_v3);
// Calculate current edge length (V₀V₂)
let current_edge_length_sq = (pos_v2 - pos_v0).magnitude_sq();
// Calculate potential new edge length (V₃V₁)
let new_edge_length_sq = (pos_v1 - pos_v3).magnitude_sq();
// todo: should we really do this?
// Only flip if new edge is shorter
if new_edge_length_sq >= current_edge_length_sq * 1.5_f64.as_() {
return None;
}
// Check if triangles are nearly coplanar and the new triangles are ok
if !Self::quad_is_convex_coplanar(
pos_v0,
pos_v1,
pos_v2,
pos_v3,
self.params.coplanar_threshold_sq,
) {
return None;
}
let valence_improvement = {
let valence_v0 = self.corner_table.valence(c0);
let valence_v1 = self.corner_table.valence(c1p);
let valence_v2 = self.corner_table.valence(c0n);
let valence_v3 = self.corner_table.valence(c0p);
// 5. Check valence improvement
let valence_before =
self.quad_valence_deviation(valence_v0, valence_v1, valence_v2, valence_v3, false);
let valence_after =
self.quad_valence_deviation(valence_v0, valence_v1, valence_v2, valence_v3, true);
// Hard constraint: don't worsen valence
if valence_after > valence_before {
return None;
}
valence_before - valence_after
};
// 6. Check triangle quality improvement
let quality_before =
self.min_triangle_quality_in_quad(pos_v0, pos_v1, pos_v2, pos_v3, false);
let quality_after = self.min_triangle_quality_in_quad(pos_v0, pos_v1, pos_v2, pos_v3, true);
let quality_improvement = quality_after - quality_before;
// If valence improves, prioritize by valence but use quality as tiebreaker
if valence_improvement > 0 {
// Score: valence_improvement * 1000 + quality_improvement
// Since quality is typically 0-1 range, this keeps valence dominant
let score = (valence_improvement as f64).as_() * 1000.0_f64.as_() + quality_improvement;
if Float::is_finite(score) {
return Some((c0, c1, score, self.corner_table.canonical_edge_hash(c0)));
} else {
return None;
}
}
// Valence neutral - only flip if quality improves enough
if quality_improvement > quality_before * (quality_threshold - S::ONE) {
return Some((
c0,
c1,
quality_improvement,
self.corner_table.canonical_edge_hash(c0),
));
}
None
}
fn min_triangle_quality_in_quad(
&self,
p0: S::Vec3Simd,
p1: S::Vec3Simd,
p2: S::Vec3Simd,
p3: S::Vec3Simd,
after_flip: bool,
) -> S {
let (tri1_quality, tri2_quality) = if after_flip {
// Triangles after flip: (v0, v2, v3) and (v1, v3, v2)
(
self.triangle_quality(p0, p2, p3),
self.triangle_quality(p1, p3, p2),
)
} else {
// Current triangles: (v0, v1, v2) and (v0, v3, v1)
(
self.triangle_quality(p0, p1, p2),
self.triangle_quality(p0, p3, p1),
)
};
Float::min(tri1_quality, tri2_quality)
}
fn triangle_quality(&self, p0: S::Vec3Simd, p1: S::Vec3Simd, p2: S::Vec3Simd) -> S {
let min_angle = self.min_angle(p0, p1, p2);
min_angle / 60.0_f64.to_radians().as_() // 60° is perfect
}
fn min_angle(&self, p0: S::Vec3Simd, p1: S::Vec3Simd, p2: S::Vec3Simd) -> S {
let v0 = p1 - p0;
let v0_magnitude = v0.magnitude();
let v1 = p2 - p1;
let v1_magnitude = v1.magnitude();
let v2 = p0 - p2;
let v2_magnitude = v2.magnitude();
let angle0: S = (v0.dot(-v2) / (v0_magnitude * v2_magnitude)).acos();
let angle1: S = (v1.dot(-v0) / (v1_magnitude * v0_magnitude)).acos();
let angle2: S = (v2.dot(-v1) / (v2_magnitude * v1_magnitude)).acos();
Float::min(Float::min(angle0, angle1), angle2)
}
}