1 /*
2 normal.cpp: Helper routines for computing vertex normals
3
4 This file is part of the implementation of
5
6 Instant Field-Aligned Meshes
7 Wenzel Jakob, Daniele Panozzo, Marco Tarini, and Olga Sorkine-Hornung
8 In ACM Transactions on Graphics (Proc. SIGGRAPH Asia 2015)
9
10 All rights reserved. Use of this source code is governed by a
11 BSD-style license that can be found in the LICENSE.txt file.
12 */
13
14 #include "normal.h"
15 #include "dedge.h"
16
17 extern void
generate_smooth_normals(const MatrixXu & F,const MatrixXf & V,MatrixXf & N,bool deterministic,const ProgressCallback & progress)18 generate_smooth_normals(const MatrixXu &F, const MatrixXf &V, MatrixXf &N,
19 bool deterministic,
20 const ProgressCallback &progress) {
21 cout << "Computing vertex normals .. ";
22 cout.flush();
23
24 std::atomic<uint32_t> badFaces(0);
25 Timer<> timer;
26 N.resize(V.rows(), V.cols());
27 N.setZero();
28
29 auto map = [&](const tbb::blocked_range<uint32_t> &range) {
30 for (uint32_t f = range.begin(); f != range.end(); ++f) {
31 Vector3f fn = Vector3f::Zero();
32 for (int i=0; i<3; ++i) {
33 Vector3f v0 = V.col(F(i, f)),
34 v1 = V.col(F((i+1)%3, f)),
35 v2 = V.col(F((i+2)%3, f)),
36 d0 = v1-v0,
37 d1 = v2-v0;
38
39 if (i == 0) {
40 fn = d0.cross(d1);
41 Float norm = fn.norm();
42 if (norm < RCPOVERFLOW) {
43 badFaces++; /* degenerate */
44 break;
45 }
46 fn /= norm;
47 }
48
49 /* "Computing Vertex Normals from Polygonal Facets"
50 by Grit Thuermer and Charles A. Wuethrich, JGT 1998, Vol 3 */
51 Float angle = fast_acos(d0.dot(d1) / std::sqrt(d0.squaredNorm() * d1.squaredNorm()));
52 for (uint32_t k=0; k<3; ++k)
53 atomicAdd(&N.coeffRef(k, F(i, f)), fn[k]*angle);
54 }
55 }
56 SHOW_PROGRESS_RANGE(range, F.cols(), "Computing vertex normals (1/2)");
57 };
58
59 tbb::blocked_range<uint32_t> range(0u, (uint32_t) F.cols(), GRAIN_SIZE);
60
61 if (!deterministic)
62 tbb::parallel_for(range, map);
63 else
64 map(range);
65
66 tbb::parallel_for(
67 tbb::blocked_range<uint32_t>(0u, (uint32_t) V.cols(), GRAIN_SIZE),
68 [&](const tbb::blocked_range<uint32_t> &range) {
69 for (uint32_t i = range.begin(); i != range.end(); ++i) {
70 Float norm = N.col(i).norm();
71 if (norm < RCPOVERFLOW) {
72 N.col(i) = Vector3f::UnitX();
73 } else {
74 N.col(i) /= norm;
75 }
76 }
77 SHOW_PROGRESS_RANGE(range, V.cols(), "Computing vertex normals (2/2)");
78 }
79 );
80
81 cout << "done. (";
82 if (badFaces > 0)
83 cout << badFaces << " degenerate faces, ";
84 cout << "took " << timeString(timer.value()) << ")" << endl;
85 }
86
generate_smooth_normals(const MatrixXu & F,const MatrixXf & V,const VectorXu & V2E,const VectorXu & E2E,const VectorXb & nonManifold,MatrixXf & N,const ProgressCallback & progress)87 void generate_smooth_normals(const MatrixXu &F, const MatrixXf &V, const VectorXu &V2E,
88 const VectorXu &E2E, const VectorXb &nonManifold,
89 MatrixXf &N, const ProgressCallback &progress) {
90 cout << "Computing vertex normals .. ";
91 cout.flush();
92
93 std::atomic<uint32_t> badFaces(0);
94 Timer<> timer;
95
96 /* Compute face normals */
97 MatrixXf Nf(3, F.cols());
98
99 tbb::parallel_for(
100 tbb::blocked_range<uint32_t>(0u, (uint32_t) F.cols(), GRAIN_SIZE),
101 [&](const tbb::blocked_range<uint32_t> &range) {
102 for (uint32_t f = range.begin(); f != range.end(); ++f) {
103 Vector3f v0 = V.col(F(0, f)),
104 v1 = V.col(F(1, f)),
105 v2 = V.col(F(2, f)),
106 n = (v1-v0).cross(v2-v0);
107 Float norm = n.norm();
108 if (norm < RCPOVERFLOW) {
109 badFaces++; /* degenerate */
110 n = Vector3f::UnitX();
111 } else {
112 n /= norm;
113 }
114 Nf.col(f) = n;
115 }
116 SHOW_PROGRESS_RANGE(range, F.cols(), "Computing vertex normals (1/2)");
117 }
118 );
119
120 N.resize(3, V.cols());
121
122 /* Finally, compute the normals */
123 tbb::parallel_for(
124 tbb::blocked_range<uint32_t>(0u, (uint32_t) V.cols(), GRAIN_SIZE),
125 [&](const tbb::blocked_range<uint32_t> &range) {
126 for (uint32_t i = range.begin(); i != range.end(); ++i) {
127 uint32_t edge = V2E[i];
128 if (nonManifold[i] || edge == INVALID) {
129 N.col(i) = Vector3f::UnitX();
130 continue;
131 }
132
133 uint32_t stop = edge;
134 Vector3f normal = Vector3f::Zero();
135 do {
136 uint32_t idx = edge % 3;
137
138 Vector3f d0 = V.col(F((idx+1)%3, edge/3)) - V.col(i);
139 Vector3f d1 = V.col(F((idx+2)%3, edge/3)) - V.col(i);
140 Float angle = fast_acos(d0.dot(d1) / std::sqrt(d0.squaredNorm() * d1.squaredNorm()));
141
142 /* "Computing Vertex Normals from Polygonal Facets"
143 by Grit Thuermer and Charles A. Wuethrich, JGT 1998, Vol 3 */
144 if (std::isfinite(angle))
145 normal += Nf.col(edge/3) * angle;
146
147 uint32_t opp = E2E[edge];
148 if (opp == INVALID)
149 break;
150
151 edge = dedge_next_3(opp);
152 } while (edge != stop);
153 Float norm = normal.norm();
154 N.col(i) = norm > RCPOVERFLOW ? Vector3f(normal / norm)
155 : Vector3f::UnitX();
156 }
157 SHOW_PROGRESS_RANGE(range, V.cols(), "Computing vertex normals (2/2)");
158 }
159 );
160
161 cout << "done. (";
162 if (badFaces > 0)
163 cout << badFaces << " degenerate faces, ";
164 cout << "took " << timeString(timer.value()) << ")" << endl;
165 }
166
generate_crease_normals(MatrixXu & F,MatrixXf & V,const VectorXu & _V2E,const VectorXu & E2E,const VectorXb boundary,const VectorXb & nonManifold,Float angleThreshold,MatrixXf & N,std::map<uint32_t,uint32_t> & creases,const ProgressCallback & progress)167 void generate_crease_normals(MatrixXu &F, MatrixXf &V, const VectorXu &_V2E,
168 const VectorXu &E2E, const VectorXb boundary,
169 const VectorXb &nonManifold, Float angleThreshold,
170 MatrixXf &N, std::map<uint32_t, uint32_t> &creases,
171 const ProgressCallback &progress) {
172 const Float dpThreshold = std::cos(angleThreshold * M_PI / 180);
173
174 cout << "Computing vertex & crease normals .. ";
175 cout.flush();
176 creases.clear();
177
178 std::atomic<uint32_t> badFaces(0), creaseVert(0), offset(V.cols());
179 Timer<> timer;
180 VectorXu V2E(_V2E);
181
182 /* Compute face normals */
183 MatrixXf Nf(3, F.cols());
184
185 tbb::parallel_for(
186 tbb::blocked_range<uint32_t>(0u, (uint32_t) F.cols(), GRAIN_SIZE),
187 [&](const tbb::blocked_range<uint32_t> &range) {
188 for (uint32_t f = range.begin(); f != range.end(); ++f) {
189 Vector3f v0 = V.col(F(0, f)),
190 v1 = V.col(F(1, f)),
191 v2 = V.col(F(2, f)),
192 n = (v1-v0).cross(v2-v0);
193 Float norm = n.norm();
194 if (norm < RCPOVERFLOW) {
195 badFaces++; /* degenerate */
196 n = Vector3f::UnitX();
197 } else {
198 n /= norm;
199 }
200 Nf.col(f) = n;
201 }
202 SHOW_PROGRESS_RANGE(range, F.cols(), "Computing vertex & crease normals (1/3)");
203 }
204 );
205
206 /* Determine how many extra vertices are needed, and adjust
207 the vertex->edge pointers so that they are located just after
208 the first crease edge */
209 tbb::parallel_for(
210 tbb::blocked_range<uint32_t>(0u, (uint32_t) V.cols(), GRAIN_SIZE),
211 [&](const tbb::blocked_range<uint32_t> &range) {
212 for (uint32_t i = range.begin(); i != range.end(); ++i) {
213 uint32_t edge = V2E[i], stop = edge;
214 if (nonManifold[i] || edge == INVALID)
215 continue;
216 uint32_t creaseEdge = INVALID, nCreaseEdges = 0;
217 bool is_boundary = boundary[i];
218 do {
219 uint32_t opp = E2E[edge];
220 if (opp == INVALID)
221 break;
222 uint32_t nextEdge = dedge_next_3(opp);
223 if (Nf.col(edge / 3).dot(Nf.col(nextEdge / 3)) < dpThreshold) {
224 nCreaseEdges++;
225 if (creaseEdge == INVALID || creaseEdge < nextEdge)
226 creaseEdge = nextEdge;
227 }
228 edge = nextEdge;
229 } while (edge != stop);
230
231 if (creaseEdge != INVALID) {
232 if (!is_boundary)
233 V2E[i] = creaseEdge;
234 creaseVert += nCreaseEdges - (is_boundary ? 0 : 1);
235 }
236 }
237 SHOW_PROGRESS_RANGE(range, V.cols(), "Computing vertex & crease normals (2/3)");
238 }
239 );
240
241 uint32_t oldSize = V.cols();
242 V.conservativeResize(3, V.cols() + creaseVert);
243 N.resize(3, V.cols());
244
245 tbb::spin_mutex mutex;
246
247 /* Finally, compute the normals */
248 tbb::parallel_for(
249 tbb::blocked_range<uint32_t>(0u, oldSize, GRAIN_SIZE),
250 [&](const tbb::blocked_range<uint32_t> &range) {
251 for (uint32_t i = range.begin(); i != range.end(); ++i) {
252 uint32_t edge = V2E[i];
253 if (nonManifold[i] || edge == INVALID) {
254 N.col(i) = Vector3f::UnitX();
255 continue;
256 }
257
258 uint32_t stop = edge, vertexID = i;
259 Vector3f normal = Vector3f::Zero();
260 do {
261 uint32_t idx = edge % 3;
262 if (vertexID != i)
263 F(idx, edge/3) = vertexID;
264
265 Vector3f d0 = V.col(F((idx+1)%3, edge/3)) - V.col(i);
266 Vector3f d1 = V.col(F((idx+2)%3, edge/3)) - V.col(i);
267 Float angle = fast_acos(d0.dot(d1) / std::sqrt(d0.squaredNorm() * d1.squaredNorm()));
268
269 /* "Computing Vertex Normals from Polygonal Facets"
270 by Grit Thuermer and Charles A. Wuethrich, JGT 1998, Vol 3 */
271 if (std::isfinite(angle))
272 normal += Nf.col(edge/3) * angle;
273
274 uint32_t opp = E2E[edge];
275 if (opp == INVALID) {
276 Float norm = normal.norm();
277 N.col(vertexID) = norm > RCPOVERFLOW ? Vector3f(normal / norm)
278 : Vector3f::UnitX();
279 break;
280 }
281
282 uint32_t nextEdge = dedge_next_3(opp);
283 if (Nf.col(edge / 3).dot(Nf.col(nextEdge / 3)) < dpThreshold ||
284 nextEdge == stop) {
285 Float norm = normal.norm();
286 N.col(vertexID) = norm > RCPOVERFLOW ? Vector3f(normal / norm)
287 : Vector3f::UnitX();
288 normal = Vector3f::Zero();
289 if (nextEdge != stop) {
290 vertexID = offset++;
291 V.col(vertexID) = V.col(i);
292 mutex.lock();
293 creases[vertexID] = i;
294 mutex.unlock();
295 }
296 }
297 edge = nextEdge;
298 } while (edge != stop);
299 }
300 SHOW_PROGRESS_RANGE(range, oldSize, "Computing vertex & crease normals (3/3)");
301 }
302 );
303
304 if (offset != (uint32_t) V.cols())
305 throw std::runtime_error("Internal error (incorrect final vertex count)!");
306
307 cout << "done. (";
308 if (badFaces > 0)
309 cout << badFaces << " degenerate faces, ";
310 if (creaseVert > 0)
311 cout << creaseVert << " crease vertices, ";
312 cout << "took " << timeString(timer.value()) << ")" << endl;
313 }
314
generate_crease_normals(const MatrixXu & F,const MatrixXf & V,const VectorXu & _V2E,const VectorXu & E2E,const VectorXb boundary,const VectorXb & nonManifold,Float angleThreshold,MatrixXf & N,std::set<uint32_t> & creases,const ProgressCallback & progress)315 void generate_crease_normals(const MatrixXu &F, const MatrixXf &V,
316 const VectorXu &_V2E, const VectorXu &E2E,
317 const VectorXb boundary,
318 const VectorXb &nonManifold, Float angleThreshold,
319 MatrixXf &N, std::set<uint32_t> &creases,
320 const ProgressCallback &progress) {
321 const Float dpThreshold = std::cos(angleThreshold * M_PI / 180);
322
323 cout << "Computing vertex & crease normals .. ";
324 cout.flush();
325 creases.clear();
326
327 Timer<> timer;
328 VectorXu V2E(_V2E);
329
330 /* Compute face normals */
331 MatrixXf Nf(3, F.cols());
332 std::atomic<uint32_t> badFaces(0);
333
334 tbb::parallel_for(
335 tbb::blocked_range<uint32_t>(0u, (uint32_t) F.cols(), GRAIN_SIZE),
336 [&](const tbb::blocked_range<uint32_t> &range) {
337 for (uint32_t f = range.begin(); f != range.end(); ++f) {
338 Vector3f v0 = V.col(F(0, f)),
339 v1 = V.col(F(1, f)),
340 v2 = V.col(F(2, f)),
341 n = (v1-v0).cross(v2-v0);
342 Float norm = n.norm();
343 if (norm < RCPOVERFLOW) {
344 badFaces++; /* degenerate */
345 n = Vector3f::UnitX();
346 } else {
347 n /= norm;
348 }
349 Nf.col(f) = n;
350 }
351 SHOW_PROGRESS_RANGE(range, F.cols(), "Computing vertex & crease normals (1/3)");
352 }
353 );
354
355 /* Determine how many extra vertices are needed, and adjust
356 the vertex->edge pointers so that they are located just after
357 the first crease edge */
358 tbb::parallel_for(
359 tbb::blocked_range<uint32_t>(0u, (uint32_t) V.cols(), GRAIN_SIZE),
360 [&](const tbb::blocked_range<uint32_t> &range) {
361 for (uint32_t i = range.begin(); i != range.end(); ++i) {
362 uint32_t edge = V2E[i], stop = edge;
363 if (nonManifold[i] || edge == INVALID)
364 continue;
365 uint32_t creaseEdge = INVALID;
366 do {
367 uint32_t opp = E2E[edge];
368 if (opp == INVALID)
369 break;
370 uint32_t nextEdge = dedge_next_3(opp);
371 if (Nf.col(edge / 3).dot(Nf.col(nextEdge / 3)) < dpThreshold) {
372 if (creaseEdge == INVALID || creaseEdge < nextEdge)
373 creaseEdge = nextEdge;
374 }
375 edge = nextEdge;
376 } while (edge != stop);
377
378 if (creaseEdge != INVALID) {
379 if (!boundary[i])
380 V2E[i] = creaseEdge;
381 }
382 }
383 SHOW_PROGRESS_RANGE(range, V.cols(), "Computing vertex & crease normals (2/3)");
384 }
385 );
386
387 N.resize(3, V.cols());
388
389 /* Finally, compute the normals */
390 tbb::spin_mutex mutex;
391 tbb::parallel_for(
392 tbb::blocked_range<uint32_t>(0u, V.cols(), GRAIN_SIZE),
393 [&](const tbb::blocked_range<uint32_t> &range) {
394 for (uint32_t i = range.begin(); i != range.end(); ++i) {
395 uint32_t edge = V2E[i];
396 if (nonManifold[i] || edge == INVALID) {
397 N.col(i) = Vector3f::UnitX();
398 continue;
399 }
400
401 uint32_t stop = edge;
402 Vector3f normal = Vector3f::Zero();
403 do {
404 uint32_t idx = edge % 3, face = edge/3;
405
406 Vector3f d0 = V.col(F((idx+1)%3, face)) - V.col(i);
407 Vector3f d1 = V.col(F((idx+2)%3, face)) - V.col(i);
408 Float angle = fast_acos(d0.dot(d1) / std::sqrt(d0.squaredNorm() * d1.squaredNorm()));
409
410 /* "Computing Vertex Normals from Polygonal Facets"
411 by Grit Thuermer and Charles A. Wuethrich, JGT 1998, Vol 3 */
412 if (std::isfinite(angle))
413 normal += Nf.col(edge/3) * angle;
414
415 uint32_t opp = E2E[edge];
416 if (opp == INVALID)
417 break;
418
419 uint32_t nextEdge = dedge_next_3(opp);
420 if (Nf.col(edge / 3).dot(Nf.col(nextEdge / 3)) < dpThreshold) {
421 mutex.lock();
422 creases.insert(i);
423 mutex.unlock();
424 break;
425 }
426
427 edge = nextEdge;
428 } while (edge != stop);
429 Float norm = normal.norm();
430 N.col(i) = norm > RCPOVERFLOW ? Vector3f(normal / norm)
431 : Vector3f::UnitX();
432 }
433 SHOW_PROGRESS_RANGE(range, V.cols(), "Computing vertex & crease normals (3/3)");
434 }
435 );
436
437 cout << "done. (";
438 if (badFaces > 0)
439 cout << badFaces << " degenerate faces, ";
440 if (!creases.empty())
441 cout << creases.size() << " crease vertices, ";
442 cout << "took " << timeString(timer.value()) << ")" << endl;
443 }
444