TTK
Loading...
Searching...
No Matches
FiberSurface.cpp
Go to the documentation of this file.
1#include <FiberSurface.h>
2
3using namespace std;
4using namespace ttk;
5
6static const float PREC_FLT{powf(10.F, -FLT_DIG)};
7static const float PREC_FLT_2{powf(10.F, -FLT_DIG + 2)};
8static const double PREC_DBL{Geometry::pow(10.0, -DBL_DIG)};
9static const double PREC_DBL_4{Geometry::pow(10.0, -DBL_DIG + 4)};
10
12 this->setDebugMsgPrefix("FiberSurface");
13}
14
16 const SimplexId &tetId,
17 const SimplexId &triangleId0,
18 const SimplexId &triangleId1,
19 const vector<vector<IntersectionTriangle>> &tetIntersections) const {
20
21 // if the two triangles have at least 2 vertices in common, they are adjacent.
22 SimplexId commonVertexNumber = 0;
23
24 for(int i = 0; i < 3; i++) {
25 std::array<double, 3> p0{};
26
27 for(int j = 0; j < 3; j++) {
28 p0[j] = tetIntersections[tetId][triangleId0].p_[i][j];
29 }
30
31 // check if this guy exists in the other triangle
32 for(int j = 0; j < 3; j++) {
33 std::array<double, 3> p1{};
34
35 bool isTheSame = true;
36 for(int k = 0; k < 3; k++) {
37 p1[k] = tetIntersections[tetId][triangleId1].p_[j][k];
38
39 if(fabs(p0[k] - p1[k]) > PREC_FLT) {
40 isTheSame = false;
41 break;
42 }
43 }
44 if(isTheSame) {
45 commonVertexNumber++;
46 break;
47 }
48 }
49 }
50
51 return commonVertexNumber;
52}
53
55 const SimplexId &tetId,
56 const SimplexId &triangleId,
57 const pair<double, double> &intersection,
58 const vector<vector<IntersectionTriangle>> &tetIntersections,
59 std::array<double, 3> &pA,
60 std::array<double, 3> &pB,
61 SimplexId &pivotVertexId,
62 bool &edgeFiber) const {
63
64 pivotVertexId = -1;
65
66 // let's check first if an edge coincide with the fiber
67 for(int i = 0; i < 3; i++) {
68 if((fabs(intersection.first
69 - tetIntersections[tetId][triangleId].uv_[i].first)
70 < PREC_DBL_4)
71 && fabs(intersection.second
72 - tetIntersections[tetId][triangleId].uv_[i].second)
73 < PREC_DBL_4
74 && fabs(intersection.first
75 - tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first)
76 < PREC_DBL_4
77 && fabs(intersection.second
78 - tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second)
79 < PREC_DBL_4) {
80 // edge 0 - 1 is on the fiber. the pivot is 2
81 pivotVertexId = (i + 2) % 3;
82 edgeFiber = true;
83 break;
84 }
85 }
86
87 // NOTE: it'd be better here to compute the actual distance
88 if(pivotVertexId == -1) {
89 for(int i = 0; i < 3; i++) {
90 if(tetIntersections[tetId][triangleId].uv_[i].first
91 > intersection.first) {
92 if((tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first
93 <= intersection.first)
94 && (tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].first
95 <= intersection.first)) {
96 pivotVertexId = i;
97 break;
98 }
99 }
100 if(tetIntersections[tetId][triangleId].uv_[i].first
101 < intersection.first) {
102 if((tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first
103 >= intersection.first)
104 && (tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].first
105 >= intersection.first)) {
106 pivotVertexId = i;
107 break;
108 }
109 }
110 }
111 }
112 if(pivotVertexId == -1) {
113 for(int i = 0; i < 3; i++) {
114 if(tetIntersections[tetId][triangleId].uv_[i].second
115 > intersection.second) {
116 if((tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second
117 <= intersection.second)
118 && (tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].second
119 <= intersection.second)) {
120 pivotVertexId = i;
121 break;
122 }
123 }
124 if(tetIntersections[tetId][triangleId].uv_[i].second
125 < intersection.second) {
126 if((tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second
127 >= intersection.second)
128 && (tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].second
129 >= intersection.second)) {
130 pivotVertexId = i;
131 break;
132 }
133 }
134 }
135 }
136
137 if(pivotVertexId == -1) {
138 // initially, several triangles can intersect.
139 // after a few iterations, a valid triangle for re-meshing may be
140 // subdivided into several sub-triangles. not all of them may
141 // intersect with a given query (pivotVertexId == -1)
142 return -1;
143 }
144
145 // compute the interpolations
146 std::array<double, 2> baryCentrics0{}, baryCentrics1{};
147 std::array<double, 2> p{}, p0{}, p1{}, p2{};
148
149 p[0] = intersection.first;
150 p[1] = intersection.second;
151
152 p0[0] = tetIntersections[tetId][triangleId].uv_[pivotVertexId].first;
153 p0[1] = tetIntersections[tetId][triangleId].uv_[pivotVertexId].second;
154
155 p1[0]
156 = tetIntersections[tetId][triangleId].uv_[(pivotVertexId + 1) % 3].first;
157 p1[1]
158 = tetIntersections[tetId][triangleId].uv_[(pivotVertexId + 1) % 3].second;
159
160 p2[0]
161 = tetIntersections[tetId][triangleId].uv_[(pivotVertexId + 2) % 3].first;
162 p2[1]
163 = tetIntersections[tetId][triangleId].uv_[(pivotVertexId + 2) % 3].second;
164
166 p0.data(), p1.data(), p.data(), baryCentrics0, 2);
168 p0.data(), p2.data(), p.data(), baryCentrics1, 2);
169
170 for(int i = 0; i < 3; i++) {
171 pA[i] = baryCentrics0[0]
172 * tetIntersections[tetId][triangleId].p_[pivotVertexId][i]
173 + baryCentrics0[1]
174 * tetIntersections[tetId][triangleId]
175 .p_[(pivotVertexId + 1) % 3][i];
176 }
177
178 for(int i = 0; i < 3; i++) {
179 pB[i] = baryCentrics1[0]
180 * tetIntersections[tetId][triangleId].p_[pivotVertexId][i]
181 + baryCentrics1[1]
182 * tetIntersections[tetId][triangleId]
183 .p_[(pivotVertexId + 2) % 3][i];
184 }
185
186 return 0;
187}
188
190 const SimplexId &tetId,
191 const SimplexId &triangleId0,
192 const SimplexId &triangleId1,
193 const SimplexId &polygonEdgeId0,
194 const SimplexId &polygonEdgeId1,
195 const std::pair<double, double> &intersection,
196 SimplexId &newVertexNumber,
197 SimplexId &newTriangleNumber,
198 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
199 std::vector<std::vector<Vertex>> &tetNewVertices) const {
200
201 SimplexId commonVertexNumber = getNumberOfCommonVertices(
202 tetId, triangleId0, triangleId1, tetIntersections);
203
204 // make sure the two triangles are not already adjacent
205 if(commonVertexNumber == 2) {
206 // NOTE: here, we used to quit with only one.
207 // however, in high res data-sets you can have triangles that intersect
208 // and share one vertex.
209 return -1;
210 }
211
212 SimplexId pivotVertexIda = -1, pivotVertexIdb = -1;
213 std::array<double, 3> p0a{}, p1a{}, p0b{}, p1b{};
214
215 // extract the fiber in both triangles and see if they match up
216 bool edgeFiber0 = false;
217 computeTriangleFiber(tetId, triangleId0, intersection, tetIntersections, p0a,
218 p1a, pivotVertexIda, edgeFiber0);
219
220 bool edgeFiber1 = false;
221 computeTriangleFiber(tetId, triangleId1, intersection, tetIntersections, p0b,
222 p1b, pivotVertexIdb, edgeFiber1);
223
224 // if((commonVertexNumber == 1)&&
225 // ((edgeFiber0)||(edgeFiber1))){
226 // // case of adjacent triangles along a vertex.
227 // // one of the two triangles has a fiber along an edge.
228 // return -2;
229 // }
230
231 if(p0a.size() != 3)
232 return -1;
233 if(p1a.size() != 3)
234 return -2;
235 if(p0b.size() != 3)
236 return -3;
237 if(p1b.size() != 3)
238 return -4;
239
240 // we need to make sure p0a and p1a are not the same (vertex case)
241 bool vertexA = false;
242 bool vertexB = false;
243 if((fabs(p0a[0] - p1a[0]) < PREC_DBL) && (fabs(p0a[1] - p1a[1]) < PREC_DBL)
244 && (fabs(p0a[2] - p1a[2]) < PREC_DBL)) {
245 vertexA = true;
246 }
247 if((fabs(p0b[0] - p1b[0]) < PREC_DBL) && (fabs(p0b[1] - p1b[1]) < PREC_DBL)
248 && (fabs(p0b[2] - p1b[2]) < PREC_DBL)) {
249 vertexB = true;
250 }
251 if((vertexA) || (vertexB)) {
252 // intersection on triangle vertices
253 return -1;
254 }
255
256 bool foundA = false, foundB = false;
257 std::array<double, 3> pA{}, pB{};
258
259 // test if p0a lies in [p0b, p1b]
260 if(Geometry::isPointOnSegment(p0a.data(), p0b.data(), p1b.data())) {
261 // p0a is in the segment
262 pA = p0a;
263 foundA = true;
264 }
265
266 // test if p1a lies in [p0b, p1b]
267 if(Geometry::isPointOnSegment(p1a.data(), p0b.data(), p1b.data())) {
268 if(!foundA) {
269 pA = p1a;
270 foundA = true;
271 } else if(!foundB) {
272 // check it's far enough from pA
273 if((fabs(pA[0] - p1a[0]) > PREC_DBL_4)
274 || (fabs(pA[1] - p1a[1]) > PREC_DBL_4)
275 || (fabs(pA[2] - p1a[2]) > PREC_DBL_4)) {
276 pB = p1a;
277 foundB = true;
278 }
279 }
280 }
281
282 // test if p0b lies in [p0a, p1a]
283 if(Geometry::isPointOnSegment(p0b.data(), p0a.data(), p1a.data())) {
284 if(!foundA) {
285 pA = p0b;
286 foundA = true;
287 } else if(!foundB) {
288 // check it's far enough from pA
289 if((fabs(pA[0] - p0b[0]) > PREC_DBL_4)
290 || (fabs(pA[1] - p0b[1]) > PREC_DBL_4)
291 || (fabs(pA[2] - p0b[2]) > PREC_DBL_4)) {
292 pB = p0b;
293 foundB = true;
294 }
295 }
296 }
297
298 // test if p1b lies in [p0a, p1a]
299 if(Geometry::isPointOnSegment(p1b.data(), p0a.data(), p1a.data())) {
300 if(!foundA) {
301 pA = p1b;
302 } else if(!foundB) {
303 // check it's far enough from pA
304 if((fabs(pA[0] - p1b[0]) > PREC_DBL_4)
305 || (fabs(pA[1] - p1b[1]) > PREC_DBL_4)
306 || (fabs(pA[2] - p1b[2]) > PREC_DBL_4)) {
307 pB = p1b;
308 }
309 }
310 }
311
312 if((!pA.size()) || (!pB.size()))
313 return -2;
314
315 if(!edgeFiber0) {
317 tetId, triangleId0, polygonEdgeId0, intersection, pA, pB, pivotVertexIda,
318 newVertexNumber, newTriangleNumber, tetIntersections, tetNewVertices);
319 }
320
321 if(!edgeFiber1) {
323 tetId, triangleId1, polygonEdgeId1, intersection, pA, pB, pivotVertexIdb,
324 newVertexNumber, newTriangleNumber, tetIntersections, tetNewVertices);
325 }
326
327 return 0;
328}
329
331 const SimplexId &tetId,
332 const SimplexId &triangleId,
333 const SimplexId &polygonEdgeId,
334 const std::pair<double, double> &intersection,
335 const std::array<double, 3> &pA,
336 const std::array<double, 3> &pB,
337 const SimplexId &pivotVertexId,
338 SimplexId &newVertexNumber,
339 SimplexId &newTriangleNumber,
340 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
341 std::vector<std::vector<Vertex>> &tetNewVertices) const {
342
343 // check if the triangle has already been intersected on that fiber
344 if((fabs(tetIntersections[tetId][triangleId].intersection_.first
345 - intersection.first)
346 < PREC_FLT)
347 && (fabs(tetIntersections[tetId][triangleId].intersection_.second
348 - intersection.second)
349 < PREC_FLT)) {
350
351 return -2;
352 }
353
354 // check if the fiber is on an edge or not (if so we stop)
355 for(int i = 0; i < 3; i++) {
356 if((fabs(intersection.first
357 - tetIntersections[tetId][triangleId].uv_[i].first)
358 < PREC_FLT)
359 && fabs(intersection.second
360 - tetIntersections[tetId][triangleId].uv_[i].second)
361 < PREC_FLT
362 && fabs(intersection.first
363 - tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first)
364 < PREC_FLT
365 && fabs(intersection.second
366 - tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second)
367 < PREC_FLT) {
368 return -3;
369 }
370 }
371
372 // 1. compute the barycentric coordinates of pA and pB
373 std::array<double, 3> barypA{}, barypB{};
375 tetIntersections[tetId][triangleId].p_[0].data(),
376 tetIntersections[tetId][triangleId].p_[1].data(),
377 tetIntersections[tetId][triangleId].p_[2].data(), pA.data(), barypA);
378
380 tetIntersections[tetId][triangleId].p_[0].data(),
381 tetIntersections[tetId][triangleId].p_[1].data(),
382 tetIntersections[tetId][triangleId].p_[2].data(), pB.data(), barypB);
383
384 // 2. between the two, find the closest point from the edge
385 // [pivotVertexId, (pivotVertexId+2)%3]
386 // that's the vertex which minimizes its coordinate [(pivotVertexId+1)%3]
387 std::array<double, 3> A = pA, B = pB;
388 std::array<double, 3> baryA = barypA, baryB = barypB;
389 if(fabs(barypB[(pivotVertexId + 1) % 3])
390 < fabs(barypA[(pivotVertexId + 1) % 3])) {
391 // let's switch the two
392 A = pB;
393 B = pA;
394 baryA = barypB;
395 baryB = barypA;
396 }
397
398 bool isAVertex = false;
399 for(int i = 0; i < 3; i++) {
400 if(fabs(baryA[i] - 1) < PREC_DBL) {
401 isAVertex = true;
402 break;
403 }
404 }
405 bool isBVertex = false;
406 for(int i = 0; i < 3; i++) {
407 if(fabs(baryB[i] - 1) < PREC_DBL) {
408 isBVertex = true;
409 break;
410 }
411 }
412 if((isAVertex) && (isBVertex))
413 return -4;
414
415 // 3. create the two new vertices A and B
416 SimplexId vertexIdA = newVertexNumber;
417 newVertexNumber++;
418 tetNewVertices[tetId].resize(tetNewVertices[tetId].size() + 1);
419 for(int i = 0; i < 3; i++)
420 tetNewVertices[tetId].back().p_[i] = A[i];
421 tetNewVertices[tetId].back().polygonEdgeId_ = polygonEdgeId;
422 tetNewVertices[tetId].back().uv_.first
423 = baryA[0] * tetIntersections[tetId][triangleId].uv_[0].first
424 + baryA[1] * tetIntersections[tetId][triangleId].uv_[1].first
425 + baryA[2] * tetIntersections[tetId][triangleId].uv_[2].first;
426 tetNewVertices[tetId].back().uv_.second
427 = baryA[0] * tetIntersections[tetId][triangleId].uv_[0].second
428 + baryA[1] * tetIntersections[tetId][triangleId].uv_[1].second
429 + baryA[2] * tetIntersections[tetId][triangleId].uv_[2].second;
430 tetNewVertices[tetId].back().t_
431 = baryA[0] * tetIntersections[tetId][triangleId].t_[0]
432 + baryA[1] * tetIntersections[tetId][triangleId].t_[1]
433 + baryA[2] * tetIntersections[tetId][triangleId].t_[2];
434 tetNewVertices[tetId].back().isIntersectionPoint_ = true;
435
436 SimplexId vertexIdB = newVertexNumber;
437 newVertexNumber++;
438 tetNewVertices[tetId].resize(tetNewVertices[tetId].size() + 1);
439 for(int i = 0; i < 3; i++)
440 tetNewVertices[tetId].back().p_[i] = B[i];
441 tetNewVertices[tetId].back().polygonEdgeId_ = polygonEdgeId;
442 tetNewVertices[tetId].back().uv_.first
443 = baryB[0] * tetIntersections[tetId][triangleId].uv_[0].first
444 + baryB[1] * tetIntersections[tetId][triangleId].uv_[1].first
445 + baryB[2] * tetIntersections[tetId][triangleId].uv_[2].first;
446 tetNewVertices[tetId].back().uv_.second
447 = baryB[0] * tetIntersections[tetId][triangleId].uv_[0].second
448 + baryB[1] * tetIntersections[tetId][triangleId].uv_[1].second
449 + baryB[2] * tetIntersections[tetId][triangleId].uv_[2].second;
450 tetNewVertices[tetId].back().t_
451 = baryB[0] * tetIntersections[tetId][triangleId].t_[0]
452 + baryB[1] * tetIntersections[tetId][triangleId].t_[1]
453 + baryB[2] * tetIntersections[tetId][triangleId].t_[2];
454 tetNewVertices[tetId].back().isIntersectionPoint_ = true;
455
456 // 4. create the triangle (without saving intersection):
457 // pivotVertexId, A, (pivotVertexId+2)%3
458 createNewIntersectionTriangle(tetId, triangleId, pivotVertexId, -vertexIdA,
459 (pivotVertexId + 2) % 3, tetNewVertices,
460 newTriangleNumber, tetIntersections);
461
462 // 5. create the triangle (saving intersection)
463 // A, B, pivotVertexId+2
464 // special case where a, b and p+2 are actually aligned
465 // we should detect a colinear triangle and test the opposite diagonal instead
466 // a, b, (pivotVertexId+1)%3
468 tetId, triangleId, -vertexIdA, -vertexIdB, (pivotVertexId + 2) % 3,
469 tetNewVertices, newTriangleNumber, tetIntersections, &intersection);
470 if(ret == -1) {
471 // flip the diagonal
473 tetId, triangleId, -vertexIdA, -vertexIdB, (pivotVertexId + 1) % 3,
474 tetNewVertices, newTriangleNumber, tetIntersections, &intersection);
475 createNewIntersectionTriangle(tetId, triangleId, (pivotVertexId + 1) % 3,
476 (pivotVertexId + 2) % 3, -vertexIdA,
477 tetNewVertices, newTriangleNumber,
478 tetIntersections, &intersection);
479 } else {
480 // 6. create the triangle (saving intersection)
481 // pivotVertexId+1, pivotVertexId+2, B
482 createNewIntersectionTriangle(tetId, triangleId, (pivotVertexId + 1) % 3,
483 (pivotVertexId + 2) % 3, -vertexIdB,
484 tetNewVertices, newTriangleNumber,
485 tetIntersections, &intersection);
486 }
487
488 // 7. create the triangle (without saving intersection)
489 // pivotVertexId, B, pivotVertexId+1
490 createNewIntersectionTriangle(tetId, triangleId, pivotVertexId, -vertexIdB,
491 (pivotVertexId + 1) % 3, tetNewVertices,
492 newTriangleNumber, tetIntersections);
493
494 // 8. edit the original triangle to
495 // p, A, B (saving intersection)
496 tetIntersections[tetId][triangleId].intersection_ = intersection;
497 tetIntersections[tetId][triangleId].vertexIds_[0]
498 = tetIntersections[tetId][triangleId].vertexIds_[pivotVertexId];
499 tetIntersections[tetId][triangleId].uv_[0]
500 = tetIntersections[tetId][triangleId].uv_[pivotVertexId];
501 tetIntersections[tetId][triangleId].t_[0]
502 = tetIntersections[tetId][triangleId].t_[pivotVertexId];
503 for(int i = 0; i < 3; i++)
504 tetIntersections[tetId][triangleId].p_[0][i]
505 = tetIntersections[tetId][triangleId].p_[pivotVertexId][i];
506 // vertexA
507 tetIntersections[tetId][triangleId].vertexIds_[1] = -vertexIdA;
508 tetIntersections[tetId][triangleId].uv_[1]
509 = tetNewVertices[tetId][vertexIdA - 1].uv_;
510 tetIntersections[tetId][triangleId].t_[1]
511 = tetNewVertices[tetId][vertexIdA - 1].t_;
512 for(int i = 0; i < 3; i++)
513 tetIntersections[tetId][triangleId].p_[1][i]
514 = tetNewVertices[tetId][vertexIdA - 1].p_[i];
515 // vertexB
516 tetIntersections[tetId][triangleId].vertexIds_[2] = -vertexIdB;
517 tetIntersections[tetId][triangleId].uv_[2]
518 = tetNewVertices[tetId][vertexIdB - 1].uv_;
519 tetIntersections[tetId][triangleId].t_[2]
520 = tetNewVertices[tetId][vertexIdB - 1].t_;
521 for(int i = 0; i < 3; i++)
522 tetIntersections[tetId][triangleId].p_[2][i]
523 = tetNewVertices[tetId][vertexIdB - 1].p_[i];
524
525 return 0;
526}
527
529 const SimplexId &tetId,
530 const SimplexId &triangleId,
531 const SimplexId &vertexId0,
532 const SimplexId &vertexId1,
533 const SimplexId &vertexId2,
534 const vector<vector<Vertex>> &tetNewVertices,
535 SimplexId &newTriangleNumber,
536 vector<vector<IntersectionTriangle>> &tetIntersections,
537 const pair<double, double> *intersection) const {
538
539 if(isIntersectionTriangleColinear(tetId, triangleId, tetIntersections,
540 tetNewVertices, vertexId0, vertexId1,
541 vertexId2)) {
542
543 return -1;
544 }
545
546 tetIntersections[tetId].resize(tetIntersections[tetId].size() + 1);
547 tetIntersections[tetId].back().caseId_
548 = tetIntersections[tetId][triangleId].caseId_;
549 tetIntersections[tetId].back().polygonEdgeId_
550 = tetIntersections[tetId][triangleId].polygonEdgeId_;
551 tetIntersections[tetId].back().triangleId_ = -newTriangleNumber;
552
553 if(intersection) {
554 tetIntersections[tetId].back().intersection_ = *intersection;
555 } else {
556 tetIntersections[tetId].back().intersection_.first = -DBL_MAX;
557 tetIntersections[tetId].back().intersection_.second = -DBL_MAX;
558 }
559
560 newTriangleNumber++;
561
562 // process the vertices
563 for(int i = 0; i < 3; i++) {
564
565 SimplexId vertexId = vertexId0;
566 if(i == 1)
567 vertexId = vertexId1;
568 if(i == 2)
569 vertexId = vertexId2;
570
571 if(vertexId >= 0) {
572 // this is not a newly created vertex
573 tetIntersections[tetId].back().vertexIds_[i]
574 = tetIntersections[tetId][triangleId].vertexIds_[vertexId];
575 tetIntersections[tetId].back().uv_[i]
576 = tetIntersections[tetId][triangleId].uv_[vertexId];
577 tetIntersections[tetId].back().t_[i]
578 = tetIntersections[tetId][triangleId].t_[vertexId];
579 for(int j = 0; j < 3; j++) {
580 tetIntersections[tetId].back().p_[i][j]
581 = tetIntersections[tetId][triangleId].p_[vertexId][j];
582 }
583 } else {
584 // this is a newly created vertex
585 tetIntersections[tetId].back().vertexIds_[i] = vertexId;
586 tetIntersections[tetId].back().uv_[i]
587 = tetNewVertices[tetId][(-vertexId) - 1].uv_;
588 tetIntersections[tetId].back().t_[i]
589 = tetNewVertices[tetId][(-vertexId) - 1].t_;
590 for(int j = 0; j < 3; j++) {
591 tetIntersections[tetId].back().p_[i][j]
592 = tetNewVertices[tetId][(-vertexId) - 1].p_[j];
593 }
594 }
595 }
596
597 return 0;
598}
599
601
602 Timer t;
603
604 vector<vector<pair<SimplexId, SimplexId>>> tetTriangles(tetNumber_);
605 vector<bool> inQueue(tetNumber_, false);
606 vector<SimplexId> tetList;
607
608 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
609 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
610 j++) {
611
612 SimplexId tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
613
614 if(!inQueue[tetId]) {
615 inQueue[tetId] = true;
616 tetList.push_back(tetId);
617 }
618
619 tetTriangles[tetId].push_back(pair<SimplexId, SimplexId>(i, j));
620 }
621 }
622
623#ifdef TTK_ENABLE_OPENMP
624#pragma omp parallel for num_threads(threadNumber_)
625#endif
626 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
627
628 SimplexId tetId = tetList[i];
629
630 if(tetTriangles[tetId].size() >= 2) {
631
632 flipEdges(tetTriangles[tetId]);
633 }
634 }
635
636 this->printMsg(
637 "Performed edge flips", 1.0, t.getElapsedTime(), this->threadNumber_);
638
639 return 0;
640}
641
643 std::vector<std::pair<SimplexId, SimplexId>> &triangles) const {
644
645 for(SimplexId it = 0; it < (SimplexId)triangles.size(); it++) {
646
647 if(it == 2)
648 break;
649
650 // avoid infinite loops
651 bool hasFlipped = false;
652
653 // go in a greedy manner
654 // for a triangle, grab its smallest angle alpha.
655 // look at all its neighbors.
656 // do the edge flip virtually and evaluate the smallest angle beta.
657 // if beta is bigger than alpha do the edge flip
658
659 // heuristic: sort the triangles in order of their minimum angle
660 // (make sure we start with bad angles)
661 vector<pair<double, pair<SimplexId, SimplexId>>> localTriangles;
662 for(SimplexId i = 0; i < (SimplexId)triangles.size(); i++) {
663 std::array<SimplexId, 3> vertexIds{};
664
665 vertexIds[0]
666 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
667 .vertexIds_[0];
668 vertexIds[1]
669 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
670 .vertexIds_[1];
671 vertexIds[2]
672 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
673 .vertexIds_[2];
674
675 std::array<double, 3> angles{};
677 (*globalVertexList_)[vertexIds[0]].p_.data(),
678 (*globalVertexList_)[vertexIds[1]].p_.data(),
679 (*globalVertexList_)[vertexIds[2]].p_.data(), angles);
680
681 double alpha = -1;
682 for(int j = 0; j < 3; j++) {
683 if((alpha < 0) || (fabs(angles[j]) < alpha))
684 alpha = fabs(angles[j]);
685 }
686
687 localTriangles.emplace_back(alpha, triangles[i]);
688 }
689
690 const auto FiberSurfaceTriangleCmp
691 = [](const pair<double, pair<SimplexId, SimplexId>> &t0,
692 const pair<double, pair<SimplexId, SimplexId>> &t1) {
693 return t0.first < t1.first;
694 };
695
696 sort(localTriangles.begin(), localTriangles.end(), FiberSurfaceTriangleCmp);
697
698 for(SimplexId i = 0; i < (SimplexId)triangles.size(); i++) {
699 triangles[i] = localTriangles[i].second;
700 }
701
702 for(SimplexId i = 0; i < (SimplexId)triangles.size(); i++) {
703
704 std::array<SimplexId, 3> vertexIds{};
705
706 vertexIds[0]
707 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
708 .vertexIds_[0];
709 vertexIds[1]
710 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
711 .vertexIds_[1];
712 vertexIds[2]
713 = (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]
714 .vertexIds_[2];
715
716 std::array<double, 3> angles{};
718 (*globalVertexList_)[vertexIds[0]].p_.data(),
719 (*globalVertexList_)[vertexIds[1]].p_.data(),
720 (*globalVertexList_)[vertexIds[2]].p_.data(), angles);
721
722 double alpha = -1;
723 for(int j = 0; j < 3; j++) {
724 if((alpha < 0) || (fabs(angles[j]) < alpha))
725 alpha = fabs(angles[j]);
726 }
727
728 // look at neighbors
729 for(SimplexId j = 0; j < (SimplexId)triangles.size(); j++) {
730 if((i != j)
731 // same polygon edge
732 && (((*polygonEdgeTriangleLists_[triangles[j].first])[triangles[j]
733 .second]
734 .polygonEdgeId_
735 == (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i]
736 .second]
737 .polygonEdgeId_))) {
738
739 SimplexId commonVertexId0 = -1;
740 SimplexId commonVertexId1 = -1;
741 SimplexId otherNonCommonVertexId = -1;
742
743 for(int k = 0; k < 3; k++) {
744 SimplexId vertexId
745 = (*polygonEdgeTriangleLists_[triangles[j].first])[triangles[j]
746 .second]
747 .vertexIds_[k];
748
749 bool isCommon = false;
750
751 for(int l = 0; l < 3; l++) {
752 if(vertexId == vertexIds[l]) {
753 if(commonVertexId0 == -1) {
754 commonVertexId0 = vertexId;
755 } else {
756 commonVertexId1 = vertexId;
757 }
758 isCommon = true;
759 }
760 }
761 if(!isCommon)
762 otherNonCommonVertexId = vertexId;
763 }
764
765 if((commonVertexId0 != -1) && (commonVertexId1 != -1)) {
766
767 if((!(*globalVertexList_)[commonVertexId0].isIntersectionPoint_)
768 && (!(*globalVertexList_)[commonVertexId1]
769 .isIntersectionPoint_)) {
770
771 SimplexId nonCommonVertexId = -1;
772 for(int k = 0; k < 3; k++) {
773 if((vertexIds[k] != commonVertexId0)
774 && (vertexIds[k] != commonVertexId1)) {
775 nonCommonVertexId = vertexIds[k];
776 break;
777 }
778 }
779
780 // we have a neighbor
781 // now we want to evaluate the angles of the following triangles
782 // nonCommonVertexId, commonVertexId0, otherNonCommonVertexId
783 // nonCommonVertexId, commonVertexId1, otherNonCommonVertexId
784 // if the min angles of these two guys is bigger than alpha, flip!
785
786 if((nonCommonVertexId != -1) && (otherNonCommonVertexId != -1)) {
787
788 std::array<double, 3> beta0angles{}, beta1angles{};
789
791 (*globalVertexList_)[nonCommonVertexId].p_.data(),
792 (*globalVertexList_)[commonVertexId0].p_.data(),
793 (*globalVertexList_)[otherNonCommonVertexId].p_.data(),
794 beta0angles);
795
796 double beta0 = -1;
797 for(int k = 0; k < 3; k++) {
798 if((beta0 < 0) || (fabs(beta0angles[j]) < beta0))
799 beta0 = fabs(beta0angles[j]);
800 }
801
803 (*globalVertexList_)[nonCommonVertexId].p_.data(),
804 (*globalVertexList_)[commonVertexId1].p_.data(),
805 (*globalVertexList_)[otherNonCommonVertexId].p_.data(),
806 beta1angles);
807
808 double beta1 = -1;
809 for(int k = 0; k < 3; k++) {
810 if((beta1 < 0) || (fabs(beta1angles[j]) < beta1))
811 beta1 = fabs(beta1angles[j]);
812 }
813
814 if((beta0 > alpha) && (beta1 > alpha)) {
815 // flip!
816
817 if(isEdgeFlippable(commonVertexId0, commonVertexId1,
818 nonCommonVertexId,
819 otherNonCommonVertexId)) {
820
821 // original triangle:
822 // nonCommonVertexId, otherNonCommonVertexId,
823 // commonVertexId0
824 (*polygonEdgeTriangleLists_[triangles[i]
825 .first])[triangles[i].second]
826 .vertexIds_[0]
827 = nonCommonVertexId;
828 (*polygonEdgeTriangleLists_[triangles[i]
829 .first])[triangles[i].second]
830 .vertexIds_[1]
831 = otherNonCommonVertexId;
832 (*polygonEdgeTriangleLists_[triangles[i]
833 .first])[triangles[i].second]
834 .vertexIds_[2]
835 = commonVertexId0;
836
837 // other triangle:
838 // nonCommonVertexId, otherNonCommonVertexId,
839 // commonVertexId1
840 (*polygonEdgeTriangleLists_[triangles[j]
841 .first])[triangles[j].second]
842 .vertexIds_[0]
843 = nonCommonVertexId;
844 (*polygonEdgeTriangleLists_[triangles[j]
845 .first])[triangles[j].second]
846 .vertexIds_[1]
847 = otherNonCommonVertexId;
848 (*polygonEdgeTriangleLists_[triangles[j]
849 .first])[triangles[j].second]
850 .vertexIds_[2]
851 = commonVertexId1;
852
853 hasFlipped = true;
854 }
855 }
856 }
857 }
858 }
859 }
860 if(hasFlipped)
861 break;
862 }
863 if(hasFlipped)
864 break;
865 }
866
867 if(!hasFlipped)
868 break;
869 }
870
871 return 0;
872}
873
875 const SimplexId &tetId,
876 const SimplexId &triangleId,
877 const vector<vector<IntersectionTriangle>> &tetIntersections,
878 pair<double, double> &extremity0,
879 pair<double, double> &extremity1) const {
880
881 std::array<double, 2> p0{}, p1{}, p{};
882 std::array<double, 2> baryCentrics{};
883 bool isInBetween = true;
884
885 // check for edges that project to points first
886 for(int i = 0; i < 3; i++) {
887
888 p[0] = tetIntersections[tetId][triangleId].uv_[i].first;
889 p[1] = tetIntersections[tetId][triangleId].uv_[i].second;
890
891 p0[0] = tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first;
892 p0[1] = tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second;
893
894 p1[0] = tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].first;
895 p1[1] = tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].second;
896
897 if((fabs(p0[0] - p1[0]) < PREC_FLT) && (fabs(p0[1] - p1[1]) < PREC_FLT)) {
898 // one edge of the triangle projects to a point
899 extremity0.first = p[0];
900 extremity0.second = p[1];
901
902 extremity1.first = p0[0];
903 extremity1.second = p0[1];
904
905 return 0;
906 }
907 }
908
909 for(int i = 0; i < 3; i++) {
910
911 p[0] = tetIntersections[tetId][triangleId].uv_[i].first;
912 p[1] = tetIntersections[tetId][triangleId].uv_[i].second;
913
914 p0[0] = tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].first;
915 p0[1] = tetIntersections[tetId][triangleId].uv_[(i + 1) % 3].second;
916
917 p1[0] = tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].first;
918 p1[1] = tetIntersections[tetId][triangleId].uv_[(i + 2) % 3].second;
919
921 p0.data(), p1.data(), p.data(), baryCentrics, 2);
922
923 isInBetween = true;
924 for(int j = 0; j < 2; j++) {
925
926 if((baryCentrics[j] < -PREC_FLT) || (baryCentrics[j] > 1 + PREC_FLT)) {
927 isInBetween = false;
928 break;
929 }
930 }
931 if(isInBetween) {
932 extremity0.first = p0[0];
933 extremity0.second = p0[1];
934
935 extremity1.first = p1[0];
936 extremity1.second = p1[1];
937
938 return 0;
939 }
940 }
941
942 return 0;
943}
944
946 const double *p1,
947 const double *p2) const {
948
949 if((p0[0] == p1[0]) && (p0[1] == p1[1]) && (p0[2] == p1[2]))
950 return true;
951
952 if((p2[0] == p1[0]) && (p2[1] == p1[1]) && (p2[2] == p1[2]))
953 return true;
954
955 if((p0[0] == p2[0]) && (p0[1] == p2[1]) && (p0[2] == p2[2]))
956 return true;
957
958 return false;
959}
960
961int FiberSurface::interpolateBasePoints(const std::array<double, 3> &p0,
962 const pair<double, double> &uv0,
963 const double &t0,
964 const std::array<double, 3> &p1,
965 const pair<double, double> &uv1,
966 const double &t1,
967 const double &t,
968 Vertex &v) const {
969
970 for(int j = 0; j < 3; j++) {
971 v.p_[j] = p0[j] + ((t - t0) / (t1 - t0)) * (p1[j] - p0[j]);
972 }
973
974 // interpolate uv
975 v.uv_.first = uv0.first + ((t - t0) / (t1 - t0)) * (uv1.first - uv0.first);
976 v.uv_.second
977 = uv0.second + ((t - t0) / (t1 - t0)) * (uv1.second - uv0.second);
978
979 v.isBasePoint_ = false;
980
981 return 0;
982}
983
985 const SimplexId &source,
986 const SimplexId &destination,
987 const SimplexId &pivotVertexId,
988 const vector<pair<SimplexId, SimplexId>> &starNeighbors) const {
989
990 // NOTE:
991 // here I should really look at the triangles' normals more than the angles...
992
993 SimplexId baseId = -1;
994 double baseAngle = 0;
995 for(SimplexId i = 0; i < (SimplexId)starNeighbors.size(); i++) {
996 if(((starNeighbors[i].first == source)
997 && (starNeighbors[i].second == destination))
998 || ((starNeighbors[i].first == destination)
999 && (starNeighbors[i].second == source))) {
1000
1001 baseAngle = Geometry::angle((*globalVertexList_)[source].p_.data(),
1002 (*globalVertexList_)[pivotVertexId].p_.data(),
1003 (*globalVertexList_)[pivotVertexId].p_.data(),
1004 (*globalVertexList_)[destination].p_.data());
1005 baseId = i;
1006 break;
1007 }
1008 }
1009
1010 for(SimplexId i = 0; i < (SimplexId)starNeighbors.size(); i++) {
1011 if(i != baseId) {
1012 if((starNeighbors[i].first == source)
1013 || (starNeighbors[i].first == destination)
1014 || (starNeighbors[i].second == source)
1015 || (starNeighbors[i].second == destination)) {
1016
1017 double localAngle = Geometry::angle(
1018 (*globalVertexList_)[starNeighbors[i].first].p_.data(),
1019 (*globalVertexList_)[pivotVertexId].p_.data(),
1020 (*globalVertexList_)[pivotVertexId].p_.data(),
1021 (*globalVertexList_)[starNeighbors[i].second].p_.data());
1022 if(localAngle + baseAngle > 0.9 * M_PI)
1023 return false;
1024 }
1025 }
1026 }
1027 return true;
1028}
1029
1031 const SimplexId &edgeVertexId1,
1032 const SimplexId &otherVertexId0,
1033 const SimplexId &otherVertexId1) const {
1034
1035 double angle0
1036 = Geometry::angle((*globalVertexList_)[edgeVertexId0].p_.data(),
1037 (*globalVertexList_)[edgeVertexId1].p_.data(),
1038 (*globalVertexList_)[edgeVertexId1].p_.data(),
1039 (*globalVertexList_)[otherVertexId0].p_.data());
1040
1041 double angle1
1042 = Geometry::angle((*globalVertexList_)[edgeVertexId0].p_.data(),
1043 (*globalVertexList_)[edgeVertexId1].p_.data(),
1044 (*globalVertexList_)[edgeVertexId1].p_.data(),
1045 (*globalVertexList_)[otherVertexId1].p_.data());
1046
1047 if(angle0 + angle1 > 0.9 * M_PI)
1048 return false;
1049
1050 // now do the angles at the other extremity of the edge.
1051 angle0 = Geometry::angle((*globalVertexList_)[edgeVertexId1].p_.data(),
1052 (*globalVertexList_)[edgeVertexId0].p_.data(),
1053 (*globalVertexList_)[edgeVertexId0].p_.data(),
1054 (*globalVertexList_)[otherVertexId0].p_.data());
1055
1056 angle1 = Geometry::angle((*globalVertexList_)[edgeVertexId1].p_.data(),
1057 (*globalVertexList_)[edgeVertexId0].p_.data(),
1058 (*globalVertexList_)[edgeVertexId0].p_.data(),
1059 (*globalVertexList_)[otherVertexId1].p_.data());
1060
1061 if(angle0 + angle1 > 0.9 * M_PI)
1062 return false;
1063
1064 return true;
1065}
1066
1067int FiberSurface::mergeEdges(const double &distanceThreshold) const {
1068
1069 Timer t;
1070
1071 // TODO: forbid edge collapse for edges that do not live on the boundary of
1072 // tetrahedrons....
1073 // this is the case when the tetId of the two triangles is the same.
1074 // THAT is not OK
1075
1076 SimplexId initVertexNumber = (*globalVertexList_).size();
1077
1078 for(SimplexId it = 0; it < (SimplexId)(*globalVertexList_).size(); it++) {
1079 // avoid infinite loops
1080
1081 // make the local list of vertex neighbors
1082 vector<vector<SimplexId>> vertexNeighbors((*globalVertexList_).size());
1083 vector<vector<SimplexId>> vertexNeighborsTets((*globalVertexList_).size());
1084 // for each triangle of the star, list of other two vertices
1085 vector<vector<pair<SimplexId, SimplexId>>> vertexTriangleNeighbors(
1086 vertexNeighbors.size());
1087
1088 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1089 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1090 j++) {
1091
1092 for(int k = 0; k < 3; k++) {
1093 SimplexId vertexId = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
1094
1095 vertexTriangleNeighbors[vertexId].resize(
1096 vertexTriangleNeighbors[vertexId].size() + 1);
1097 vertexTriangleNeighbors[vertexId].back().first = -1;
1098 vertexTriangleNeighbors[vertexId].back().second = -1;
1099 vertexNeighborsTets[vertexId].push_back(
1100 (*polygonEdgeTriangleLists_[i])[j].tetId_);
1101
1102 for(int l = 0; l < 3; l++) {
1103 if(l != k) {
1104 SimplexId otherVertexId
1105 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[l];
1106
1107 if(vertexTriangleNeighbors[vertexId].back().first == -1) {
1108 vertexTriangleNeighbors[vertexId].back().first = otherVertexId;
1109 } else {
1110 vertexTriangleNeighbors[vertexId].back().second = otherVertexId;
1111 }
1112
1113 bool isIn = false;
1114 for(SimplexId m = 0;
1115 m < (SimplexId)vertexNeighbors[vertexId].size(); m++) {
1116 if(vertexNeighbors[vertexId][m] == otherVertexId) {
1117 isIn = true;
1118 break;
1119 }
1120 }
1121 if(!isIn) {
1122 vertexNeighbors[vertexId].push_back(otherVertexId);
1123 }
1124 }
1125 }
1126 }
1127 }
1128 }
1129
1130 bool hasMerged = false;
1131
1132#ifdef TTK_ENABLE_OPENMP
1133#pragma omp parallel for num_threads(threadNumber_)
1134#endif
1135 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1136
1137 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1138 j++) {
1139
1140 SimplexId minimizer = -1;
1141 double minDistance = -1;
1142
1143 // find the smallest edge on this triangle
1144 for(int k = 0; k < 3; k++) {
1145
1146 SimplexId vertexId0
1147 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
1148 SimplexId vertexId1
1149 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[(k + 1) % 3];
1150
1151 bool areAlreadySnapped = true;
1152 for(int l = 0; l < 3; l++) {
1153 if((*globalVertexList_)[vertexId0].p_[l]
1154 != (*globalVertexList_)[vertexId1].p_[l]) {
1155
1156 areAlreadySnapped = false;
1157 break;
1158 }
1159 }
1160
1161 if(((*globalVertexList_)[vertexId0].isBasePoint_)
1162 && ((*globalVertexList_)[vertexId1].isBasePoint_)) {
1163 areAlreadySnapped = true;
1164 }
1165
1166 if(!areAlreadySnapped) {
1167 double distance
1168 = Geometry::distance((*globalVertexList_)[vertexId0].p_.data(),
1169 (*globalVertexList_)[vertexId1].p_.data());
1170
1171 if((minDistance == -1) || (distance < minDistance)) {
1172 minDistance = distance;
1173 minimizer = k;
1174 }
1175 }
1176 }
1177
1178 if((minDistance != -1) && (minDistance < distanceThreshold)) {
1179 SimplexId vertexId0
1180 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[minimizer];
1181 SimplexId vertexId1 = (*polygonEdgeTriangleLists_[i])[j]
1182 .vertexIds_[(minimizer + 1) % 3];
1183
1184 // find the number of common neighbors
1185 vector<SimplexId> commonNeighbors;
1186 for(SimplexId k = 0; k < (SimplexId)vertexNeighbors[vertexId0].size();
1187 k++) {
1188 for(SimplexId l = 0;
1189 l < (SimplexId)vertexNeighbors[vertexId1].size(); l++) {
1190 if(vertexNeighbors[vertexId0][k]
1191 == vertexNeighbors[vertexId1][l]) {
1192 commonNeighbors.push_back(vertexNeighbors[vertexId0][k]);
1193 }
1194 }
1195 }
1196
1197 if(commonNeighbors.size() == 2) {
1198
1199 // we may create extra non-manifold triangles otherwise
1200 // plus we don't collapse edges that are strictly inside a tet
1201 // (only collsapse those on the boundary between two tets)
1202
1203 SimplexId source = vertexId0, destination = vertexId1;
1204
1205 if((*globalVertexList_)[destination].isBasePoint_) {
1206 source = vertexId1;
1207 destination = vertexId0;
1208 }
1209
1210 if(!(*globalVertexList_)[destination].isBasePoint_) {
1211
1212 // now check the angles
1213 bool isCollapsible = true;
1214 for(SimplexId k = 0; k < (SimplexId)commonNeighbors.size(); k++) {
1216 source, destination, commonNeighbors[k],
1217 vertexTriangleNeighbors[commonNeighbors[k]])) {
1218 isCollapsible = false;
1219 break;
1220 }
1221 }
1222
1223 // different tets?
1224 SimplexId tetId0 = -1, tetId1 = -1;
1225 for(SimplexId k = 0;
1226 k < (SimplexId)vertexNeighborsTets[source].size(); k++) {
1227 if((vertexTriangleNeighbors[source][k].first == destination)
1228 || (vertexTriangleNeighbors[source][k].second
1229 == destination)) {
1230 if(tetId0 == -1) {
1231 tetId0 = vertexNeighborsTets[source][k];
1232 } else {
1233 tetId1 = vertexNeighborsTets[source][k];
1234 }
1235 }
1236 }
1237
1238 if(tetId0 == tetId1) {
1239 isCollapsible = false;
1240 }
1241
1242#ifdef TTK_ENABLE_OPENMP
1243#pragma omp critical
1244#endif
1245 if(isCollapsible) {
1246 for(int k = 0; k < 3; k++) {
1247 (*globalVertexList_)[destination].p_[k]
1248 = (*globalVertexList_)[source].p_[k];
1249 }
1250
1251 hasMerged = true;
1252 }
1253 }
1254 }
1255 }
1256 }
1257 }
1258
1259 if(!hasMerged)
1260 break;
1261
1262 // now update the vertices
1263 mergeVertices(0);
1264 }
1265
1266 this->printMsg(
1267 "Performed edge collapses", 1.0, t.getElapsedTime(), this->threadNumber_);
1268 this->printMsg(std::vector<std::vector<std::string>>{
1269 {"#Vertices removed",
1270 std::to_string(initVertexNumber - (*globalVertexList_).size())}});
1271
1272 return 0;
1273}
1274
1275int FiberSurface::mergeVertices(const double &distanceThreshold) const {
1276
1277 Timer t;
1278
1279 vector<Vertex> tmpList = (*globalVertexList_);
1280
1281 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1282 tmpList[i].localId_ = i;
1283 }
1284
1285 const auto FiberSurfaceVertexComparisonX
1286 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1287 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1288 // let's consider x coordinates are equal
1289 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1290 // let's consider y coordinates are equal
1291 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1292 // let's consider z coordinates are equal
1293 // NOTE: the local Id should be sufficient
1294 return v0.globalId_ < v1.globalId_;
1295 } else
1296 return v0.p_[2] < v1.p_[2];
1297 } else
1298 return v0.p_[1] < v1.p_[1];
1299 } else {
1300 return v0.p_[0] < v1.p_[0];
1301 }
1302 };
1303
1304 const auto FiberSurfaceVertexComparisonY
1305 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1306 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1307 // let's consider y coordinates are equal
1308 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1309 // let's consider z coordinates are equal
1310 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1311 // let's consider x coordinates are equal
1312 // NOTE: the local Id should be sufficient
1313 return v0.globalId_ < v1.globalId_;
1314 } else
1315 return v0.p_[0] < v1.p_[0];
1316 } else
1317 return v0.p_[2] < v1.p_[2];
1318 } else {
1319 return v0.p_[1] < v1.p_[1];
1320 }
1321 };
1322
1323 const auto FiberSurfaceVertexComparisonZ
1324 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1325 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1326 // let's consider z coordinates are equal
1327 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1328 // let's consider x coordinates are equal
1329 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1330 // let's consider y coordinates are equal
1331 // NOTE: the local Id should be sufficient
1332 return v0.globalId_ < v1.globalId_;
1333 } else
1334 return v0.p_[1] < v1.p_[1];
1335 } else
1336 return v0.p_[0] < v1.p_[0];
1337 } else {
1338 return v0.p_[2] < v1.p_[2];
1339 }
1340 };
1341
1342 // now do a parallel sort
1343 SimplexId uniqueVertexNumber = 0;
1344 for(int k = 0; k < 3; k++) {
1345
1346 switch(k) {
1347 case 0:
1348 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonX);
1349 break;
1350
1351 case 1:
1352 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonY);
1353 break;
1354
1355 case 2:
1356 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonZ);
1357 break;
1358 }
1359
1360 // NOTE:
1361 // ethaneDiol dataset
1362 // 0.019224 in parallel (2 cores HT)
1363 // 0.041353 in sequential (speedup: x2.15, perfect (HT))
1364
1365 // now merge the thing
1366 // 1. Identity duplicates
1367 // NOTE: order is important, so no parallelism
1368 // NOTE: points are represented with single precision (floats) so use that
1369 // as a limit for distances
1370 uniqueVertexNumber = 0;
1371 double distance = 0;
1372 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1373
1374 bool canMerge = false;
1375
1376 if(i) {
1377 distance
1378 = Geometry::distance(tmpList[i].p_.data(), tmpList[i - 1].p_.data());
1379
1380 if(distance <= distanceThreshold) {
1381
1382 // one of the two is -1 (interior vertex)
1383 // if((tmpList[i - 1].meshEdge_.first == -1)
1384 // ||(tmpList[i].meshEdge_.first == -1)){
1385 // canMerge = true;
1386 // }
1387 //
1388 // // one vertex in common
1389 // if((tmpList[i].meshEdge_.first ==
1390 // tmpList[i - 1].meshEdge_.first)
1391 // ||
1392 // (tmpList[i].meshEdge_.first ==
1393 // tmpList[i - 1].meshEdge_.second)){
1394 // canMerge = true;
1395 // }
1396 //
1397 // // one vertex in common
1398 // if((tmpList[i].meshEdge_.second ==
1399 // tmpList[i - 1].meshEdge_.first)
1400 // ||
1401 // (tmpList[i].meshEdge_.second ==
1402 // tmpList[i - 1].meshEdge_.second)){
1403 // canMerge = true;
1404 // }
1405
1406 // NOTE: still some bugs here in terms of manifoldness.
1407
1408 if((tmpList[i].meshEdge_.first == tmpList[i - 1].meshEdge_.first)
1409 && (tmpList[i].meshEdge_.second
1410 == tmpList[i - 1].meshEdge_.second)) {
1411 canMerge = true;
1412 }
1413 if((tmpList[i].meshEdge_.first == tmpList[i - 1].meshEdge_.second)
1414 && (tmpList[i].meshEdge_.second
1415 == tmpList[i - 1].meshEdge_.first)) {
1416 canMerge = true;
1417 }
1418
1419 // canMerge = true;
1420 }
1421
1422 if(canMerge) {
1423 tmpList[i].globalId_ = tmpList[i - 1].globalId_;
1424 if((tmpList[i].isBasePoint_) || (tmpList[i - 1].isBasePoint_)) {
1425 tmpList[i].isBasePoint_ = true;
1426 tmpList[i - 1].isBasePoint_ = true;
1427 }
1428 if((tmpList[i].isIntersectionPoint_)
1429 || (tmpList[i - 1].isIntersectionPoint_)) {
1430 tmpList[i].isIntersectionPoint_ = true;
1431 tmpList[i - 1].isIntersectionPoint_ = true;
1432 }
1433 for(int j = 0; j < 3; j++) {
1434 tmpList[i].p_[j] = tmpList[i - 1].p_[j];
1435 }
1436
1437 // update meshEdge...
1438 if((tmpList[i - 1].meshEdge_.first == -1)
1439 && (tmpList[i].meshEdge_.first != -1)) {
1440 tmpList[i - 1].meshEdge_ = tmpList[i].meshEdge_;
1441 }
1442 if((tmpList[i].meshEdge_.first == -1)
1443 && (tmpList[i - 1].meshEdge_.first != -1)) {
1444 tmpList[i].meshEdge_ = tmpList[i - 1].meshEdge_;
1445 }
1446 }
1447 }
1448
1449 if((!i) || (!canMerge)) {
1450 tmpList[i].globalId_ = uniqueVertexNumber;
1451 uniqueVertexNumber++;
1452 }
1453
1454 (*globalVertexList_)[tmpList[i].localId_].globalId_
1455 = tmpList[i].globalId_;
1456 }
1457 }
1458
1459 // update the triangles
1460 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1461 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1462 j++) {
1463 for(int k = 0; k < 3; k++) {
1464 (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k]
1466 .vertexIds_[k]]
1467 .globalId_;
1468 }
1469 }
1470 }
1471
1472 // 2. create the actual global list
1473 // NOTE: order is no longer important, we can do this in parallel.
1474 (*globalVertexList_).resize(uniqueVertexNumber);
1475
1476 // duplicate vertices but only a subset have the right edge information
1477 // (because they have actually been computed on the edge)
1478 SimplexId lastId = -1;
1479 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1480
1481 if(lastId != -1) {
1482 if((tmpList[lastId].globalId_ == tmpList[i].globalId_)
1483 && (tmpList[i].meshEdge_.first != -1)) {
1484
1485 (*globalVertexList_)[tmpList[lastId].globalId_].meshEdge_
1486 = tmpList[i].meshEdge_;
1487 (*globalVertexList_)[tmpList[lastId].globalId_].isBasePoint_ = true;
1488 lastId = -1;
1489 }
1490 }
1491
1492 if((!i) || (tmpList[i].globalId_ != tmpList[i - 1].globalId_)) {
1493 (*globalVertexList_)[tmpList[i].globalId_] = tmpList[i];
1494 if((*globalVertexList_)[tmpList[i].globalId_].meshEdge_.first == -1) {
1495 lastId = i;
1496 } else {
1497 lastId = -1;
1498 }
1499 }
1500 }
1501
1502 // 3. update the 2-sheets, ignore zero-area triangles
1503 // and free their vertex memory
1504 // NOTE: order is no longer important, parallel
1505 vector<vector<bool>> keepTriangle(polygonEdgeTriangleLists_.size());
1506 vector<vector<FiberSurface::Triangle>> tmpTriangleLists(
1508 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1509 tmpTriangleLists[i] = (*polygonEdgeTriangleLists_[i]);
1510 keepTriangle[i].resize((*polygonEdgeTriangleLists_[i]).size(), true);
1511 }
1512
1513#ifdef TTK_ENABLE_OPENMP
1514#pragma omp parallel for num_threads(threadNumber_)
1515#endif
1516 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1517 for(SimplexId j = 0; j < (SimplexId)tmpTriangleLists[i].size(); j++) {
1518 for(int k = 0; k < 3; k++) {
1519 if((k)
1520 && (tmpTriangleLists[i][j].vertexIds_[k]
1521 == tmpTriangleLists[i][j].vertexIds_[(k - 1)])) {
1522 keepTriangle[i][j] = false;
1523 }
1524 }
1525 if(tmpTriangleLists[i][j].vertexIds_[0]
1526 == tmpTriangleLists[i][j].vertexIds_[2]) {
1527 keepTriangle[i][j] = false;
1528 }
1529 }
1530
1531 // now copy triangles over with non zero-area triangles
1532 // NOTE: no need to re-allocate the memory, we know we are not going to use
1533 // more.
1534 (*polygonEdgeTriangleLists_[i]).clear();
1535 for(SimplexId j = 0; j < (SimplexId)tmpTriangleLists[i].size(); j++) {
1536
1537 if(keepTriangle[i][j]) {
1538 (*polygonEdgeTriangleLists_[i]).push_back(tmpTriangleLists[i][j]);
1539 }
1540 }
1541 }
1542
1543 this->printMsg(
1544 "Output made manifold", 1.0, t.getElapsedTime(), this->threadNumber_);
1545
1546 return 0;
1547}
1548
1549int FiberSurface::snapToBasePoint(const vector<vector<double>> &basePoints,
1550 const vector<pair<double, double>> &uv,
1551 const vector<double> &t,
1552 Vertex &v) const {
1553
1555 return -1;
1556
1557 SimplexId minimizer = 0;
1558 double minDistance = -1;
1559
1560 for(SimplexId i = 0; i < (SimplexId)basePoints.size(); i++) {
1561 double distance = Geometry::distance(basePoints[i].data(), v.p_.data());
1562 if((minDistance < 0) || (distance < minDistance)) {
1563 minDistance = distance;
1564 minimizer = i;
1565 }
1566 }
1567
1568 if(minDistance < pointSnappingThreshold_) {
1569 for(int i = 0; i < 3; i++) {
1570 v.p_[i] = basePoints[minimizer][i];
1571 }
1572 v.uv_ = uv[minimizer];
1573 v.t_ = t[minimizer];
1574 v.isBasePoint_ = true;
1575 }
1576
1577 return 0;
1578}
1579
1581
1582 vector<bool> inQueue(tetNumber_, false);
1583 vector<vector<pair<SimplexId, SimplexId>>> tetTriangles(tetNumber_);
1584 vector<SimplexId> tetList;
1585
1586 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1587
1588 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1589 j++) {
1590
1591 SimplexId tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
1592
1593 if(!inQueue[tetId]) {
1594 inQueue[tetId] = true;
1595 tetList.push_back(tetId);
1596 }
1597
1598 tetTriangles[tetId].push_back(pair<SimplexId, SimplexId>(i, j));
1599 }
1600 }
1601
1602#ifdef TTK_ENABLE_OPENMP
1603#pragma omp parallel for num_threads(threadNumber_)
1604#endif
1605 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
1606 snapVertexBarycentrics(tetList[i], tetTriangles[tetList[i]]);
1607 }
1608
1609 mergeVertices(0);
1610
1611 return 0;
1612}
1613
1615 const SimplexId &tetId,
1616 const std::vector<std::pair<SimplexId, SimplexId>> &triangles) const {
1617
1618 for(SimplexId i = 0; i < (SimplexId)triangles.size(); i++) {
1619
1620 Triangle *t = &(
1621 (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]);
1622
1623 for(int j = 0; j < 3; j++) {
1624 SimplexId vertexId = t->vertexIds_[j];
1625
1626 // check for each triangle of the tet
1627 double minimum = -DBL_MAX;
1628 std::array<double, 3> minBarycentrics{};
1629 std::array<SimplexId, 3> minimizer{};
1630
1631 for(int k = 0; k < 2; k++) {
1632 for(int l = k + 1; l < 3; l++) {
1633 for(int m = l + 1; m < 4; m++) {
1634
1635 SimplexId vertexId0 = tetList_[5 * tetId + 1 + k];
1636 SimplexId vertexId1 = tetList_[5 * tetId + 1 + l];
1637 SimplexId vertexId2 = tetList_[5 * tetId + 1 + m];
1638
1639 std::array<double, 3> p0{}, p1{}, p2{};
1640
1641 for(int n = 0; n < 3; n++) {
1642 p0[n] = pointSet_[3 * vertexId0 + n];
1643 p1[n] = pointSet_[3 * vertexId1 + n];
1644 p2[n] = pointSet_[3 * vertexId2 + n];
1645 }
1646
1647 std::array<double, 3> barycentrics{};
1649 p0.data(), p1.data(), p2.data(),
1650 (*globalVertexList_)[vertexId].p_.data(), barycentrics);
1651
1652 if((barycentrics[0] != -1.0) && (barycentrics[1] != -1.0)
1653 && (barycentrics[2] != -1.0)) {
1654 // vertexId lies in the current triangle
1655
1656 double localMin = -1.0;
1657 for(int n = 0; n < 3; n++) {
1658 if((localMin == -1.0) || (fabs(barycentrics[n]) < localMin)) {
1659 localMin = fabs(barycentrics[n]);
1660 }
1661 }
1662
1663 if((minimum == -DBL_MAX) || (localMin < minimum)) {
1664 minimum = localMin;
1665 minBarycentrics = barycentrics;
1666 minimizer[0] = vertexId0;
1667 minimizer[1] = vertexId1;
1668 minimizer[2] = vertexId2;
1669 }
1670 }
1671 }
1672 }
1673 }
1674
1675 if((minimum != -DBL_MAX) && (minimum < PREC_FLT_2)) {
1676 double sum = 0;
1677 int numberOfZeros = 0;
1678 for(int k = 0; k < 3; k++) {
1679 if(minBarycentrics[k] < PREC_FLT_2) {
1680 minBarycentrics[k] = 0;
1681 numberOfZeros++;
1682 }
1683 sum += minBarycentrics[k];
1684 }
1685
1686 sum = (1 - sum) / numberOfZeros;
1687
1688 for(int k = 0; k < 3; k++) {
1689 if(minBarycentrics[k] >= PREC_FLT_2) {
1690 minBarycentrics[k] += sum;
1691 }
1692 }
1693
1694 // now do the re-location
1695#ifdef TTK_ENABLE_OPENMP
1696#pragma omp critical
1697 for(int k = 0; k < 3; k++) {
1698 (*globalVertexList_)[vertexId].p_[k]
1699 = minBarycentrics[0] * ((double)pointSet_[3 * minimizer[0] + k])
1700 + minBarycentrics[1] * ((double)pointSet_[3 * minimizer[1] + k])
1701 + minBarycentrics[2] * ((double)pointSet_[3 * minimizer[2] + k]);
1702 }
1703#endif
1704 }
1705 }
1706 }
1707
1708 return 0;
1709}
#define M_PI
Definition: Os.h:50
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
int interpolateBasePoints(const std::array< double, 3 > &p0, const std::pair< double, double > &uv0, const double &t0, const std::array< double, 3 > &p1, const std::pair< double, double > &uv1, const double &t1, const double &t, Vertex &v) const
double pointSnappingThreshold_
Definition: FiberSurface.h:491
std::vector< std::vector< Triangle > * > polygonEdgeTriangleLists_
Definition: FiberSurface.h:498
int getNumberOfCommonVertices(const SimplexId &tetId, const SimplexId &triangleId0, const SimplexId &triangleId1, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections) const
int mergeEdges(const double &distanceThreshold) const
int getTriangleRangeExtremities(const SimplexId &tetId, const SimplexId &triangleId, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::pair< double, double > &extremity0, std::pair< double, double > &extremity1) const
int flipEdges() const
int snapVertexBarycentrics() const
bool isEdgeAngleCollapsible(const SimplexId &source, const SimplexId &destination, const SimplexId &pivotVertexId, const std::vector< std::pair< SimplexId, SimplexId > > &starNeighbors) const
const float * pointSet_
Definition: FiberSurface.h:484
bool isIntersectionTriangleColinear(const SimplexId &tetId, const SimplexId &triangleId, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, const std::vector< std::vector< Vertex > > &tetNewVertices, const SimplexId &vertexId0, const SimplexId &vertexId1, const SimplexId &vertexId2) const
Definition: FiberSurface.h:430
bool hasDuplicatedVertices(const double *p0, const double *p1, const double *p2) const
int mergeVertices(const double &distanceThreshold) const
std::vector< Vertex > * globalVertexList_
Definition: FiberSurface.h:496
int computeTriangleFiber(const SimplexId &tetId, const SimplexId &triangleId, const std::pair< double, double > &intersection, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::array< double, 3 > &pA, std::array< double, 3 > &pB, SimplexId &pivotVertexId, bool &edgeFiber) const
const SimplexId * tetList_
Definition: FiberSurface.h:485
int snapToBasePoint(const std::vector< std::vector< double > > &basePoints, const std::vector< std::pair< double, double > > &uv, const std::vector< double > &t, Vertex &v) const
int computeTriangleIntersection(const SimplexId &tetId, const SimplexId &triangleId0, const SimplexId &triangleId1, const SimplexId &polygonEdgeId0, const SimplexId &polygonEdgeId1, const std::pair< double, double > &intersection, SimplexId &newVertexNumber, SimplexId &newTriangleNumber, std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::vector< std::vector< Vertex > > &tetNewVertices) const
bool isEdgeFlippable(const SimplexId &edgeVertexId0, const SimplexId &edgeVertexId1, const SimplexId &otherVertexId0, const SimplexId &otherVertexId1) const
SimplexId tetNumber_
Definition: FiberSurface.h:482
int createNewIntersectionTriangle(const SimplexId &tetId, const SimplexId &triangleId, const SimplexId &vertexId0, const SimplexId &vertexId1, const SimplexId &vertexId2, const std::vector< std::vector< Vertex > > &tetNewVertices, SimplexId &newTriangleNumber, std::vector< std::vector< IntersectionTriangle > > &tetIntersections, const std::pair< double, double > *intersection=nullptr) const
double getElapsedTime()
Definition: Timer.h:15
bool isPointOnSegment(const T &x, const T &y, const T &xA, const T &yA, const T &xB, const T &yB)
Definition: Geometry.cpp:425
int computeBarycentricCoordinates(const T *p0, const T *p1, const T *p, std::array< T, 2 > &baryCentrics, const int &dimension=3)
Definition: Geometry.cpp:124
T1 pow(const T1 val, const T2 n)
Definition: Geometry.h:411
T angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1)
Definition: Geometry.cpp:13
int computeTriangleAngles(const T *p0, const T *p1, const T *p2, std::array< T, 3 > &angles)
Definition: Geometry.cpp:290
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition: Geometry.cpp:344
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
std::array< SimplexId, 3 > vertexIds_
Definition: FiberSurface.h:64
std::pair< double, double > uv_
Definition: FiberSurface.h:60
std::array< double, 3 > p_
Definition: FiberSurface.h:58