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
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
//! Algorithms for subdividing a mesh.
use std::{
collections::HashMap,
};
use lina::Point3;
use crate::{
prelude::*,
cast, hsize,
map::{DenseMap, SparseMap, set::DenseSet},
util::{PrimitiveFloat, Pos3Like},
};
/// The √3 subdivision algorithm.
///
///
/// # References
///
/// Kobbelt, Leif. "√ 3-subdivision." Proceedings of the 27th annual conference
/// on Computer graphics and interactive techniques. 2000.
///
/// <https://www.graphics.rwth-aachen.de/media/papers/sqrt31.pdf>
pub fn sqrt3<MeshT, MapT, ScalarT>(
mesh: &mut MeshT,
vertex_positions: &mut MapT,
num_iterations: u32,
)
where
MeshT: TriMesh + EdgeMesh + MeshMut + EdgeAdj, // TODO: doesn't need to be tri in the beginning
MapT: PropStoreMut<VertexHandle>,
MapT::Target: Pos3Like<Scalar = ScalarT>,
ScalarT: PrimitiveFloat,
{
for i in 0..num_iterations {
sqrt3_impl(mesh, vertex_positions, i % 2 == 1);
}
}
// `split_boundary = true` has the condition that all boundary faces must only
// have one boundary edge!
#[inline(never)]
fn sqrt3_impl<MeshT, MapT, ScalarT>(
mesh: &mut MeshT,
vertex_positions: &mut MapT,
split_boundary: bool,
)
where
MeshT: TriMesh + EdgeMesh + MeshMut + EdgeAdj, // TODO: doesn't need to be tri in the beginning
MapT: PropStoreMut<VertexHandle>,
MapT::Target: Pos3Like<Scalar = ScalarT>,
ScalarT: PrimitiveFloat,
{
// Helper macro to create literal values of type `ScalarT`
macro_rules! lit {
($x:literal) => (cast::lossless::<f32, ScalarT>($x));
}
/// A simple memoization helper to avoid calculating `alpha` all the time
/// as it's only dependent on the valence of a vertex. The number of
/// different valences is usually pretty low.
struct AlphaCache<ScalarT: PrimitiveFloat> {
low: [ScalarT; 10],
// TODO: use a faster hash function
rest: HashMap<hsize, ScalarT>,
}
impl<ScalarT: PrimitiveFloat> AlphaCache<ScalarT> {
fn alpha_for(valence: hsize) -> ScalarT {
(
lit!(4.0) - lit!(2.0) * (
lit!(2.0) * ScalarT::PI() / cast::lossy::<_, ScalarT>(valence)
).cos()
) / lit!(9.0)
}
fn new() -> Self {
Self {
low: [
lit!(0.0),
Self::alpha_for(1),
Self::alpha_for(2),
Self::alpha_for(3),
Self::alpha_for(4),
Self::alpha_for(5),
Self::alpha_for(6),
Self::alpha_for(7),
Self::alpha_for(8),
Self::alpha_for(9),
],
rest: HashMap::new(),
}
}
fn get(&mut self, valence: hsize) -> ScalarT {
if valence < 10 {
self.low[valence as usize]
} else {
self.get_from_hashmap(valence)
}
}
#[cold]
#[inline(never)]
fn get_from_hashmap(&mut self, valence: hsize) -> ScalarT {
*self.rest.entry(valence).or_insert_with(|| Self::alpha_for(valence))
}
}
// ----- (1) Calculate positions for new boundary vertices -----------------------------------
// Remember the original edges of the mesh and calculate boundary vertex
// positions if we will split the boundaries. We have to do this now
// because later we will already change the topology and can't properly
// calculate the positions anymore.
let mut old_edges = DenseSet::with_capacity(mesh.num_edges());
let mut new_boundary_points = SparseMap::new();
for e in mesh.edges() {
if e.is_boundary() {
if split_boundary {
// Here we prepare the positions for the two new vertices that
// will be created by splitting the boundary face later on.
//
// ll x rr ll x rr
// \ / \ / \ /| |\ /
// \ / \ / \ / | | \ /
// \ / F \ / \ / / \ \ /
// \ / \ / \ / | | \ /
// l --------- r l - a - b - r
// e
//
let [l, r] = e.endpoints();
let ll = l.adjacent_edges()
.find(|en| en.is_boundary() && en.handle() != e.handle())
.map(|en| en.opposite_endpoint_of(l.handle()))
.unwrap();
let rr = r.adjacent_edges()
.find(|en| en.is_boundary() && en.handle() != e.handle())
.map(|en| en.opposite_endpoint_of(r.handle()))
.unwrap();
let pos_l = vertex_positions[l.handle()].to_point3().to_vec();
let pos_r = vertex_positions[r.handle()].to_point3().to_vec();
let pos_ll = vertex_positions[ll.handle()].to_point3().to_vec();
let pos_rr = vertex_positions[rr.handle()].to_point3().to_vec();
let pos_a = Point3::origin()
+ (pos_r * lit!(10.0) + pos_l * lit!(16.0) + pos_ll) / lit!(27.0);
let pos_b = Point3::origin()
+ (pos_l * lit!(10.0) + pos_r * lit!(16.0) + pos_rr) / lit!(27.0);
new_boundary_points.insert(e.handle(), [pos_a, pos_b]);
}
} else {
old_edges.insert(e.handle());
}
}
// ----- (2) Calculate new positions for old vertices ----------------------------------------
// We have to calculate a new position for all already existing vertices
// (except boundary vertices!). To do that we need their old positions, so
// we have no choice but making a copy. We write the new positions into
// this copy and only write them back into `vertex_positions` at the very
// end, since the calculation of new vertex points also relies on the old
// positions.
let mut alphas = AlphaCache::new();
let new_positions = mesh.vertices().map(|v| {
let vh = v.handle();
let old_pos = vertex_positions[vh];
let new_pos = if v.is_boundary() {
if !split_boundary || v.is_isolated() {
old_pos
} else {
// Instead of taking the centroid of all adjacent vertices,
// only the two adjacent boundary vertices are used. We checked
// that this vertex is boundary and not isolated, so it has at
// least two adjacent boundary vertices. If it has more, it's a
// multi-blade vertex which we don't allow (right now).
let mut neighbors = v.adjacent_edges()
.filter(|e| e.is_boundary())
.map(|e| e.opposite_endpoint_of(vh));
let neighbor_a = neighbors.next().unwrap();
let neighbor_b = neighbors.next().unwrap();
if neighbors.next().is_some() {
// TODO: add that to doc comment
panic!(
"encountered multi-blade vertex {:?} in sqrt3 (these are not allowed)",
vh,
);
}
let pos_a = vertex_positions[neighbor_a.handle()];
let pos_b = vertex_positions[neighbor_b.handle()];
// The neighbor vertices are weighted with 4/27 and the old
// vertex with 19/27. This is from formula (9) in chapter 5.
MapT::Target::from_coords(
(lit!(4.0) * (pos_a.x() + pos_b.x()) + lit!(19.0) * old_pos.x()) / lit!(27.0),
(lit!(4.0) * (pos_a.y() + pos_b.y()) + lit!(19.0) * old_pos.y()) / lit!(27.0),
(lit!(4.0) * (pos_a.z() + pos_b.z()) + lit!(19.0) * old_pos.z()) / lit!(27.0),
)
}
} else {
// Count the number of neighbors and calculate the centroid of all
// neighbors.
let mut valence = 0;
let centroid = v.adjacent_vertices()
.inspect(|_| valence += 1)
.map(|v| vertex_positions[v.handle()])
.centroid()
.unwrap(); // we checked the vertex is not a boundary vertex
let alpha = alphas.get(valence);
// Lerp between `old_pos` and `centroid` by the amount `alpha`
MapT::Target::from_coords(
(lit!(1.0) - alpha) * old_pos.x() + alpha * centroid.x(),
(lit!(1.0) - alpha) * old_pos.y() + alpha * centroid.y(),
(lit!(1.0) - alpha) * old_pos.z() + alpha * centroid.z(),
)
};
(vh, new_pos)
}).collect::<DenseMap<_, _>>();
// ----- (3) Split faces and calc new vertex positions ---------------------------------------
// We create a new vertex per face by splitting each face into three new
// ones. First we can reserve a bunch of memory, because we know exactly
// how many elements we will end up with.
vertex_positions.reserve(mesh.num_faces());
mesh.reserve_for_vertices(mesh.num_faces());
mesh.reserve_for_faces(mesh.num_faces() * 3);
let mut it = mesh.face_handles_mut();
while let Some(fh) = it.next() {
let f = it.mesh().get_ref(fh);
if split_boundary {
// Check if the face has a boundary edge. We know that if it has,
// there is only one such edge (function contract).
let boundary_edge = f.adjacent_edges().find(|e| e.is_boundary());
if let Some(boundary_edge) = boundary_edge {
//
// ll x rr ll x rr
// \ / \ / \ /| |\ /
// \ / \ / \ / | | \ /
// \ / F \ / \ / / \ \ /
// \ / \ / \ / | | \ /
// l --------- r l - a - b - r
//
// boundary_edge
//
let [_, r] = boundary_edge.endpoints();
let rh = r.handle();
let eh = boundary_edge.handle();
// Retrieve the new points we calculated before
let [pos_a, pos_b] = new_boundary_points[eh];
// Split the boundary edge/face once
let res = it.mesh().split_edge_with_faces(eh);
let va = res.vertex;
// Decide what edge we need to split next.
let other_edge = res.replacement_edges
.iter()
.copied()
.find(|&new_edge| it.mesh().endpoints_of_edge(new_edge).contains(&rh))
.unwrap();
// Split another edge
let res = it.mesh().split_edge_with_faces(other_edge);
let vb = res.vertex;
vertex_positions.insert(va, pos_a.convert());
vertex_positions.insert(vb, pos_b.convert());
continue;
}
}
// The position of the new vertex is just the centroid of the
// face's vertices. We can unwrap because the face always has three
// vertices.
let point_pos = f
.adjacent_vertices()
.map(|v| vertex_positions[v.handle()])
.centroid()
.unwrap();
// Split face and set position of the midpoint.
let vh = it.mesh().split_face(fh);
vertex_positions.insert(vh, point_pos);
}
// ----- (4) Flip all old edges and set relaxed vertex positions -----------------------------
for eh in old_edges.handles() {
mesh.flip_edge(eh);
}
for vh in new_positions.handles() {
vertex_positions[vh] = new_positions[vh];
}
}