1 #include "dedge.hpp"
2 #include "config.hpp"
3
4 #include <atomic>
5 #include <fstream>
6 #include <iostream>
7 #include <set>
8 #include <vector>
9 #include "compare-key.hpp"
10 #ifdef WITH_TBB
11 #include "tbb/tbb.h"
12 #endif
13 namespace qflow {
14
dedge_prev(int e,int deg)15 inline int dedge_prev(int e, int deg) { return (e % deg == 0u) ? e + (deg - 1) : e - 1; }
16
atomicCompareAndExchange(volatile int * v,uint32_t newValue,int oldValue)17 inline bool atomicCompareAndExchange(volatile int* v, uint32_t newValue, int oldValue) {
18 #if defined(_WIN32)
19 return _InterlockedCompareExchange(reinterpret_cast<volatile long*>(v), (long)newValue,
20 (long)oldValue) == (long)oldValue;
21 #else
22 return __sync_bool_compare_and_swap(v, oldValue, newValue);
23 #endif
24 }
25
26 const int INVALID = -1;
27
28 #undef max
29 #undef min
compute_direct_graph(MatrixXd & V,MatrixXi & F,VectorXi & V2E,VectorXi & E2E,VectorXi & boundary,VectorXi & nonManifold)30 bool compute_direct_graph(MatrixXd& V, MatrixXi& F, VectorXi& V2E, VectorXi& E2E,
31 VectorXi& boundary, VectorXi& nonManifold) {
32 V2E.resize(V.cols());
33 V2E.setConstant(INVALID);
34
35 uint32_t deg = F.rows();
36 std::vector<std::pair<uint32_t, uint32_t>> tmp(F.size());
37
38 #ifdef WITH_TBB
39 tbb::parallel_for(
40 tbb::blocked_range<uint32_t>(0u, (uint32_t)F.cols(), GRAIN_SIZE),
41 [&](const tbb::blocked_range<uint32_t>& range) {
42 for (uint32_t f = range.begin(); f != range.end(); ++f) {
43 for (uint32_t i = 0; i < deg; ++i) {
44 uint32_t idx_cur = F(i, f), idx_next = F((i + 1) % deg, f),
45 edge_id = deg * f + i;
46 if (idx_cur >= V.cols() || idx_next >= V.cols())
47 throw std::runtime_error(
48 "Mesh data contains an out-of-bounds vertex reference!");
49 if (idx_cur == idx_next) continue;
50
51 tmp[edge_id] = std::make_pair(idx_next, INVALID);
52 if (!atomicCompareAndExchange(&V2E[idx_cur], edge_id, INVALID)) {
53 uint32_t idx = V2E[idx_cur];
54 while (!atomicCompareAndExchange((int*)&tmp[idx].second, edge_id, INVALID))
55 idx = tmp[idx].second;
56 }
57 }
58 }
59 });
60 #else
61 for (int f = 0; f < F.cols(); ++f) {
62 for (unsigned int i = 0; i < deg; ++i) {
63 unsigned int idx_cur = F(i, f), idx_next = F((i + 1) % deg, f), edge_id = deg * f + i;
64 if (idx_cur >= V.cols() || idx_next >= V.cols())
65 throw std::runtime_error("Mesh data contains an out-of-bounds vertex reference!");
66 if (idx_cur == idx_next) continue;
67
68 tmp[edge_id] = std::make_pair(idx_next, -1);
69 if (V2E[idx_cur] == -1)
70 V2E[idx_cur] = edge_id;
71 else {
72 unsigned int idx = V2E[idx_cur];
73 while (tmp[idx].second != -1) {
74 idx = tmp[idx].second;
75 }
76 tmp[idx].second = edge_id;
77 }
78 }
79 }
80 #endif
81
82 nonManifold.resize(V.cols());
83 nonManifold.setConstant(false);
84
85 E2E.resize(F.cols() * deg);
86 E2E.setConstant(INVALID);
87
88 #ifdef WITH_OMP
89 #pragma omp parallel for
90 #endif
91 for (int f = 0; f < F.cols(); ++f) {
92 for (uint32_t i = 0; i < deg; ++i) {
93 uint32_t idx_cur = F(i, f), idx_next = F((i + 1) % deg, f), edge_id_cur = deg * f + i;
94
95 if (idx_cur == idx_next) continue;
96
97 uint32_t it = V2E[idx_next], edge_id_opp = INVALID;
98 while (it != INVALID) {
99 if (tmp[it].first == idx_cur) {
100 if (edge_id_opp == INVALID) {
101 edge_id_opp = it;
102 } else {
103 nonManifold[idx_cur] = true;
104 nonManifold[idx_next] = true;
105 edge_id_opp = INVALID;
106 break;
107 }
108 }
109 it = tmp[it].second;
110 }
111
112 if (edge_id_opp != INVALID && edge_id_cur < edge_id_opp) {
113 E2E[edge_id_cur] = edge_id_opp;
114 E2E[edge_id_opp] = edge_id_cur;
115 }
116 }
117 }
118 std::atomic<uint32_t> nonManifoldCounter(0), boundaryCounter(0), isolatedCounter(0);
119
120 boundary.resize(V.cols());
121 boundary.setConstant(false);
122
123 /* Detect boundary regions of the mesh and adjust vertex->edge pointers*/
124 #ifdef WITH_OMP
125 #pragma omp parallel for
126 #endif
127 for (int i = 0; i < V.cols(); ++i) {
128 uint32_t edge = V2E[i];
129 if (edge == INVALID) {
130 isolatedCounter++;
131 continue;
132 }
133 if (nonManifold[i]) {
134 nonManifoldCounter++;
135 V2E[i] = INVALID;
136 continue;
137 }
138
139 /* Walk backwards to the first boundary edge (if any) */
140 uint32_t start = edge, v2e = INVALID;
141 do {
142 v2e = std::min(v2e, edge);
143 uint32_t prevEdge = E2E[dedge_prev(edge, deg)];
144 if (prevEdge == INVALID) {
145 /* Reached boundary -- update the vertex->edge link */
146 v2e = edge;
147 boundary[i] = true;
148 boundaryCounter++;
149 break;
150 }
151 edge = prevEdge;
152 } while (edge != start);
153 V2E[i] = v2e;
154 }
155 #ifdef LOG_OUTPUT
156 printf("counter triangle %d %d\n", (int)boundaryCounter, (int)nonManifoldCounter);
157 #endif
158 return true;
159 std::vector<std::vector<int>> vert_to_edges(V2E.size());
160 for (int i = 0; i < F.cols(); ++i) {
161 for (int j = 0; j < 3; ++j) {
162 int v = F(j, i);
163 vert_to_edges[v].push_back(i * 3 + j);
164 }
165 }
166 std::vector<int> colors(F.cols() * 3, -1);
167 bool update = false;
168 int num_v = V.cols();
169 std::map<int, int> new_vertices;
170 for (int i = 0; i < vert_to_edges.size(); ++i) {
171 int num_color = 0;
172 for (int j = 0; j < vert_to_edges[i].size(); ++j) {
173 int deid0 = vert_to_edges[i][j];
174 if (colors[deid0] == -1) {
175 int deid = deid0;
176 do {
177 colors[deid] = num_color;
178 if (num_color != 0) F(deid % 3, deid / 3) = num_v;
179 deid = deid / 3 * 3 + (deid + 2) % 3;
180 deid = E2E[deid];
181 } while (deid != deid0);
182 num_color += 1;
183 if (num_color > 1) {
184 update = true;
185 new_vertices[num_v] = i;
186 num_v += 1;
187 }
188 }
189 }
190 }
191 if (update) {
192 V.conservativeResize(3, num_v);
193 for (auto& p : new_vertices) {
194 V.col(p.first) = V.col(p.second);
195 }
196 return false;
197 }
198 return true;
199 }
200
compute_direct_graph_quad(std::vector<Vector3d> & V,std::vector<Vector4i> & F,std::vector<int> & V2E,std::vector<int> & E2E,VectorXi & boundary,VectorXi & nonManifold)201 void compute_direct_graph_quad(std::vector<Vector3d>& V, std::vector<Vector4i>& F, std::vector<int>& V2E, std::vector<int>& E2E, VectorXi& boundary, VectorXi& nonManifold) {
202 V2E.clear();
203 E2E.clear();
204 boundary = VectorXi();
205 nonManifold = VectorXi();
206 V2E.resize(V.size(), INVALID);
207
208 uint32_t deg = 4;
209 std::vector<std::pair<uint32_t, uint32_t>> tmp(F.size() * deg);
210
211 #ifdef WITH_TBB
212 tbb::parallel_for(
213 tbb::blocked_range<uint32_t>(0u, (uint32_t)F.size(), GRAIN_SIZE),
214 [&](const tbb::blocked_range<uint32_t>& range) {
215 for (uint32_t f = range.begin(); f != range.end(); ++f) {
216 for (uint32_t i = 0; i < deg; ++i) {
217 uint32_t idx_cur = F[f][i], idx_next = F[f][(i + 1) % deg],
218 edge_id = deg * f + i;
219 if (idx_cur >= V.size() || idx_next >= V.size())
220 throw std::runtime_error(
221 "Mesh data contains an out-of-bounds vertex reference!");
222 if (idx_cur == idx_next) continue;
223
224 tmp[edge_id] = std::make_pair(idx_next, INVALID);
225 if (!atomicCompareAndExchange(&V2E[idx_cur], edge_id, INVALID)) {
226 uint32_t idx = V2E[idx_cur];
227 while (!atomicCompareAndExchange((int*)&tmp[idx].second, edge_id, INVALID))
228 idx = tmp[idx].second;
229 }
230 }
231 }
232 });
233 #else
234 for (int f = 0; f < F.size(); ++f) {
235 for (unsigned int i = 0; i < deg; ++i) {
236 unsigned int idx_cur = F[f][i], idx_next = F[f][(i + 1) % deg], edge_id = deg * f + i;
237 if (idx_cur >= V.size() || idx_next >= V.size())
238 throw std::runtime_error("Mesh data contains an out-of-bounds vertex reference!");
239 if (idx_cur == idx_next) continue;
240 tmp[edge_id] = std::make_pair(idx_next, -1);
241 if (V2E[idx_cur] == -1) {
242 V2E[idx_cur] = edge_id;
243 }
244 else {
245 unsigned int idx = V2E[idx_cur];
246 while (tmp[idx].second != -1) {
247 idx = tmp[idx].second;
248 }
249 tmp[idx].second = edge_id;
250 }
251 }
252 }
253 #endif
254 nonManifold.resize(V.size());
255 nonManifold.setConstant(false);
256
257 E2E.resize(F.size() * deg, INVALID);
258
259 #ifdef WITH_OMP
260 #pragma omp parallel for
261 #endif
262 for (int f = 0; f < F.size(); ++f) {
263 for (uint32_t i = 0; i < deg; ++i) {
264 uint32_t idx_cur = F[f][i], idx_next = F[f][(i + 1) % deg], edge_id_cur = deg * f + i;
265
266 if (idx_cur == idx_next) continue;
267
268 uint32_t it = V2E[idx_next], edge_id_opp = INVALID;
269 while (it != INVALID) {
270 if (tmp[it].first == idx_cur) {
271 if (edge_id_opp == INVALID) {
272 edge_id_opp = it;
273 } else {
274 nonManifold[idx_cur] = true;
275 nonManifold[idx_next] = true;
276 edge_id_opp = INVALID;
277 break;
278 }
279 }
280 it = tmp[it].second;
281 }
282
283 if (edge_id_opp != INVALID && edge_id_cur < edge_id_opp) {
284 E2E[edge_id_cur] = edge_id_opp;
285 E2E[edge_id_opp] = edge_id_cur;
286 }
287 }
288 }
289 std::atomic<uint32_t> nonManifoldCounter(0), boundaryCounter(0), isolatedCounter(0);
290
291 boundary.resize(V.size());
292 boundary.setConstant(false);
293
294 /* Detect boundary regions of the mesh and adjust vertex->edge pointers*/
295 #ifdef WITH_OMP
296 #pragma omp parallel for
297 #endif
298 for (int i = 0; i < V.size(); ++i) {
299 uint32_t edge = V2E[i];
300 if (edge == INVALID) {
301 isolatedCounter++;
302 continue;
303 }
304 if (nonManifold[i]) {
305 nonManifoldCounter++;
306 V2E[i] = INVALID;
307 continue;
308 }
309
310 /* Walk backwards to the first boundary edge (if any) */
311 uint32_t start = edge, v2e = INVALID;
312 do {
313 v2e = std::min(v2e, edge);
314 uint32_t prevEdge = E2E[dedge_prev(edge, deg)];
315 if (prevEdge == INVALID) {
316 /* Reached boundary -- update the vertex->edge link */
317 v2e = edge;
318 boundary[i] = true;
319 boundaryCounter++;
320 break;
321 }
322 edge = prevEdge;
323 } while (edge != start);
324 V2E[i] = v2e;
325 }
326 #ifdef LOG_OUTPUT
327 printf("counter %d %d\n", (int)boundaryCounter, (int)nonManifoldCounter);
328 #endif
329 }
330
remove_nonmanifold(std::vector<Vector4i> & F,std::vector<Vector3d> & V)331 void remove_nonmanifold(std::vector<Vector4i>& F, std::vector<Vector3d>& V) {
332 typedef std::pair<uint32_t, uint32_t> Edge;
333
334 int degree = 4;
335 std::map<uint32_t, std::map<uint32_t, std::pair<uint32_t, uint32_t>>> irregular;
336 std::vector<std::set<int>> E(V.size());
337 std::vector<std::set<int>> VF(V.size());
338
339 auto kill_face_single = [&](uint32_t f) {
340 if (F[f][0] == INVALID) return;
341 for (int i = 0; i < degree; ++i) E[F[f][i]].erase(F[f][(i + 1) % degree]);
342 F[f].setConstant(INVALID);
343 };
344
345 auto kill_face = [&](uint32_t f) {
346 if (degree == 4 && F[f][2] == F[f][3]) {
347 auto it = irregular.find(F[f][2]);
348 if (it != irregular.end()) {
349 for (auto& item : it->second) {
350 kill_face_single(item.second.second);
351 }
352 }
353 }
354 kill_face_single(f);
355 };
356
357 uint32_t nm_edge = 0, nm_vert = 0;
358
359 for (uint32_t f = 0; f < (uint32_t)F.size(); ++f) {
360 if (F[f][0] == INVALID) continue;
361 if (degree == 4 && F[f][2] == F[f][3]) {
362 /* Special handling of irregular faces */
363 irregular[F[f][2]][F[f][0]] = std::make_pair(F[f][1], f);
364 continue;
365 }
366
367 bool nonmanifold = false;
368 for (uint32_t e = 0; e < degree; ++e) {
369 uint32_t v0 = F[f][e], v1 = F[f][(e + 1) % degree], v2 = F[f][(e + 2) % degree];
370 if (E[v0].find(v1) != E[v0].end() || (degree == 4 && E[v0].find(v2) != E[v0].end()))
371 nonmanifold = true;
372 }
373
374 if (nonmanifold) {
375 nm_edge++;
376 F[f].setConstant(INVALID);
377 continue;
378 }
379
380 for (uint32_t e = 0; e < degree; ++e) {
381 uint32_t v0 = F[f][e], v1 = F[f][(e + 1) % degree], v2 = F[f][(e + 2) % degree];
382
383 E[v0].insert(v1);
384 if (degree == 4) E[v0].insert(v2);
385 VF[v0].insert(f);
386 }
387 }
388
389 std::vector<Edge> edges;
390 for (auto item : irregular) {
391 bool nonmanifold = false;
392 auto face = item.second;
393 edges.clear();
394
395 uint32_t cur = face.begin()->first, stop = cur;
396 while (true) {
397 uint32_t pred = cur;
398 cur = face[cur].first;
399 uint32_t next = face[cur].first, it = 0;
400 while (true) {
401 ++it;
402 if (next == pred) break;
403 if (E[cur].find(next) != E[cur].end() && it == 1) nonmanifold = true;
404 edges.push_back(Edge(cur, next));
405 next = face[next].first;
406 }
407 if (cur == stop) break;
408 }
409
410 if (nonmanifold) {
411 nm_edge++;
412 for (auto& i : item.second) F[i.second.second].setConstant(INVALID);
413 continue;
414 } else {
415 for (auto e : edges) {
416 E[e.first].insert(e.second);
417
418 for (auto e2 : face) VF[e.first].insert(e2.second.second);
419 }
420 }
421 }
422
423 /* Check vertices */
424 std::set<uint32_t> v_marked, v_unmarked, f_adjacent;
425
426 std::function<void(uint32_t)> dfs = [&](uint32_t i) {
427 v_marked.insert(i);
428 v_unmarked.erase(i);
429
430 for (uint32_t f : VF[i]) {
431 if (f_adjacent.find(f) == f_adjacent.end()) /* if not part of adjacent face */
432 continue;
433 for (uint32_t j = 0; j < degree; ++j) {
434 uint32_t k = F[f][j];
435 if (v_unmarked.find(k) == v_unmarked.end() || /* if not unmarked OR */
436 v_marked.find(k) != v_marked.end()) /* if already marked */
437 continue;
438 dfs(k);
439 }
440 }
441 };
442
443 for (uint32_t i = 0; i < (uint32_t)V.size(); ++i) {
444 v_marked.clear();
445 v_unmarked.clear();
446 f_adjacent.clear();
447
448 for (uint32_t f : VF[i]) {
449 if (F[f][0] == INVALID) continue;
450
451 for (uint32_t k = 0; k < degree; ++k) v_unmarked.insert(F[f][k]);
452
453 f_adjacent.insert(f);
454 }
455
456 if (v_unmarked.empty()) continue;
457 v_marked.insert(i);
458 v_unmarked.erase(i);
459
460 dfs(*v_unmarked.begin());
461
462 if (v_unmarked.size() > 0) {
463 nm_vert++;
464 for (uint32_t f : f_adjacent) kill_face(f);
465 }
466 }
467
468 if (nm_vert > 0 || nm_edge > 0) {
469 std::cout << "Non-manifold elements: vertices=" << nm_vert << ", edges=" << nm_edge
470 << std::endl;
471 }
472 uint32_t nFaces = 0, nFacesOrig = F.size();
473 for (uint32_t f = 0; f < (uint32_t)F.size(); ++f) {
474 if (F[f][0] == INVALID) continue;
475 if (nFaces != f) {
476 F[nFaces] = F[f];
477 }
478 ++nFaces;
479 }
480
481 if (nFacesOrig != nFaces) {
482 F.resize(nFaces);
483 std::cout << "Faces reduced from " << nFacesOrig << " -> " << nFaces << std::endl;
484 }
485 }
486
487 } // namespace qflow
488