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
11FiberSurface::FiberSurface() {
12 this->setDebugMsgPrefix("FiberSurface");
13}
14
15int FiberSurface::getNumberOfCommonVertices(
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
54int FiberSurface::computeTriangleFiber(
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
165 Geometry::computeBarycentricCoordinates(
166 p0.data(), p1.data(), p.data(), baryCentrics0, 2);
167 Geometry::computeBarycentricCoordinates(
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
189int FiberSurface::computeTriangleIntersection(
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 const 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) {
316 computeTriangleIntersection(
317 tetId, triangleId0, polygonEdgeId0, intersection, pA, pB, pivotVertexIda,
318 newVertexNumber, newTriangleNumber, tetIntersections, tetNewVertices);
319 }
320
321 if(!edgeFiber1) {
322 computeTriangleIntersection(
323 tetId, triangleId1, polygonEdgeId1, intersection, pA, pB, pivotVertexIdb,
324 newVertexNumber, newTriangleNumber, tetIntersections, tetNewVertices);
325 }
326
327 return 0;
328}
329
330int FiberSurface::computeTriangleIntersection(
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{};
374 Geometry::computeBarycentricCoordinates(
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
379 Geometry::computeBarycentricCoordinates(
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 const 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 const 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
467 SimplexId const ret = createNewIntersectionTriangle(
468 tetId, triangleId, -vertexIdA, -vertexIdB, (pivotVertexId + 2) % 3,
469 tetNewVertices, newTriangleNumber, tetIntersections, &intersection);
470 if(ret == -1) {
471 // flip the diagonal
472 createNewIntersectionTriangle(
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
528int FiberSurface::createNewIntersectionTriangle(
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
600int FiberSurface::flipEdges() const {
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 const tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
613
614 if(!inQueue[tetId]) {
615 inQueue[tetId] = true;
616 tetList.push_back(tetId);
617 }
618
619 tetTriangles[tetId].emplace_back(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 const 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
642int FiberSurface::flipEdges(
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{};
676 Geometry::computeTriangleAngles(
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{};
717 Geometry::computeTriangleAngles(
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 const 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
790 Geometry::computeTriangleAngles(
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
802 Geometry::computeTriangleAngles(
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
874int FiberSurface::getTriangleRangeExtremities(
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
920 Geometry::computeBarycentricCoordinates(
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
945bool FiberSurface::hasDuplicatedVertices(const double *p0,
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
984bool FiberSurface::isEdgeAngleCollapsible(
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 const 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
1030bool FiberSurface::isEdgeFlippable(const SimplexId &edgeVertexId0,
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 const 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 const vertexId
1094 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
1095
1096 vertexTriangleNeighbors[vertexId].resize(
1097 vertexTriangleNeighbors[vertexId].size() + 1);
1098 vertexTriangleNeighbors[vertexId].back().first = -1;
1099 vertexTriangleNeighbors[vertexId].back().second = -1;
1100 vertexNeighborsTets[vertexId].push_back(
1101 (*polygonEdgeTriangleLists_[i])[j].tetId_);
1102
1103 for(int l = 0; l < 3; l++) {
1104 if(l != k) {
1105 SimplexId const otherVertexId
1106 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[l];
1107
1108 if(vertexTriangleNeighbors[vertexId].back().first == -1) {
1109 vertexTriangleNeighbors[vertexId].back().first = otherVertexId;
1110 } else {
1111 vertexTriangleNeighbors[vertexId].back().second = otherVertexId;
1112 }
1113
1114 bool isIn = false;
1115 for(SimplexId m = 0;
1116 m < (SimplexId)vertexNeighbors[vertexId].size(); m++) {
1117 if(vertexNeighbors[vertexId][m] == otherVertexId) {
1118 isIn = true;
1119 break;
1120 }
1121 }
1122 if(!isIn) {
1123 vertexNeighbors[vertexId].push_back(otherVertexId);
1124 }
1125 }
1126 }
1127 }
1128 }
1129 }
1130
1131 bool hasMerged = false;
1132
1133#ifdef TTK_ENABLE_OPENMP
1134#pragma omp parallel for num_threads(threadNumber_)
1135#endif
1136 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1137
1138 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1139 j++) {
1140
1141 SimplexId minimizer = -1;
1142 double minDistance = -1;
1143
1144 // find the smallest edge on this triangle
1145 for(int k = 0; k < 3; k++) {
1146
1147 SimplexId const vertexId0
1148 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
1149 SimplexId const vertexId1
1150 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[(k + 1) % 3];
1151
1152 bool areAlreadySnapped = true;
1153 for(int l = 0; l < 3; l++) {
1154 if((*globalVertexList_)[vertexId0].p_[l]
1155 != (*globalVertexList_)[vertexId1].p_[l]) {
1156
1157 areAlreadySnapped = false;
1158 break;
1159 }
1160 }
1161
1162 if(((*globalVertexList_)[vertexId0].isBasePoint_)
1163 && ((*globalVertexList_)[vertexId1].isBasePoint_)) {
1164 areAlreadySnapped = true;
1165 }
1166
1167 if(!areAlreadySnapped) {
1168 double const distance
1169 = Geometry::distance((*globalVertexList_)[vertexId0].p_.data(),
1170 (*globalVertexList_)[vertexId1].p_.data());
1171
1172 if((minDistance == -1) || (distance < minDistance)) {
1173 minDistance = distance;
1174 minimizer = k;
1175 }
1176 }
1177 }
1178
1179 if((minDistance != -1) && (minDistance < distanceThreshold)) {
1180 SimplexId const vertexId0
1181 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[minimizer];
1182 SimplexId const vertexId1 = (*polygonEdgeTriangleLists_[i])[j]
1183 .vertexIds_[(minimizer + 1) % 3];
1184
1185 // find the number of common neighbors
1186 vector<SimplexId> commonNeighbors;
1187 for(SimplexId k = 0; k < (SimplexId)vertexNeighbors[vertexId0].size();
1188 k++) {
1189 for(SimplexId l = 0;
1190 l < (SimplexId)vertexNeighbors[vertexId1].size(); l++) {
1191 if(vertexNeighbors[vertexId0][k]
1192 == vertexNeighbors[vertexId1][l]) {
1193 commonNeighbors.push_back(vertexNeighbors[vertexId0][k]);
1194 }
1195 }
1196 }
1197
1198 if(commonNeighbors.size() == 2) {
1199
1200 // we may create extra non-manifold triangles otherwise
1201 // plus we don't collapse edges that are strictly inside a tet
1202 // (only collsapse those on the boundary between two tets)
1203
1204 SimplexId source = vertexId0, destination = vertexId1;
1205
1206 if((*globalVertexList_)[destination].isBasePoint_) {
1207 source = vertexId1;
1208 destination = vertexId0;
1209 }
1210
1211 if(!(*globalVertexList_)[destination].isBasePoint_) {
1212
1213 // now check the angles
1214 bool isCollapsible = true;
1215 for(SimplexId k = 0; k < (SimplexId)commonNeighbors.size(); k++) {
1216 if(!isEdgeAngleCollapsible(
1217 source, destination, commonNeighbors[k],
1218 vertexTriangleNeighbors[commonNeighbors[k]])) {
1219 isCollapsible = false;
1220 break;
1221 }
1222 }
1223
1224 // different tets?
1225 SimplexId tetId0 = -1, tetId1 = -1;
1226 for(SimplexId k = 0;
1227 k < (SimplexId)vertexNeighborsTets[source].size(); k++) {
1228 if((vertexTriangleNeighbors[source][k].first == destination)
1229 || (vertexTriangleNeighbors[source][k].second
1230 == destination)) {
1231 if(tetId0 == -1) {
1232 tetId0 = vertexNeighborsTets[source][k];
1233 } else {
1234 tetId1 = vertexNeighborsTets[source][k];
1235 }
1236 }
1237 }
1238
1239 if(tetId0 == tetId1) {
1240 isCollapsible = false;
1241 }
1242
1243#ifdef TTK_ENABLE_OPENMP
1244#pragma omp critical
1245#endif
1246 if(isCollapsible) {
1247 for(int k = 0; k < 3; k++) {
1248 (*globalVertexList_)[destination].p_[k]
1249 = (*globalVertexList_)[source].p_[k];
1250 }
1251
1252 hasMerged = true;
1253 }
1254 }
1255 }
1256 }
1257 }
1258 }
1259
1260 if(!hasMerged)
1261 break;
1262
1263 // now update the vertices
1264 mergeVertices(0);
1265 }
1266
1267 this->printMsg(
1268 "Performed edge collapses", 1.0, t.getElapsedTime(), this->threadNumber_);
1269 this->printMsg(std::vector<std::vector<std::string>>{
1270 {"#Vertices removed",
1271 std::to_string(initVertexNumber - (*globalVertexList_).size())}});
1272
1273 return 0;
1274}
1275
1276int FiberSurface::mergeVertices(const double &distanceThreshold) const {
1277
1278 Timer t;
1279
1280 vector<Vertex> tmpList = (*globalVertexList_);
1281
1282 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1283 tmpList[i].localId_ = i;
1284 }
1285
1286 const auto FiberSurfaceVertexComparisonX
1287 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1288 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1289 // let's consider x coordinates are equal
1290 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1291 // let's consider y coordinates are equal
1292 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1293 // let's consider z coordinates are equal
1294 // NOTE: the local Id should be sufficient
1295 return v0.globalId_ < v1.globalId_;
1296 } else
1297 return v0.p_[2] < v1.p_[2];
1298 } else
1299 return v0.p_[1] < v1.p_[1];
1300 } else {
1301 return v0.p_[0] < v1.p_[0];
1302 }
1303 };
1304
1305 const auto FiberSurfaceVertexComparisonY
1306 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1307 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1308 // let's consider y coordinates are equal
1309 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1310 // let's consider z coordinates are equal
1311 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1312 // let's consider x coordinates are equal
1313 // NOTE: the local Id should be sufficient
1314 return v0.globalId_ < v1.globalId_;
1315 } else
1316 return v0.p_[0] < v1.p_[0];
1317 } else
1318 return v0.p_[2] < v1.p_[2];
1319 } else {
1320 return v0.p_[1] < v1.p_[1];
1321 }
1322 };
1323
1324 const auto FiberSurfaceVertexComparisonZ
1325 = [](const FiberSurface::Vertex &v0, const FiberSurface::Vertex &v1) {
1326 if(fabs(v0.p_[2] - v1.p_[2]) < PREC_DBL) {
1327 // let's consider z coordinates are equal
1328 if(fabs(v0.p_[0] - v1.p_[0]) < PREC_DBL) {
1329 // let's consider x coordinates are equal
1330 if(fabs(v0.p_[1] - v1.p_[1]) < PREC_DBL) {
1331 // let's consider y coordinates are equal
1332 // NOTE: the local Id should be sufficient
1333 return v0.globalId_ < v1.globalId_;
1334 } else
1335 return v0.p_[1] < v1.p_[1];
1336 } else
1337 return v0.p_[0] < v1.p_[0];
1338 } else {
1339 return v0.p_[2] < v1.p_[2];
1340 }
1341 };
1342
1343 // now do a parallel sort
1344 SimplexId uniqueVertexNumber = 0;
1345 for(int k = 0; k < 3; k++) {
1346
1347 switch(k) {
1348 case 0:
1349 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonX);
1350 break;
1351
1352 case 1:
1353 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonY);
1354 break;
1355
1356 case 2:
1357 sort(tmpList.begin(), tmpList.end(), FiberSurfaceVertexComparisonZ);
1358 break;
1359 }
1360
1361 // NOTE:
1362 // ethaneDiol dataset
1363 // 0.019224 in parallel (2 cores HT)
1364 // 0.041353 in sequential (speedup: x2.15, perfect (HT))
1365
1366 // now merge the thing
1367 // 1. Identity duplicates
1368 // NOTE: order is important, so no parallelism
1369 // NOTE: points are represented with single precision (floats) so use that
1370 // as a limit for distances
1371 uniqueVertexNumber = 0;
1372 double distance = 0;
1373 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1374
1375 bool canMerge = false;
1376
1377 if(i) {
1378 distance
1379 = Geometry::distance(tmpList[i].p_.data(), tmpList[i - 1].p_.data());
1380
1381 if(distance <= distanceThreshold) {
1382
1383 // one of the two is -1 (interior vertex)
1384 // if((tmpList[i - 1].meshEdge_.first == -1)
1385 // ||(tmpList[i].meshEdge_.first == -1)){
1386 // canMerge = true;
1387 // }
1388 //
1389 // // one vertex in common
1390 // if((tmpList[i].meshEdge_.first ==
1391 // tmpList[i - 1].meshEdge_.first)
1392 // ||
1393 // (tmpList[i].meshEdge_.first ==
1394 // tmpList[i - 1].meshEdge_.second)){
1395 // canMerge = true;
1396 // }
1397 //
1398 // // one vertex in common
1399 // if((tmpList[i].meshEdge_.second ==
1400 // tmpList[i - 1].meshEdge_.first)
1401 // ||
1402 // (tmpList[i].meshEdge_.second ==
1403 // tmpList[i - 1].meshEdge_.second)){
1404 // canMerge = true;
1405 // }
1406
1407 // NOTE: still some bugs here in terms of manifoldness.
1408
1409 if((tmpList[i].meshEdge_.first == tmpList[i - 1].meshEdge_.first)
1410 && (tmpList[i].meshEdge_.second
1411 == tmpList[i - 1].meshEdge_.second)) {
1412 canMerge = true;
1413 }
1414 if((tmpList[i].meshEdge_.first == tmpList[i - 1].meshEdge_.second)
1415 && (tmpList[i].meshEdge_.second
1416 == tmpList[i - 1].meshEdge_.first)) {
1417 canMerge = true;
1418 }
1419
1420 // canMerge = true;
1421 }
1422
1423 if(canMerge) {
1424 tmpList[i].globalId_ = tmpList[i - 1].globalId_;
1425 if((tmpList[i].isBasePoint_) || (tmpList[i - 1].isBasePoint_)) {
1426 tmpList[i].isBasePoint_ = true;
1427 tmpList[i - 1].isBasePoint_ = true;
1428 }
1429 if((tmpList[i].isIntersectionPoint_)
1430 || (tmpList[i - 1].isIntersectionPoint_)) {
1431 tmpList[i].isIntersectionPoint_ = true;
1432 tmpList[i - 1].isIntersectionPoint_ = true;
1433 }
1434 for(int j = 0; j < 3; j++) {
1435 tmpList[i].p_[j] = tmpList[i - 1].p_[j];
1436 }
1437
1438 // update meshEdge...
1439 if((tmpList[i - 1].meshEdge_.first == -1)
1440 && (tmpList[i].meshEdge_.first != -1)) {
1441 tmpList[i - 1].meshEdge_ = tmpList[i].meshEdge_;
1442 }
1443 if((tmpList[i].meshEdge_.first == -1)
1444 && (tmpList[i - 1].meshEdge_.first != -1)) {
1445 tmpList[i].meshEdge_ = tmpList[i - 1].meshEdge_;
1446 }
1447 }
1448 }
1449
1450 if((!i) || (!canMerge)) {
1451 tmpList[i].globalId_ = uniqueVertexNumber;
1452 uniqueVertexNumber++;
1453 }
1454
1455 (*globalVertexList_)[tmpList[i].localId_].globalId_
1456 = tmpList[i].globalId_;
1457 }
1458 }
1459
1460 // update the triangles
1461 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1462 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1463 j++) {
1464 for(int k = 0; k < 3; k++) {
1465 (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k]
1466 = (*globalVertexList_)[(*polygonEdgeTriangleLists_[i])[j]
1467 .vertexIds_[k]]
1468 .globalId_;
1469 }
1470 }
1471 }
1472
1473 // 2. create the actual global list
1474 // NOTE: order is no longer important, we can do this in parallel.
1475 (*globalVertexList_).resize(uniqueVertexNumber);
1476
1477 // duplicate vertices but only a subset have the right edge information
1478 // (because they have actually been computed on the edge)
1479 SimplexId lastId = -1;
1480 for(SimplexId i = 0; i < (SimplexId)tmpList.size(); i++) {
1481
1482 if(lastId != -1) {
1483 if((tmpList[lastId].globalId_ == tmpList[i].globalId_)
1484 && (tmpList[i].meshEdge_.first != -1)) {
1485
1486 (*globalVertexList_)[tmpList[lastId].globalId_].meshEdge_
1487 = tmpList[i].meshEdge_;
1488 (*globalVertexList_)[tmpList[lastId].globalId_].isBasePoint_ = true;
1489 lastId = -1;
1490 }
1491 }
1492
1493 if((!i) || (tmpList[i].globalId_ != tmpList[i - 1].globalId_)) {
1494 (*globalVertexList_)[tmpList[i].globalId_] = tmpList[i];
1495 if((*globalVertexList_)[tmpList[i].globalId_].meshEdge_.first == -1) {
1496 lastId = i;
1497 } else {
1498 lastId = -1;
1499 }
1500 }
1501 }
1502
1503 // 3. update the 2-sheets, ignore zero-area triangles
1504 // and free their vertex memory
1505 // NOTE: order is no longer important, parallel
1506 vector<vector<bool>> keepTriangle(polygonEdgeTriangleLists_.size());
1507 vector<vector<FiberSurface::Triangle>> tmpTriangleLists(
1508 polygonEdgeTriangleLists_.size());
1509 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1510 tmpTriangleLists[i] = (*polygonEdgeTriangleLists_[i]);
1511 keepTriangle[i].resize((*polygonEdgeTriangleLists_[i]).size(), true);
1512 }
1513
1514#ifdef TTK_ENABLE_OPENMP
1515#pragma omp parallel for num_threads(threadNumber_)
1516#endif
1517 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1518 for(SimplexId j = 0; j < (SimplexId)tmpTriangleLists[i].size(); j++) {
1519 for(int k = 0; k < 3; k++) {
1520 if((k)
1521 && (tmpTriangleLists[i][j].vertexIds_[k]
1522 == tmpTriangleLists[i][j].vertexIds_[(k - 1)])) {
1523 keepTriangle[i][j] = false;
1524 }
1525 }
1526 if(tmpTriangleLists[i][j].vertexIds_[0]
1527 == tmpTriangleLists[i][j].vertexIds_[2]) {
1528 keepTriangle[i][j] = false;
1529 }
1530 }
1531
1532 // now copy triangles over with non zero-area triangles
1533 // NOTE: no need to re-allocate the memory, we know we are not going to use
1534 // more.
1535 (*polygonEdgeTriangleLists_[i]).clear();
1536 for(SimplexId j = 0; j < (SimplexId)tmpTriangleLists[i].size(); j++) {
1537
1538 if(keepTriangle[i][j]) {
1539 (*polygonEdgeTriangleLists_[i]).push_back(tmpTriangleLists[i][j]);
1540 }
1541 }
1542 }
1543
1544 this->printMsg(
1545 "Output made manifold", 1.0, t.getElapsedTime(), this->threadNumber_);
1546
1547 return 0;
1548}
1549
1550int FiberSurface::snapToBasePoint(const vector<vector<double>> &basePoints,
1551 const vector<pair<double, double>> &uv,
1552 const vector<double> &t,
1553 Vertex &v) const {
1554
1555 if(!pointSnappingThreshold_)
1556 return -1;
1557
1558 SimplexId minimizer = 0;
1559 double minDistance = -1;
1560
1561 for(SimplexId i = 0; i < (SimplexId)basePoints.size(); i++) {
1562 double const distance
1563 = Geometry::distance(basePoints[i].data(), v.p_.data());
1564 if((minDistance < 0) || (distance < minDistance)) {
1565 minDistance = distance;
1566 minimizer = i;
1567 }
1568 }
1569
1570 if(minDistance < pointSnappingThreshold_) {
1571 for(int i = 0; i < 3; i++) {
1572 v.p_[i] = basePoints[minimizer][i];
1573 }
1574 v.uv_ = uv[minimizer];
1575 v.t_ = t[minimizer];
1576 v.isBasePoint_ = true;
1577 }
1578
1579 return 0;
1580}
1581
1582int FiberSurface::snapVertexBarycentrics() const {
1583
1584 vector<bool> inQueue(tetNumber_, false);
1585 vector<vector<pair<SimplexId, SimplexId>>> tetTriangles(tetNumber_);
1586 vector<SimplexId> tetList;
1587
1588 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1589
1590 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1591 j++) {
1592
1593 SimplexId const tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
1594
1595 if(!inQueue[tetId]) {
1596 inQueue[tetId] = true;
1597 tetList.push_back(tetId);
1598 }
1599
1600 tetTriangles[tetId].emplace_back(i, j);
1601 }
1602 }
1603
1604#ifdef TTK_ENABLE_OPENMP
1605#pragma omp parallel for num_threads(threadNumber_)
1606#endif
1607 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
1608 snapVertexBarycentrics(tetList[i], tetTriangles[tetList[i]]);
1609 }
1610
1611 mergeVertices(0);
1612
1613 return 0;
1614}
1615
1616int FiberSurface::snapVertexBarycentrics(
1617 const SimplexId &tetId,
1618 const std::vector<std::pair<SimplexId, SimplexId>> &triangles) const {
1619
1620 for(SimplexId i = 0; i < (SimplexId)triangles.size(); i++) {
1621
1622 Triangle *t = &(
1623 (*polygonEdgeTriangleLists_[triangles[i].first])[triangles[i].second]);
1624
1625 for(int j = 0; j < 3; j++) {
1626 SimplexId const vertexId = t->vertexIds_[j];
1627
1628 // check for each triangle of the tet
1629 double minimum = -DBL_MAX;
1630 std::array<double, 3> minBarycentrics{};
1631 std::array<SimplexId, 3> minimizer{};
1632
1633 for(int k = 0; k < 2; k++) {
1634 for(int l = k + 1; l < 3; l++) {
1635 for(int m = l + 1; m < 4; m++) {
1636
1637 SimplexId const vertexId0 = tetList_[5 * tetId + 1 + k];
1638 SimplexId const vertexId1 = tetList_[5 * tetId + 1 + l];
1639 SimplexId const vertexId2 = tetList_[5 * tetId + 1 + m];
1640
1641 std::array<double, 3> p0{}, p1{}, p2{};
1642
1643 for(int n = 0; n < 3; n++) {
1644 p0[n] = pointSet_[3 * vertexId0 + n];
1645 p1[n] = pointSet_[3 * vertexId1 + n];
1646 p2[n] = pointSet_[3 * vertexId2 + n];
1647 }
1648
1649 std::array<double, 3> barycentrics{};
1650 Geometry::computeBarycentricCoordinates(
1651 p0.data(), p1.data(), p2.data(),
1652 (*globalVertexList_)[vertexId].p_.data(), barycentrics);
1653
1654 if((barycentrics[0] != -1.0) && (barycentrics[1] != -1.0)
1655 && (barycentrics[2] != -1.0)) {
1656 // vertexId lies in the current triangle
1657
1658 double localMin = -1.0;
1659 for(int n = 0; n < 3; n++) {
1660 if((localMin == -1.0) || (fabs(barycentrics[n]) < localMin)) {
1661 localMin = fabs(barycentrics[n]);
1662 }
1663 }
1664
1665 if((minimum == -DBL_MAX) || (localMin < minimum)) {
1666 minimum = localMin;
1667 minBarycentrics = barycentrics;
1668 minimizer[0] = vertexId0;
1669 minimizer[1] = vertexId1;
1670 minimizer[2] = vertexId2;
1671 }
1672 }
1673 }
1674 }
1675 }
1676
1677 if((minimum != -DBL_MAX) && (minimum < PREC_FLT_2)) {
1678 double sum = 0;
1679 int numberOfZeros = 0;
1680 for(int k = 0; k < 3; k++) {
1681 if(minBarycentrics[k] < PREC_FLT_2) {
1682 minBarycentrics[k] = 0;
1683 numberOfZeros++;
1684 }
1685 sum += minBarycentrics[k];
1686 }
1687
1688 sum = (1 - sum) / numberOfZeros;
1689
1690 for(int k = 0; k < 3; k++) {
1691 if(minBarycentrics[k] >= PREC_FLT_2) {
1692 minBarycentrics[k] += sum;
1693 }
1694 }
1695
1696 // now do the re-location
1697#ifdef TTK_ENABLE_OPENMP
1698#pragma omp critical
1699 for(int k = 0; k < 3; k++) {
1700 (*globalVertexList_)[vertexId].p_[k]
1701 = minBarycentrics[0] * ((double)pointSet_[3 * minimizer[0] + k])
1702 + minBarycentrics[1] * ((double)pointSet_[3 * minimizer[1] + k])
1703 + minBarycentrics[2] * ((double)pointSet_[3 * minimizer[2] + k]);
1704 }
1705#endif
1706 }
1707 }
1708 }
1709
1710 return 0;
1711}
M_PI
#define M_PI
Definition Os.h:50
ttk::Debug::setDebugMsgPrefix
void setDebugMsgPrefix(const std::string &prefix)
Definition Debug.h:364
ttk::FiberSurface::interpolateBasePoints
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
Definition FiberSurface.cpp:961
ttk::FiberSurface::pointSnappingThreshold_
double pointSnappingThreshold_
Definition FiberSurface.h:491
ttk::FiberSurface::polygonEdgeTriangleLists_
std::vector< std::vector< Triangle > * > polygonEdgeTriangleLists_
Definition FiberSurface.h:498
ttk::FiberSurface::getNumberOfCommonVertices
int getNumberOfCommonVertices(const SimplexId &tetId, const SimplexId &triangleId0, const SimplexId &triangleId1, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections) const
Definition FiberSurface.cpp:15
ttk::FiberSurface::mergeEdges
int mergeEdges(const double &distanceThreshold) const
Definition FiberSurface.cpp:1067
ttk::FiberSurface::getTriangleRangeExtremities
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
Definition FiberSurface.cpp:874
ttk::FiberSurface::flipEdges
int flipEdges() const
Definition FiberSurface.cpp:600
ttk::FiberSurface::snapVertexBarycentrics
int snapVertexBarycentrics() const
Definition FiberSurface.cpp:1582
ttk::FiberSurface::isEdgeAngleCollapsible
bool isEdgeAngleCollapsible(const SimplexId &source, const SimplexId &destination, const SimplexId &pivotVertexId, const std::vector< std::pair< SimplexId, SimplexId > > &starNeighbors) const
Definition FiberSurface.cpp:984
ttk::FiberSurface::pointSet_
const float * pointSet_
Definition FiberSurface.h:484
ttk::FiberSurface::isIntersectionTriangleColinear
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
ttk::FiberSurface::hasDuplicatedVertices
bool hasDuplicatedVertices(const double *p0, const double *p1, const double *p2) const
Definition FiberSurface.cpp:945
ttk::FiberSurface::mergeVertices
int mergeVertices(const double &distanceThreshold) const
Definition FiberSurface.cpp:1276
ttk::FiberSurface::globalVertexList_
std::vector< Vertex > * globalVertexList_
Definition FiberSurface.h:496
ttk::FiberSurface::computeTriangleFiber
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
Definition FiberSurface.cpp:54
ttk::FiberSurface::tetList_
const SimplexId * tetList_
Definition FiberSurface.h:485
ttk::FiberSurface::snapToBasePoint
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
Definition FiberSurface.cpp:1550
ttk::FiberSurface::computeTriangleIntersection
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
Definition FiberSurface.cpp:189
ttk::FiberSurface::isEdgeFlippable
bool isEdgeFlippable(const SimplexId &edgeVertexId0, const SimplexId &edgeVertexId1, const SimplexId &otherVertexId0, const SimplexId &otherVertexId1) const
Definition FiberSurface.cpp:1030
ttk::FiberSurface::createNewIntersectionTriangle
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
Definition FiberSurface.cpp:528
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
ttk::Geometry::isPointOnSegment
bool isPointOnSegment(const T &x, const T &y, const T &xA, const T &yA, const T &xB, const T &yB)
Definition Geometry.cpp:443
ttk::Geometry::computeBarycentricCoordinates
int computeBarycentricCoordinates(const T *p0, const T *p1, const T *p, std::array< T, 2 > &baryCentrics, const int &dimension=3)
Definition Geometry.cpp:142
ttk::Geometry::pow
T1 pow(const T1 val, const T2 n)
Definition Geometry.h:456
ttk::Geometry::angle
T angle(const T *vA0, const T *vA1, const T *vB0, const T *vB1)
Definition Geometry.cpp:13
ttk::Geometry::computeTriangleAngles
int computeTriangleAngles(const T *p0, const T *p1, const T *p2, std::array< T, 3 > &angles)
Definition Geometry.cpp:308
ttk::Geometry::distance
T distance(const T *p0, const T *p1, const int &dimension=3)
Definition Geometry.cpp:362
ttk
The Topology ToolKit.
Definition AbstractTriangulation.h:51
ttk::SimplexId
int SimplexId
Identifier type for simplices of any dimension.
Definition DataTypes.h:22
ttk::FiberSurface::Triangle::vertexIds_
std::array< SimplexId, 3 > vertexIds_
Definition FiberSurface.h:64
ttk::FiberSurface::Vertex::uv_
std::pair< double, double > uv_
Definition FiberSurface.h:60
ttk::FiberSurface::Vertex::p_
std::array< double, 3 > p_
Definition FiberSurface.h:58
printMsg
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)
Generated on Thu Sep 26 2024 06:50:18 for TTK by doxygen 1.9.8