TTK
Loading...
Searching...
No Matches
FiberSurface.h
Go to the documentation of this file.
1
31
32#pragma once
33
34// standard includes
35#include <array>
36#include <queue>
37
38// base code includes
39#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
40#include <RangeDrivenOctree.h>
41#endif
42
43#include <Debug.h>
44#include <Geometry.h>
45#include <Triangulation.h>
46
47namespace ttk {
48
49 class FiberSurface : virtual public Debug {
50
51 public:
52 struct Vertex {
55 // TODO also encode the vertex ids of the triangle of the input mesh
56 // where this point has been computed (for constrained triangulation)
57 std::pair<SimplexId, SimplexId> meshEdge_{};
58 std::array<double, 3> p_{};
59 double t_{};
60 std::pair<double, double> uv_{};
61 };
62
63 struct Triangle {
64 std::array<SimplexId, 3> vertexIds_{};
66 };
67
69
70#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
71 template <class dataTypeU, class dataTypeV, typename triangulationType>
72 inline int buildOctree(const triangulationType *const triangulation);
73#endif
74
75 template <class dataTypeU, class dataTypeV, typename triangulationType>
76 inline int computeContour(const std::pair<double, double> &rangePoint0,
77 const std::pair<double, double> &rangePoint1,
78 const std::vector<SimplexId> &seedTetList,
79 const triangulationType *const triangulation,
80 const SimplexId &polygonEdgeId = 0) const;
81
82 template <class dataTypeU, class dataTypeV>
83 inline int computeContour(
84 const std::vector<std::pair<std::pair<double, double>,
85 std::pair<double, double>>> &edgeList,
86 const std::vector<SimplexId> &seedTetList,
87 const std::vector<SimplexId> *edgeIdList = nullptr) const;
88
89 template <class dataTypeU, class dataTypeV, typename triangulationType>
90 inline int computeSurface(const std::pair<double, double> &rangePoint0,
91 const std::pair<double, double> &rangePoint1,
92 const triangulationType *const triangulation,
93 const SimplexId &polygonEdgeId = 0) const;
94
95 template <class dataTypeU, class dataTypeV, typename triangulationType>
96 inline int computeSurface(const triangulationType *const triangulation);
97
98#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
99 template <class dataTypeU, class dataTypeV, typename triangulationType>
100 inline int
101 computeSurfaceWithOctree(const std::pair<double, double> &rangePoint0,
102 const std::pair<double, double> &rangePoint1,
103 const triangulationType *const triangulation,
104 const SimplexId &polygonEdgeId) const;
105#endif
106
107 template <class dataTypeU, class dataTypeV>
108 inline int finalize(const bool &mergeDuplicatedVertices = false,
109 const bool &removeSmallEdges = false,
110 const bool &edgeFlips = false,
111 const bool &intersectionRemesh = false);
112
113#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
114 inline void flushOctree() {
115 octree_.flush();
116 }
117#endif
118
119 template <class dataTypeU, class dataTypeV, typename triangulationType>
120 inline int processTetrahedron(const SimplexId &tetId,
121 const std::pair<double, double> &rangePoint0,
122 const std::pair<double, double> &rangePoint1,
123 const triangulationType *const triangulation,
124 const SimplexId &polygonEdgeId = 0) const;
125
126 inline int setGlobalVertexList(std::vector<Vertex> *globalList) {
127 globalVertexList_ = globalList;
128 return 0;
129 }
130
131 inline int setInputField(const void *uField, const void *vField) {
132
133 uField_ = uField;
134 vField_ = vField;
135
136 return 0;
137 }
138
139 inline int setPointMerging(const bool &onOff) {
140 pointSnapping_ = onOff;
141 return 0;
142 }
143
144 inline int setPointMergingThreshold(const double &threshold) {
145 pointSnappingThreshold_ = threshold;
146 return 0;
147 }
148
149 inline int setPointNumber(const SimplexId &number) {
150 pointNumber_ = number;
151 return 0;
152 }
153
154 inline int setPointSet(const float *pointSet) {
155 pointSet_ = pointSet;
156 return 0;
157 }
158
159 inline int setPolygon(
160 const std::vector<std::pair<std::pair<double, double>,
161 std::pair<double, double>>> *polygon) {
162 polygon_ = polygon;
163 return 0;
164 }
165
166 inline int setPolygonEdgeNumber(const SimplexId &polygonEdgeNumber) {
167 polygonEdgeNumber_ = polygonEdgeNumber;
168 polygonEdgeVertexLists_.resize(polygonEdgeNumber, nullptr);
169 polygonEdgeTriangleLists_.resize(polygonEdgeNumber, nullptr);
170 return 0;
171 }
172
173 inline int setTetList(const SimplexId *tetList) {
174 tetList_ = tetList;
175 return 0;
176 }
177
178 inline int
179 setTetNeighbors(const std::vector<std::vector<SimplexId>> *tetNeighbors) {
180 tetNeighbors_ = tetNeighbors;
181 return 0;
182 }
183
184 inline int setTetNumber(const SimplexId &tetNumber) {
185 tetNumber_ = tetNumber;
186 return 0;
187 }
188
189 inline int setTriangleList(const SimplexId &polygonEdgeId,
190 std::vector<Triangle> *triangleList) {
191
192#ifndef TTK_ENABLE_KAMIKAZE
193 if((polygonEdgeId >= 0)
194 && (polygonEdgeId < (SimplexId)polygonEdgeTriangleLists_.size()))
195#endif
196 polygonEdgeTriangleLists_[polygonEdgeId] = triangleList;
197
198 return 0;
199 }
200
201 inline void
203 // for breadth-first search traversals
204 if(triangulation != nullptr) {
205 triangulation->preconditionCellNeighbors();
206 }
207 }
208
209 inline int setVertexList(const SimplexId &polygonEdgeId,
210 std::vector<Vertex> *vertexList) {
211
212#ifndef TTK_ENABLE_KAMIKAZE
213 if((polygonEdgeId >= 0)
214 && (polygonEdgeId < (SimplexId)polygonEdgeVertexLists_.size()))
215#endif
216 polygonEdgeVertexLists_[polygonEdgeId] = vertexList;
217
218 return 0;
219 }
220
221 protected:
224 // use negative values for new triangles
227 // use negative values for new vertices
228 std::array<SimplexId, 3> vertexIds_{};
229 std::array<std::pair<double, double>, 3> uv_{};
230 std::array<double, 3> t_{};
231 std::array<std::array<double, 3>, 3> p_{};
232 std::pair<double, double> intersection_{};
233 };
234
235 template <class dataTypeU, class dataTypeV, typename triangulationType>
236 inline int computeBaseTriangle(
237 const SimplexId &tetId,
238 const SimplexId &localEdgeId0,
239 const double &t0,
240 const double &u0,
241 const double &v0,
242 const SimplexId &localEdgeId1,
243 const double &t1,
244 const double &u1,
245 const double &v1,
246 const SimplexId &localEdgeId2,
247 const double &t2,
248 const double &u2,
249 const double &v2,
250 std::array<std::array<double, 3>, 3> &basePoints,
251 std::array<std::pair<double, double>, 3> &basePointProjections,
252 std::array<double, 3> &basePointParameterization,
253 std::array<std::pair<SimplexId, SimplexId>, 3> &baseEdges,
254 const triangulationType *const triangulation) const;
255
256 template <class dataTypeU, class dataTYpeV, typename triangulationType>
257 inline int computeCase0(const SimplexId &polygonEdgeId,
258 const SimplexId &tetId,
259 const SimplexId &localEdgeId0,
260 const double &t0,
261 const double &u0,
262 const double &v0,
263 const SimplexId &localEdgeId1,
264 const double &t1,
265 const double &u1,
266 const double &v1,
267 const SimplexId &localEdgeId2,
268 const double &t2,
269 const double &u2,
270 const double &v2,
271 const triangulationType *const triangulation) const;
272
273 template <class dataTypeU, class dataTYpeV, typename triangulationType>
274 inline int computeCase1(const SimplexId &polygonEdgeId,
275 const SimplexId &tetId,
276 const SimplexId &localEdgeId0,
277 const double &t0,
278 const double &u0,
279 const double &v0,
280 const SimplexId &localEdgeId1,
281 const double &t1,
282 const double &u1,
283 const double &v1,
284 const SimplexId &localEdgeId2,
285 const double &t2,
286 const double &u2,
287 const double &v2,
288 const triangulationType *const triangulation) const;
289
290 template <class dataTypeU, class dataTYpeV, typename triangulationType>
291 inline int computeCase2(const SimplexId &polygonEdgeId,
292 const SimplexId &tetId,
293 const SimplexId &localEdgeId0,
294 const double &t0,
295 const double &u0,
296 const double &v0,
297 const SimplexId &localEdgeId1,
298 const double &t1,
299 const double &u1,
300 const double &v1,
301 const SimplexId &localEdgeId2,
302 const double &t2,
303 const double &u2,
304 const double &v2,
305 const triangulationType *const triangulation) const;
306
307 template <class dataTypeU, class dataTYpeV, typename triangulationType>
308 inline int computeCase3(const SimplexId &polygonEdgeId,
309 const SimplexId &tetId,
310 const SimplexId &localEdgeId0,
311 const double &t0,
312 const double &u0,
313 const double &v0,
314 const SimplexId &localEdgeId1,
315 const double &t1,
316 const double &u1,
317 const double &v1,
318 const SimplexId &localEdgeId2,
319 const double &t2,
320 const double &u2,
321 const double &v2,
322 const triangulationType *const triangulation) const;
323
324 template <class dataTypeU, class dataTYpeV, typename triangulationType>
325 inline int computeCase4(const SimplexId &polygonEdgeId,
326 const SimplexId &tetId,
327 const SimplexId &localEdgeId0,
328 const double &t0,
329 const double &u0,
330 const double &v0,
331 const SimplexId &localEdgeId1,
332 const double &t1,
333 const double &u1,
334 const double &v1,
335 const SimplexId &localEdgeId2,
336 const double &t2,
337 const double &u2,
338 const double &v2,
339 const triangulationType *const triangulation) const;
340
342 const SimplexId &tetId,
343 const SimplexId &triangleId,
344 const std::pair<double, double> &intersection,
345 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
346 std::array<double, 3> &pA,
347 std::array<double, 3> &pB,
348 SimplexId &pivotVertexId,
349 bool &edgeFiber) const;
350
352 const SimplexId &tetId,
353 const SimplexId &triangleId0,
354 const SimplexId &triangleId1,
355 const SimplexId &polygonEdgeId0,
356 const SimplexId &polygonEdgeId1,
357 const std::pair<double, double> &intersection,
358 SimplexId &newVertexNumber,
359 SimplexId &newTriangleNumber,
360 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
361 std::vector<std::vector<Vertex>> &tetNewVertices) const;
362
364 const SimplexId &tetId,
365 const SimplexId &triangleId,
366 const SimplexId &polygonEdgeId,
367 const std::pair<double, double> &intersection,
368 const std::array<double, 3> &pA,
369 const std::array<double, 3> &pB,
370 const SimplexId &pivotVertexId,
371 SimplexId &newVertexNumber,
372 SimplexId &newTriangleNumber,
373 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
374 std::vector<std::vector<Vertex>> &tetNewVertices) const;
375
377 const SimplexId &tetId,
378 const SimplexId &triangleId,
379 const SimplexId &vertexId0,
380 const SimplexId &vertexId1,
381 const SimplexId &vertexId2,
382 const std::vector<std::vector<Vertex>> &tetNewVertices,
383 SimplexId &newTriangleNumber,
384 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
385 const std::pair<double, double> *intersection = nullptr) const;
386
387 int flipEdges() const;
388
389 int
390 flipEdges(std::vector<std::pair<SimplexId, SimplexId>> &triangles) const;
391
393 const SimplexId &tetId,
394 const SimplexId &triangleId0,
395 const SimplexId &triangleId1,
396 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections)
397 const;
398
400 const SimplexId &tetId,
401 const SimplexId &triangleId,
402 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
403 std::pair<double, double> &extremity0,
404 std::pair<double, double> &extremity1) const;
405
406 bool hasDuplicatedVertices(const double *p0,
407 const double *p1,
408 const double *p2) const;
409
410 int interpolateBasePoints(const std::array<double, 3> &p0,
411 const std::pair<double, double> &uv0,
412 const double &t0,
413 const std::array<double, 3> &p1,
414 const std::pair<double, double> &uv1,
415 const double &t1,
416 const double &t,
417 Vertex &v) const;
418
420 const SimplexId &source,
421 const SimplexId &destination,
422 const SimplexId &pivotVertexId,
423 const std::vector<std::pair<SimplexId, SimplexId>> &starNeighbors) const;
424
425 bool isEdgeFlippable(const SimplexId &edgeVertexId0,
426 const SimplexId &edgeVertexId1,
427 const SimplexId &otherVertexId0,
428 const SimplexId &otherVertexId1) const;
429
431 const SimplexId &tetId,
432 const SimplexId &triangleId,
433 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
434 const std::vector<std::vector<Vertex>> &tetNewVertices,
435 const SimplexId &vertexId0,
436 const SimplexId &vertexId1,
437 const SimplexId &vertexId2) const {
438
439 std::array<std::array<double, 3>, 3> points{};
440 for(int i = 0; i < 3; i++) {
441 SimplexId vertexId = vertexId0;
442 if(i == 1)
443 vertexId = vertexId1;
444 if(i == 2)
445 vertexId = vertexId2;
446
447 if(vertexId >= 0) {
448 for(int j = 0; j < 3; j++) {
449 points[i][j] = tetIntersections[tetId][triangleId].p_[vertexId][j];
450 }
451 } else {
452 for(int j = 0; j < 3; j++) {
453 points[i][j] = tetNewVertices[tetId][(-vertexId) - 1].p_[j];
454 }
455 }
456 }
457
459 points[0].data(), points[1].data(), points[2].data());
460 }
461
462 int mergeEdges(const double &distanceThreshold) const;
463
464 int mergeVertices(const double &distanceThreshold) const;
465
466 template <class dataTypeU, class dataTypeV>
467 inline int remeshIntersections() const;
468
469 int snapToBasePoint(const std::vector<std::vector<double>> &basePoints,
470 const std::vector<std::pair<double, double>> &uv,
471 const std::vector<double> &t,
472 Vertex &v) const;
473
474 int snapVertexBarycentrics() const;
475
477 const SimplexId &tetId,
478 const std::vector<std::pair<SimplexId, SimplexId>> &triangles) const;
479
480 bool pointSnapping_{false};
481
483 const void *uField_{}, *vField_{};
484 const float *pointSet_{};
486 const std::vector<std::vector<SimplexId>> *tetNeighbors_{};
487 std::array<SimplexId, 12> edgeImplicitEncoding_{
488 0, 1, 0, 2, 0, 3, 3, 1, 2, 1, 2, 3};
489
492
493 const std::vector<std::pair<std::pair<double, double>,
494 std::pair<double, double>>> *polygon_{};
495
496 std::vector<Vertex> *globalVertexList_{};
497 std::vector<std::vector<Vertex> *> polygonEdgeVertexLists_{};
498 std::vector<std::vector<Triangle> *> polygonEdgeTriangleLists_{};
499
500#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
501 RangeDrivenOctree octree_{};
502#endif
503 };
504} // namespace ttk
505
506#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
507template <class dataTypeU, class dataTypeV, typename triangulationType>
508inline int
509 ttk::FiberSurface::buildOctree(const triangulationType *const triangulation) {
510
511 if(!uField_)
512 return -1;
513 if(!vField_)
514 return -2;
515
516 if(octree_.empty()) {
517
518 octree_.setDebugLevel(debugLevel_);
519 octree_.setThreadNumber(threadNumber_);
520 if(!triangulation) {
521 octree_.setCellList(tetList_);
522 octree_.setCellNumber(tetNumber_);
523 octree_.setPointList(pointSet_);
524 octree_.setVertexNumber(pointNumber_);
525 }
526 octree_.setRange(uField_, vField_);
527
528 octree_.build<dataTypeU, dataTypeV>(triangulation);
529 }
530
531 return 0;
532}
533#endif
534
535template <class dataTypeU, class dataTypeV, typename triangulationType>
537 const SimplexId &tetId,
538 const SimplexId &localEdgeId0,
539 const double &t0,
540 const double &u0,
541 const double &v0,
542 const SimplexId &localEdgeId1,
543 const double &t1,
544 const double &u1,
545 const double &v1,
546 const SimplexId &localEdgeId2,
547 const double &t2,
548 const double &u2,
549 const double &v2,
550 std::array<std::array<double, 3>, 3> &basePoints,
551 std::array<std::pair<double, double>, 3> &basePointProjections,
552 std::array<double, 3> &basePointParameterization,
553 std::array<std::pair<SimplexId, SimplexId>, 3> &baseEdges,
554 const triangulationType *const triangulation) const {
555
556 for(int i = 0; i < 3; i++) {
557
558 SimplexId vertexId0 = 0, vertexId1 = 0;
559
560 switch(i) {
561
562 case 0:
563 if(!triangulation) {
564 vertexId0
565 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId0]];
566 vertexId1 = tetList_[5 * tetId + 1
567 + edgeImplicitEncoding_[2 * localEdgeId0 + 1]];
568 } else {
569 triangulation->getCellVertex(
570 tetId, edgeImplicitEncoding_[2 * localEdgeId0], vertexId0);
571 triangulation->getCellVertex(
572 tetId, edgeImplicitEncoding_[2 * localEdgeId0 + 1], vertexId1);
573 }
574 basePointProjections[i].first = u0;
575 basePointProjections[i].second = v0;
576 basePointParameterization[i] = t0;
577 break;
578
579 case 1:
580 if(!triangulation) {
581 vertexId0
582 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId1]];
583 vertexId1 = tetList_[5 * tetId + 1
584 + edgeImplicitEncoding_[2 * localEdgeId1 + 1]];
585 } else {
586 triangulation->getCellVertex(
587 tetId, edgeImplicitEncoding_[2 * localEdgeId1], vertexId0);
588 triangulation->getCellVertex(
589 tetId, edgeImplicitEncoding_[2 * localEdgeId1 + 1], vertexId1);
590 }
591 basePointProjections[i].first = u1;
592 basePointProjections[i].second = v1;
593 basePointParameterization[i] = t1;
594 break;
595
596 case 2:
597 if(!triangulation) {
598 vertexId0
599 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId2]];
600 vertexId1 = tetList_[5 * tetId + 1
601 + edgeImplicitEncoding_[2 * localEdgeId2 + 1]];
602 } else {
603 triangulation->getCellVertex(
604 tetId, edgeImplicitEncoding_[2 * localEdgeId2], vertexId0);
605 triangulation->getCellVertex(
606 tetId, edgeImplicitEncoding_[2 * localEdgeId2 + 1], vertexId1);
607 }
608 basePointProjections[i].first = u2;
609 basePointProjections[i].second = v2;
610 basePointParameterization[i] = t2;
611 break;
612 }
613
614 std::array<double, 2> baryCentrics{};
615 std::array<double, 2> p0{}, p1{}, p{};
616 p0[0] = ((const dataTypeU *)uField_)[vertexId0];
617 p0[1] = ((const dataTypeV *)vField_)[vertexId0];
618 p1[0] = ((const dataTypeU *)uField_)[vertexId1];
619 p1[1] = ((const dataTypeV *)vField_)[vertexId1];
620 p[0] = basePointProjections[i].first;
621 p[1] = basePointProjections[i].second;
623 p0.data(), p1.data(), p.data(), baryCentrics, 2);
624
625 std::array<float, 3> pA{}, pB{};
626 if(triangulation) {
627 triangulation->getVertexPoint(vertexId0, pA[0], pA[1], pA[2]);
628 triangulation->getVertexPoint(vertexId1, pB[0], pB[1], pB[2]);
629 }
630
631 for(int j = 0; j < 3; j++) {
632
633 double c0, c1;
634 if(!triangulation) {
635 c0 = pointSet_[3 * vertexId0 + j];
636 c1 = pointSet_[3 * vertexId1 + j];
637 } else {
638 c0 = pA[j];
639 c1 = pB[j];
640 }
641
642 basePoints[i][j] = baryCentrics[0] * c0 + baryCentrics[1] * c1;
643 }
644
645 if(vertexId0 < vertexId1) {
646 baseEdges[i] = std::pair<SimplexId, SimplexId>(vertexId0, vertexId1);
647 } else {
648 baseEdges[i] = std::pair<SimplexId, SimplexId>(vertexId1, vertexId0);
649 }
650 }
651
652 return 0;
653}
654
655template <class dataTypeU, class dataTypeV, typename triangulationType>
657 const SimplexId &polygonEdgeId,
658 const SimplexId &tetId,
659 const SimplexId &localEdgeId0,
660 const double &t0,
661 const double &u0,
662 const double &v0,
663 const SimplexId &localEdgeId1,
664 const double &t1,
665 const double &u1,
666 const double &v1,
667 const SimplexId &localEdgeId2,
668 const double &t2,
669 const double &u2,
670 const double &v2,
671 const triangulationType *const triangulation) const {
672
673 // that one's easy, make just one triangle
674 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
675
676 // alloc 1 more triangle
677 (*polygonEdgeTriangleLists_[polygonEdgeId])
678 .resize((*polygonEdgeTriangleLists_[polygonEdgeId]).size() + 1);
679 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().tetId_ = tetId;
680 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().caseId_ = 0;
681 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().polygonEdgeId_
682 = polygonEdgeId;
683
684 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[0] = vertexId;
685 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[1]
686 = vertexId + 1;
687 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[2]
688 = vertexId + 2;
689
690 // alloc 3 more vertices
691 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 3);
692 for(int i = 0; i < 3; i++) {
693 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
694 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
695 = false;
696 }
697
698 // get the vertex coordinates
699 for(int i = 0; i < 3; i++) {
700
701 SimplexId vertexId0 = 0, vertexId1 = 0;
702
703 switch(i) {
704 case 0:
705 if(!triangulation) {
706 vertexId0
707 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId0]];
708 vertexId1 = tetList_[5 * tetId + 1
709 + edgeImplicitEncoding_[2 * localEdgeId0 + 1]];
710 } else {
711 triangulation->getCellVertex(
712 tetId, edgeImplicitEncoding_[2 * localEdgeId0], vertexId0);
713 triangulation->getCellVertex(
714 tetId, edgeImplicitEncoding_[2 * localEdgeId0 + 1], vertexId1);
715 }
716 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u0;
717 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v0;
718 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t0;
719 break;
720
721 case 1:
722 if(!triangulation) {
723 vertexId0
724 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId1]];
725 vertexId1 = tetList_[5 * tetId + 1
726 + edgeImplicitEncoding_[2 * localEdgeId1 + 1]];
727 } else {
728 triangulation->getCellVertex(
729 tetId, edgeImplicitEncoding_[2 * localEdgeId1], vertexId0);
730 triangulation->getCellVertex(
731 tetId, edgeImplicitEncoding_[2 * localEdgeId1 + 1], vertexId1);
732 }
733 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u1;
734 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v1;
735 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t1;
736 break;
737
738 case 2:
739 if(!triangulation) {
740 vertexId0
741 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId2]];
742 vertexId1 = tetList_[5 * tetId + 1
743 + edgeImplicitEncoding_[2 * localEdgeId2 + 1]];
744 } else {
745 triangulation->getCellVertex(
746 tetId, edgeImplicitEncoding_[2 * localEdgeId2], vertexId0);
747 triangulation->getCellVertex(
748 tetId, edgeImplicitEncoding_[2 * localEdgeId2 + 1], vertexId1);
749 }
750 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u2;
751 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v2;
752 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t2;
753 break;
754 }
755
756 std::array<double, 2> baryCentrics{};
757 std::array<double, 2> p0{}, p1{}, p{};
758 p0[0] = ((const dataTypeU *)uField_)[vertexId0];
759 p0[1] = ((const dataTypeV *)vField_)[vertexId0];
760 p1[0] = ((const dataTypeU *)uField_)[vertexId1];
761 p1[1] = ((const dataTypeV *)vField_)[vertexId1];
762 p[0] = (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first;
763 p[1] = (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second;
765 p0.data(), p1.data(), p.data(), baryCentrics, 2);
766
767 float pA[3], pB[3];
768 if(triangulation) {
769 triangulation->getVertexPoint(vertexId0, pA[0], pA[1], pA[2]);
770 triangulation->getVertexPoint(vertexId1, pB[0], pB[1], pB[2]);
771 }
772
773 for(int j = 0; j < 3; j++) {
774
775 double c0, c1;
776 if(!triangulation) {
777 c0 = pointSet_[3 * vertexId0 + j];
778 c1 = pointSet_[3 * vertexId1 + j];
779 } else {
780 c0 = pA[j];
781 c1 = pB[j];
782 }
783
784 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
785 = baryCentrics[0] * c0 + baryCentrics[1] * c1;
786 }
787
788 if(vertexId0 < vertexId1)
789 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
790 = std::pair<SimplexId, SimplexId>(vertexId0, vertexId1);
791 else
792 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
793 = std::pair<SimplexId, SimplexId>(vertexId1, vertexId0);
794 }
795
796 // return the number of created vertices
797 return 3;
798}
799
800template <class dataTypeU, class dataTypeV, typename triangulationType>
802 const SimplexId &polygonEdgeId,
803 const SimplexId &tetId,
804 const SimplexId &localEdgeId0,
805 const double &t0,
806 const double &u0,
807 const double &v0,
808 const SimplexId &localEdgeId1,
809 const double &t1,
810 const double &u1,
811 const double &v1,
812 const SimplexId &localEdgeId2,
813 const double &t2,
814 const double &u2,
815 const double &v2,
816 const triangulationType *const triangulation) const {
817
818 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
819
820 // alloc 5 more vertices
821 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 5);
822 for(int i = 0; i < 5; i++) {
823 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
824 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
825 = false;
826 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
827 = std::pair<SimplexId, SimplexId>(-1, -1);
828 }
829
830 // alloc 3 more triangles
831 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
832 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 3);
833
834 for(int i = 0; i < 3; i++) {
835
836 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
837 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 1;
838 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
839 = polygonEdgeId;
840
841 switch(i) {
842 case 0:
843 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
844 .vertexIds_[0]
845 = vertexId;
846 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
847 .vertexIds_[1]
848 = vertexId + 1;
849 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
850 .vertexIds_[2]
851 = vertexId + 2;
852 break;
853 case 1:
854 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
855 .vertexIds_[0]
856 = vertexId + 1;
857 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
858 .vertexIds_[1]
859 = vertexId + 2;
860 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
861 .vertexIds_[2]
862 = vertexId + 3;
863 break;
864 case 2:
865 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
866 .vertexIds_[0]
867 = vertexId + 2;
868 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
869 .vertexIds_[1]
870 = vertexId + 3;
871 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
872 .vertexIds_[2]
873 = vertexId + 4;
874 break;
875 }
876 }
877
878 // compute the base triangle vertices like in case 1
879 std::array<std::array<double, 3>, 3> basePoints{};
880 std::array<std::pair<double, double>, 3> basePointProjections{};
881 std::array<double, 3> basePointParameterization{};
882 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
883
884 computeBaseTriangle<dataTypeU, dataTypeV>(
885 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
886 u2, v2, basePoints, basePointProjections, basePointParameterization,
887 baseEdges, triangulation);
888
889 // find the pivot vertex for this case
890 SimplexId pivotVertexId = -1;
891
892 if((t0 >= 0) && (t0 <= 1)) {
893 pivotVertexId = 0;
894 }
895 if((t1 >= 0) && (t1 <= 1)) {
896 pivotVertexId = 1;
897 }
898 if((t2 >= 0) && (t2 <= 1)) {
899 pivotVertexId = 2;
900 }
901
902 // now get the vertex coordinates
903 for(int i = 0; i < 5; i++) {
904
905 SimplexId vertexId0 = -1, vertexId1 = -1;
906 double t{};
907
908 if(!i) {
909 // just take the pivot vertex
910 for(int j = 0; j < 3; j++) {
911 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
912 = basePoints[pivotVertexId][j];
913 }
914
915 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
916 = basePointParameterization[pivotVertexId];
917 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
918 = basePointProjections[pivotVertexId];
919 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
920 = baseEdges[pivotVertexId];
921 } else {
922
923 switch(i) {
924
925 case 1:
926 // interpolation between pivotVertexId and (pivotVertexId-1)%3
927 // if pivot is positive: interpolate at 1
928 // if not interpolate at 0
929 vertexId0 = pivotVertexId;
930 vertexId1 = (pivotVertexId + 2) % 3;
931 if(basePointParameterization[(pivotVertexId + 2) % 3] > 1) {
932 t = 1;
933 } else {
934 t = 0;
935 }
936 break;
937
938 case 2:
939 // interpolation between pivotVertexId and (pivotVertexId+1)%3
940 // if pivot is positive: interpolate at 1
941 // if not interpolate at 0
942 vertexId0 = pivotVertexId;
943 vertexId1 = (pivotVertexId + 1) % 3;
944 if(basePointParameterization[(pivotVertexId + 1) % 3] > 1) {
945 t = 1;
946 } else {
947 t = 0;
948 }
949 break;
950
951 case 3:
952 vertexId0 = (pivotVertexId + 2) % 3;
953 vertexId1 = (pivotVertexId + 1) % 3;
954 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
955 t = 0;
956 } else {
957 t = 1;
958 }
959 break;
960
961 case 4:
962 vertexId0 = (pivotVertexId + 2) % 3;
963 vertexId1 = (pivotVertexId + 1) % 3;
964 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
965 t = 1;
966 } else {
967 t = 0;
968 }
969 break;
970 }
971
972 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
973
974 interpolateBasePoints(
975 basePoints[vertexId0], basePointProjections[vertexId0],
976 basePointParameterization[vertexId0], basePoints[vertexId1],
977 basePointProjections[vertexId1], basePointParameterization[vertexId1],
978 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
979 // snapToBasePoint(
980 // basePoints, basePointProjections, basePointParameterization,
981 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
982 }
983 }
984
985 // return the number of created vertices
986 return 5;
987}
988
989template <class dataTypeU, class dataTypeV, typename triangulationType>
991 const SimplexId &polygonEdgeId,
992 const SimplexId &tetId,
993 const SimplexId &localEdgeId0,
994 const double &t0,
995 const double &u0,
996 const double &v0,
997 const SimplexId &localEdgeId1,
998 const double &t1,
999 const double &u1,
1000 const double &v1,
1001 const SimplexId &localEdgeId2,
1002 const double &t2,
1003 const double &u2,
1004 const double &v2,
1005 const triangulationType *const triangulation) const {
1006
1007 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1008
1009 // alloc 4 more vertices
1010 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1011 for(int i = 0; i < 4; i++) {
1012 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1013 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1014 = false;
1015 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1016 = std::pair<SimplexId, SimplexId>(-1, -1);
1017 }
1018
1019 // alloc 2 more triangles
1020 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1021 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1022
1023 for(int i = 0; i < 2; i++) {
1024
1025 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1026 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 2;
1027 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1028 = polygonEdgeId;
1029
1030 if(!i) {
1031 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1032 = vertexId;
1033 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1034 = vertexId + 1;
1035 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1036 = vertexId + 2;
1037 } else {
1038 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1039 = vertexId + 1;
1040 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1041 = vertexId + 3;
1042 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1043 = vertexId + 2;
1044 }
1045 }
1046
1047 // compute the base triangle vertices like in case 1
1048 std::array<std::array<double, 3>, 3> basePoints{};
1049 std::array<std::pair<double, double>, 3> basePointProjections{};
1050 std::array<double, 3> basePointParameterization{};
1051 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1052
1053 computeBaseTriangle<dataTypeU, dataTypeV>(
1054 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1055 u2, v2, basePoints, basePointProjections, basePointParameterization,
1056 baseEdges, triangulation);
1057
1058 // find the pivot for this case
1059 bool isPivotPositive = false;
1060 SimplexId pivotVertexId = -1;
1061
1062 if(((t0 < 0) && ((t1 < 0) || (t2 < 0)))
1063 || ((t1 < 0) && ((t0 < 0) || (t2 < 0)))
1064 || ((t2 < 0) && ((t1 < 0) || (t0 < 0)))) {
1065 isPivotPositive = true;
1066 }
1067 if(isPivotPositive) {
1068 if(t0 >= 1)
1069 pivotVertexId = 0;
1070 } else {
1071 if(t0 <= 0)
1072 pivotVertexId = 0;
1073 }
1074 if(isPivotPositive) {
1075 if(t1 >= 1)
1076 pivotVertexId = 1;
1077 } else {
1078 if(t1 <= 0)
1079 pivotVertexId = 1;
1080 }
1081 if(isPivotPositive) {
1082 if(t2 >= 1)
1083 pivotVertexId = 2;
1084 } else {
1085 if(t2 <= 0)
1086 pivotVertexId = 2;
1087 }
1088
1089 // now get the vertex coordinates
1090 for(int i = 0; i < 4; i++) {
1091
1092 SimplexId vertexId0 = -1, vertexId1 = -1;
1093 double t{};
1094
1095 switch(i) {
1096
1097 case 0:
1098 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1099 // if pivot is positive: interpolate at 1
1100 // if not interpolate at 0
1101 vertexId0 = pivotVertexId;
1102 vertexId1 = (pivotVertexId + 2) % 3;
1103 if(isPivotPositive)
1104 t = 1;
1105 else
1106 t = 0;
1107 break;
1108
1109 case 1:
1110 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1111 // if pivot is positive: interpolate at 1
1112 // if not interpolate at 0
1113 vertexId0 = pivotVertexId;
1114 vertexId1 = (pivotVertexId + 1) % 3;
1115 if(isPivotPositive)
1116 t = 1;
1117 else
1118 t = 0;
1119 break;
1120
1121 case 2:
1122 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1123 // if pivot is positive: interpolate at 0
1124 // if not interpolate at 1
1125 vertexId0 = pivotVertexId;
1126 vertexId1 = (pivotVertexId + 2) % 3;
1127 if(isPivotPositive)
1128 t = 0;
1129 else
1130 t = 1;
1131 break;
1132
1133 case 3:
1134 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1135 // if pivot is positive: interpolate at 0
1136 // if not interpolate at 1
1137 vertexId0 = pivotVertexId;
1138 vertexId1 = (pivotVertexId + 1) % 3;
1139 if(isPivotPositive)
1140 t = 0;
1141 else
1142 t = 1;
1143 break;
1144 }
1145
1146 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1147
1148 interpolateBasePoints(
1149 basePoints[vertexId0], basePointProjections[vertexId0],
1150 basePointParameterization[vertexId0], basePoints[vertexId1],
1151 basePointProjections[vertexId1], basePointParameterization[vertexId1], t,
1152 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1153 // snapToBasePoint(
1154 // basePoints, basePointProjections, basePointParameterization,
1155 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1156 }
1157
1158 // return the number of created vertices
1159 return 4;
1160}
1161
1162template <class dataTypeU, class dataTypeV, typename triangulationType>
1164 const SimplexId &polygonEdgeId,
1165 const SimplexId &tetId,
1166 const SimplexId &localEdgeId0,
1167 const double &t0,
1168 const double &u0,
1169 const double &v0,
1170 const SimplexId &localEdgeId1,
1171 const double &t1,
1172 const double &u1,
1173 const double &v1,
1174 const SimplexId &localEdgeId2,
1175 const double &t2,
1176 const double &u2,
1177 const double &v2,
1178 const triangulationType *const triangulation) const {
1179
1180 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1181
1182 // alloc 3 more vertices
1183 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 3);
1184 for(int i = 0; i < 3; i++) {
1185 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1186 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1187 = false;
1188 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1189 = std::pair<SimplexId, SimplexId>(-1, -1);
1190 }
1191
1192 // alloc 1 more triangle
1193 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1194 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 1);
1195
1196 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].tetId_ = tetId;
1197 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].caseId_ = 3;
1198 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().polygonEdgeId_
1199 = polygonEdgeId;
1200
1201 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[0]
1202 = vertexId;
1203 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[1]
1204 = vertexId + 1;
1205 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[2]
1206 = vertexId + 2;
1207
1208 // compute the base triangle vertices like in case 1
1209 std::array<std::array<double, 3>, 3> basePoints{};
1210 std::array<std::pair<double, double>, 3> basePointProjections{};
1211 std::array<double, 3> basePointParameterization{};
1212 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1213
1214 computeBaseTriangle<dataTypeU, dataTypeV>(
1215 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1216 u2, v2, basePoints, basePointProjections, basePointParameterization,
1217 baseEdges, triangulation);
1218
1219 // now find the pivot
1220 bool isPivotPositive = false;
1221 SimplexId pivotVertexId = -1;
1222
1223 if((t0 <= 1) && (t0 >= 0)) {
1224 pivotVertexId = 0;
1225 } else {
1226 if(t0 < 0)
1227 isPivotPositive = true;
1228 }
1229 if((t1 <= 1) && (t1 >= 0)) {
1230 pivotVertexId = 1;
1231 } else {
1232 if(t1 < 0)
1233 isPivotPositive = true;
1234 }
1235 if((t2 <= 1) && (t2 >= 0)) {
1236 pivotVertexId = 2;
1237 } else {
1238 if(t2 < 0)
1239 isPivotPositive = true;
1240 }
1241
1242 // now get the vertex coordinates
1243 for(int i = 0; i < 3; i++) {
1244
1245 SimplexId vertexId0 = -1, vertexId1 = -1;
1246 double t{};
1247
1248 if(!i) {
1249 // special case of the pivot vertex
1250 for(int j = 0; j < 3; j++) {
1251 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
1252 = basePoints[pivotVertexId][j];
1253 }
1254
1255 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
1256 = basePointParameterization[pivotVertexId];
1257 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
1258 = basePointProjections[pivotVertexId];
1259 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
1260 = baseEdges[pivotVertexId];
1261 } else {
1262 if(i == 1) {
1263 // interpolation between pivotVertexId and pivotVertexId+1
1264 // if pivot is positive, interpolate at value 0
1265 // else interpolate at value 1
1266 vertexId0 = pivotVertexId;
1267 vertexId1 = (pivotVertexId + 1) % 3;
1268 if(isPivotPositive)
1269 t = 0;
1270 else
1271 t = 1;
1272 }
1273 if(i == 2) {
1274 // interpolation between pivotVertexId and pivotVertexId-1
1275 // if pivot is positive, interpolate at value 0
1276 // else interpolate at value 1
1277 vertexId0 = pivotVertexId;
1278 vertexId1 = (pivotVertexId + 2) % 3;
1279 if(isPivotPositive)
1280 t = 0;
1281 else
1282 t = 1;
1283 }
1284
1285 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1286
1287 interpolateBasePoints(
1288 basePoints[vertexId0], basePointProjections[vertexId0],
1289 basePointParameterization[vertexId0], basePoints[vertexId1],
1290 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1291 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1292 // snapToBasePoint(
1293 // basePoints, basePointProjections, basePointParameterization,
1294 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1295 }
1296 }
1297
1298 // return the number of created vertices
1299 return 3;
1300}
1301
1302template <class dataTypeU, class dataTypeV, typename triangulationType>
1304 const SimplexId &polygonEdgeId,
1305 const SimplexId &tetId,
1306 const SimplexId &localEdgeId0,
1307 const double &t0,
1308 const double &u0,
1309 const double &v0,
1310 const SimplexId &localEdgeId1,
1311 const double &t1,
1312 const double &u1,
1313 const double &v1,
1314 const SimplexId &localEdgeId2,
1315 const double &t2,
1316 const double &u2,
1317 const double &v2,
1318 const triangulationType *const triangulation) const {
1319
1320 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1321
1322 // alloc 4 more vertices
1323 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1324 for(int i = 0; i < 4; i++) {
1325 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1326 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1327 = false;
1328 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1329 = std::pair<SimplexId, SimplexId>(-1, -1);
1330 }
1331
1332 // alloc 2 more triangles
1333 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1334 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1335
1336 for(int i = 0; i < 2; i++) {
1337
1338 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1339 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 4;
1340 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1341 = polygonEdgeId;
1342
1343 if(!i) {
1344 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1345 = vertexId;
1346 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1347 = vertexId + 1;
1348 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1349 = vertexId + 2;
1350 } else {
1351 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1352 = vertexId + 1;
1353 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1354 = vertexId + 3;
1355 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1356 = vertexId + 2;
1357 }
1358 }
1359
1360 // compute the base triangle vertices like in case 1
1361 std::array<std::array<double, 3>, 3> basePoints{};
1362 std::array<std::pair<double, double>, 3> basePointProjections{};
1363 std::array<double, 3> basePointParameterization{};
1364 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1365
1366 computeBaseTriangle<dataTypeU, dataTypeV>(
1367 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1368 u2, v2, basePoints, basePointProjections, basePointParameterization,
1369 baseEdges, triangulation);
1370
1371 // find the pivot vertex for this case
1372 bool isPivotPositive = false;
1373 SimplexId pivotVertexId = -1;
1374
1375 if(t0 > 1) {
1376 pivotVertexId = 0;
1377 isPivotPositive = true;
1378 } else if(t0 < 0) {
1379 pivotVertexId = 0;
1380 isPivotPositive = false;
1381 }
1382
1383 if(t1 > 1) {
1384 pivotVertexId = 1;
1385 isPivotPositive = true;
1386 } else if(t1 < 0) {
1387 pivotVertexId = 1;
1388 isPivotPositive = false;
1389 }
1390 if(t2 > 1) {
1391 pivotVertexId = 2;
1392 isPivotPositive = true;
1393 } else if(t2 < 0) {
1394 pivotVertexId = 2;
1395 isPivotPositive = false;
1396 }
1397
1398 // now get the vertex coordinates depending on the case
1399 for(int i = 0; i < 4; i++) {
1400
1401 SimplexId vertexId0 = -1, vertexId1 = -1;
1402 double t{};
1403
1404 if(i < 2) {
1405 if(!i) {
1406 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1407 // if pivot is positive: interpolate at 1
1408 // if not interpolate at 0
1409 vertexId0 = pivotVertexId;
1410 vertexId1 = (pivotVertexId + 2) % 3;
1411 if(isPivotPositive)
1412 t = 1;
1413 else
1414 t = 0;
1415 } else {
1416 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1417 // if pivot is positive: interpolate at 1
1418 // if not interpolate at 0
1419 vertexId0 = pivotVertexId;
1420 vertexId1 = (pivotVertexId + 1) % 3;
1421 if(isPivotPositive)
1422 t = 1;
1423 else
1424 t = 0;
1425 }
1426
1427 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1428
1429 interpolateBasePoints(
1430 basePoints[vertexId0], basePointProjections[vertexId0],
1431 basePointParameterization[vertexId0], basePoints[vertexId1],
1432 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1433 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1434 // snapToBasePoint(
1435 // basePoints, basePointProjections, basePointParameterization,
1436 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1437
1438 } else {
1439 if(i == 2) {
1440 // take (pivotVertexId-1)%3
1441 for(int j = 0; j < 3; j++) {
1442 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1443 = basePoints[(pivotVertexId + 2) % 3][j];
1444 }
1445
1446 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1447 = basePointParameterization[(pivotVertexId + 2) % 3];
1448 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1449 = basePointProjections[(pivotVertexId + 2) % 3];
1450 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1451 = baseEdges[(pivotVertexId + 2) % 3];
1452 } else {
1453 // take (pivtoVertexId+1)%3
1454 for(int j = 0; j < 3; j++) {
1455 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1456 = basePoints[(pivotVertexId + 1) % 3][j];
1457 }
1458 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1459 = basePointParameterization[(pivotVertexId + 1) % 3];
1460 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1461 = basePointProjections[(pivotVertexId + 1) % 3];
1462 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1463 = baseEdges[(pivotVertexId + 1) % 3];
1464 }
1465 }
1466 }
1467
1468 // return the number of created vertices
1469 return 4;
1470}
1471
1472template <class dataTypeU, class dataTypeV, typename triangulationType>
1474 const std::pair<double, double> &rangePoint0,
1475 const std::pair<double, double> &rangePoint1,
1476 const std::vector<SimplexId> &seedTetList,
1477 const triangulationType *const triangulation,
1478 const SimplexId &polygonEdgeId) const {
1479
1480#ifndef TTK_ENABLE_KAMIKAZE
1481 if(!uField_)
1482 return -4;
1483 if(!vField_)
1484 return -5;
1485 if(!triangulation)
1486 return -6;
1487 if(!polygonEdgeNumber_)
1488 return -7;
1489 if(!globalVertexList_)
1490 return -8;
1491#endif
1492
1493 std::vector<bool> visitedTets(triangulation->getNumberOfCells(), false);
1494 std::queue<SimplexId> tetQueue;
1495
1496 // init the queue
1497 for(SimplexId i = 0; i < (SimplexId)seedTetList.size(); i++) {
1498 tetQueue.push(seedTetList[i]);
1499 }
1500
1501 SimplexId createdVertices = 0;
1502
1503 do {
1504
1505 SimplexId tetId = tetQueue.front();
1506 tetQueue.pop();
1507
1508 if(!visitedTets[tetId]) {
1509
1510 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1511 tetId, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1512
1513 if(createdVertices) {
1514 // only propagate if we created a triangle
1515 SimplexId tetNeighborNumber
1516 = triangulation->getCellNeighborNumber(tetId);
1517
1518 for(SimplexId i = 0; i < tetNeighborNumber; i++) {
1519
1520 SimplexId neighborId = -1;
1521 triangulation->getCellNeighbor(tetId, i, neighborId);
1522
1523 if(!visitedTets[neighborId])
1524 tetQueue.push(neighborId);
1525 }
1526 }
1527
1528 visitedTets[tetId] = true;
1529 }
1530
1531 } while(tetQueue.size());
1532
1533 return 0;
1534}
1535
1536template <class dataTypeU, class dataTypeV>
1538 const std::vector<
1539 std::pair<std::pair<double, double>, std::pair<double, double>>> &edgeList,
1540 const std::vector<SimplexId> &seedTetList,
1541 const std::vector<SimplexId> *edgeIdList) const {
1542
1543#ifndef TTK_ENABLE_KAMIKAZE
1544 if(!tetNeighbors_)
1545 return -1;
1546 if(!tetNumber_)
1547 return -2;
1548 if(!tetList_)
1549 return -3;
1550 if(!uField_)
1551 return -4;
1552 if(!vField_)
1553 return -5;
1554 if(!pointSet_)
1555 return -6;
1556 if(!polygonEdgeNumber_)
1557 return -7;
1558 if(!globalVertexList_)
1559 return -8;
1560#endif
1561
1562 std::vector<bool> visitedTets(tetNumber_, false);
1563 std::queue<SimplexId> tetQueue;
1564
1565 // init the queue
1566 for(SimplexId i = 0; i < (SimplexId)seedTetList.size(); i++) {
1567 tetQueue.push(seedTetList[i]);
1568 }
1569
1570 SimplexId createdVertices = 0;
1571
1572 do {
1573
1574 SimplexId tetId = tetQueue.front();
1575 tetQueue.pop();
1576
1577 if(!visitedTets[tetId]) {
1578
1579 std::vector<std::vector<SimplexId>> threadedTetQueue(edgeList.size());
1580#ifdef TTK_ENABLE_OPENMP
1581#pragma omp parallel for num_threads(threadNumber_)
1582#endif
1583 for(SimplexId i = 0; i < (SimplexId)edgeList.size(); i++) {
1584
1585 SimplexId polygonEdgeId = 0;
1586
1587 if(edgeIdList) {
1588 polygonEdgeId = (*edgeIdList)[i];
1589 }
1590
1591 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1592 tetId, edgeList[i].first, edgeList[i].second, polygonEdgeId);
1593
1594 if(createdVertices) {
1595 // only propagate if we created a triangle
1596 for(SimplexId j = 0; j < (SimplexId)(*tetNeighbors_)[tetId].size();
1597 j++) {
1598 if(!visitedTets[(*tetNeighbors_)[tetId][j]]) {
1599 threadedTetQueue[i].push_back((*tetNeighbors_)[tetId][j]);
1600 }
1601 }
1602 }
1603 }
1604
1605 visitedTets[tetId] = true;
1606
1607 for(SimplexId i = 0; i < (SimplexId)threadedTetQueue.size(); i++) {
1608 for(SimplexId j = 0; j < (SimplexId)threadedTetQueue[i].size(); j++) {
1609 tetQueue.push(threadedTetQueue[i][j]);
1610 }
1611 }
1612 }
1613
1614 } while(tetQueue.size());
1615
1616 return 0;
1617}
1618
1619template <class dataTypeU, class dataTypeV, typename triangulationType>
1621 const std::pair<double, double> &rangePoint0,
1622 const std::pair<double, double> &rangePoint1,
1623 const triangulationType *const triangulation,
1624 const SimplexId &polygonEdgeId) const {
1625
1626#ifndef TTK_ENABLE_KAMIKAZE
1627 if((!tetNumber_) && (!triangulation))
1628 return -1;
1629 if((!tetList_) && (!triangulation))
1630 return -2;
1631 if(!uField_)
1632 return -3;
1633 if(!vField_)
1634 return -4;
1635 if((!pointSet_) && (!triangulation))
1636 return -5;
1637 if(!polygonEdgeNumber_)
1638 return -6;
1639 if(!globalVertexList_)
1640 return -7;
1641#endif
1642
1643 SimplexId tetNumber = tetNumber_;
1644
1645 if(triangulation) {
1646 tetNumber = triangulation->getNumberOfCells();
1647 }
1648
1649#ifdef TTK_ENABLE_OPENMP
1650#pragma omp parallel for num_threads(threadNumber_)
1651#endif
1652 for(SimplexId i = 0; i < tetNumber; i++) {
1653
1654 processTetrahedron<dataTypeU, dataTypeV>(
1655 i, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1656 }
1657
1658 return 0;
1659}
1660
1661template <class dataTypeU, class dataTypeV, typename triangulationType>
1663 const triangulationType *const triangulation) {
1664
1665#ifndef TTK_ENABLE_KAMIKAZE
1666 if((!tetNumber_) && (!triangulation))
1667 return -1;
1668 if((!tetList_) && (!triangulation))
1669 return -2;
1670 if(!uField_)
1671 return -3;
1672 if(!vField_)
1673 return -4;
1674 if((!pointSet_) && (!triangulation))
1675 return -5;
1676 if(!polygon_)
1677 return -6;
1678 if(polygonEdgeNumber_ != (SimplexId)polygon_->size())
1679 return -7;
1680 if(!globalVertexList_)
1681 return -8;
1682#endif
1683
1684 Timer t;
1685
1686#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1687 if(!octree_.empty()) {
1688
1689#ifdef TTK_ENABLE_OPENMP
1690#pragma omp parallel for num_threads(threadNumber_)
1691#endif
1692 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1693
1694 computeSurfaceWithOctree<dataTypeU, dataTypeV>(
1695 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1696 }
1697 } else {
1698 // regular extraction (the octree has not been computed)
1699#ifdef TTK_ENABLE_OPENMP
1700#pragma omp parallel for num_threads(threadNumber_)
1701#endif
1702 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1703 computeSurface<dataTypeU, dataTypeV>(
1704 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1705 }
1706 }
1707
1708#else
1709#ifdef TTK_ENABLE_OPENMP
1710#pragma omp parallel for num_threads(threadNumber_)
1711#endif
1712 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1713
1714 computeSurface<dataTypeU, dataTypeV>(
1715 (*polygon_)[i].first, (*polygon_)[i].second, i);
1716 }
1717#endif
1718
1719 finalize<dataTypeU, dataTypeV>(pointSnapping_, false, false, false);
1720
1721 this->printMsg("Extracted", 1.0, t.getElapsedTime(), this->threadNumber_);
1722
1723 return 0;
1724}
1725
1726#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1727template <class dataTypeU, class dataTypeV, typename triangulationType>
1728inline int ttk::FiberSurface::computeSurfaceWithOctree(
1729 const std::pair<double, double> &rangePoint0,
1730 const std::pair<double, double> &rangePoint1,
1731 const triangulationType *const triangulation,
1732 const SimplexId &polygonEdgeId) const {
1733
1734#ifndef TTK_ENABLE_KAMIKAZE
1735 if((!tetNumber_) && (!triangulation))
1736 return -1;
1737 if((!tetList_) && (!triangulation))
1738 return -2;
1739 if(!uField_)
1740 return -3;
1741 if(!vField_)
1742 return -4;
1743 if((!pointSet_) && (!triangulation))
1744 return -5;
1745 if(!polygonEdgeNumber_)
1746 return -6;
1747 if(!globalVertexList_)
1748 return -7;
1749#endif
1750
1751 std::vector<SimplexId> tetList;
1752 octree_.rangeSegmentQuery(rangePoint0, rangePoint1, tetList);
1753
1754#ifdef TTK_ENABLE_OPENMP
1755#pragma omp parallel for num_threads(threadNumber_)
1756#endif
1757 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
1758 processTetrahedron<dataTypeU, dataTypeV>(
1759 tetList[i], rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1760 }
1761
1762 return 0;
1763}
1764#endif
1765
1766template <class dataTypeU, class dataTypeV>
1767int ttk::FiberSurface::finalize(const bool &mergeDuplicatedVertices,
1768 const bool &removeSmallEdges,
1769 const bool &edgeFlips,
1770 const bool &intersectionRemesh) {
1771
1772 // make only one vertex list
1773 SimplexId fiberSurfaceVertexNumber = 0;
1774 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1775 fiberSurfaceVertexNumber += (*polygonEdgeVertexLists_[i]).size();
1776 }
1777
1778 (*globalVertexList_).resize(fiberSurfaceVertexNumber);
1779 fiberSurfaceVertexNumber = 0;
1780 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1781 for(SimplexId j = 0; j < (SimplexId)polygonEdgeVertexLists_[i]->size();
1782 j++) {
1783 (*polygonEdgeVertexLists_[i])[j].polygonEdgeId_ = i;
1784 (*polygonEdgeVertexLists_[i])[j].localId_ = j;
1785 (*polygonEdgeVertexLists_[i])[j].globalId_ = fiberSurfaceVertexNumber;
1786 (*globalVertexList_)[fiberSurfaceVertexNumber]
1787 = (*polygonEdgeVertexLists_[i])[j];
1788 fiberSurfaceVertexNumber++;
1789 }
1790 }
1791 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1792 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1793 j++) {
1794 for(int k = 0; k < 3; k++) {
1795 (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k]
1796 = (*polygonEdgeVertexLists_[i])[(*polygonEdgeTriangleLists_[i])[j]
1797 .vertexIds_[k]]
1798 .globalId_;
1799 }
1800 }
1801 }
1802
1803 // NOTE:
1804 // 1) in a first step, make everything work in a POST process.
1805 // this is what really matters for tet gen.
1806 // 2) probably come back and apply in a pre-process (more degenerate
1807 // intersection cases).
1808 if(intersectionRemesh) {
1809 remeshIntersections<dataTypeU, dataTypeV>();
1810 }
1811
1812 if((mergeDuplicatedVertices) || (removeSmallEdges)) {
1813 // we need to have a complex to perform the edge collapses
1814 mergeVertices(pointSnappingThreshold_);
1815 }
1816
1817 if(edgeFlips)
1818 flipEdges();
1819
1820 if(removeSmallEdges)
1821 mergeEdges(edgeCollapseThreshold_);
1822
1823 // now we can release the memory for the threaded vertices
1824 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1825 polygonEdgeVertexLists_[i]->clear();
1826 }
1827
1828 return 0;
1829}
1830
1831template <class dataTypeU, class dataTypeV, typename triangulationType>
1833 const SimplexId &tetId,
1834 const std::pair<double, double> &rangePoint0,
1835 const std::pair<double, double> &rangePoint1,
1836 const triangulationType *const triangulation,
1837 const SimplexId &polygonEdgeId) const {
1838
1839 double rangeEdge[2];
1840 rangeEdge[0] = rangePoint0.first - rangePoint1.first;
1841 rangeEdge[1] = rangePoint0.second - rangePoint1.second;
1842
1843 double rangeNormal[2];
1844 rangeNormal[0] = -rangeEdge[1];
1845 rangeNormal[1] = rangeEdge[0];
1846
1847 const double prec_dbl = Geometry::powInt(10.0, -DBL_DIG);
1848
1849 // 1. compute the distance to the range line carrying the saddleEdge
1850 SimplexId upperNumber = 0;
1851 SimplexId lowerNumber = 0;
1852 SimplexId equalVertexLocalId = -1;
1853 double d[4];
1854 for(int i = 0; i < 4; i++) {
1855
1856 SimplexId vertexId = 0;
1857 if(!triangulation) {
1858 vertexId = tetList_[5 * tetId + 1 + i];
1859 } else {
1860 triangulation->getCellVertex(tetId, i, vertexId);
1861 }
1862
1863 double projectedVertex[2];
1864 projectedVertex[0] = ((const dataTypeU *)uField_)[vertexId];
1865 projectedVertex[1] = ((const dataTypeV *)vField_)[vertexId];
1866
1867 double vertexRangeEdge[2];
1868 vertexRangeEdge[0] = projectedVertex[0] - rangePoint0.first;
1869 vertexRangeEdge[1] = projectedVertex[1] - rangePoint0.second;
1870
1871 d[i] = vertexRangeEdge[0] * rangeNormal[0]
1872 + vertexRangeEdge[1] * rangeNormal[1];
1873
1874 if(fabs(d[i]) < prec_dbl)
1875 d[i] = 0;
1876
1877 if(d[i] > 0)
1878 upperNumber++;
1879 if(d[i] < 0)
1880 lowerNumber++;
1881
1882 if(d[i] == 0) {
1883 equalVertexLocalId = i;
1884 }
1885 }
1886
1887 // 2. compute the base triangle(s)
1888 if(!((upperNumber == 0) || (lowerNumber == 0))) {
1889
1890 // the fiber surface is passing through this tetrahedron.
1891 std::array<SimplexId, 2> triangleEdgeNumbers{};
1892 std::array<std::array<SimplexId, 3>, 2> triangleEdges{
1893 {{-1, -1, -1}, {-1, -1, -1}}};
1894
1895 // implicit edge encoding
1896 // 0: O-1
1897 // 1: 0-2
1898 // 2: 0-3
1899 // 3: 3-1 [order!]
1900 // 4: 2-1 [order!]
1901 // 5: 2-3
1902 SimplexId edgeCounter = 0;
1903 for(int i = 0; i < 4; i++) {
1904
1905 SimplexId jStart = i + 1;
1906 SimplexId jEnd = 4;
1907 SimplexId jStep = 1;
1908
1909 if(i == 1) {
1910 // special ordering here
1911 // (to facilitate the creation of valid base triangles)
1912 // any two consecutive edges shall share a vertex
1913 jStart = 3;
1914 jEnd = i;
1915 jStep = -1;
1916 }
1917
1918 for(SimplexId j = jStart; j != jEnd; j += jStep) {
1919
1920 if(((d[i] > 0) && (d[j] < 0)) || ((d[i] < 0) && (d[j] > 0))) {
1921
1922 // the edge is crossed by a base triangle
1923 if(triangleEdgeNumbers[0] == 3) {
1924 triangleEdges[1][triangleEdgeNumbers[1]] = edgeCounter;
1925 triangleEdgeNumbers[1]++;
1926 } else {
1927 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1928 triangleEdgeNumbers[0]++;
1929 }
1930 }
1931
1932 if((d[i] == 0) && (d[j] == 0)) {
1933 // special case of a seed tet containing the jacobi edge.
1934 // the entire edge is on the fiber surface.
1935 // let's put this edge twice.
1936 // NOTE: in such a case, we're producing only one triangle.
1937 // so we just need to add that to the first triangle.
1938 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1939 triangleEdgeNumbers[0]++;
1940
1941 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1942 triangleEdgeNumbers[0]++;
1943 }
1944
1945 edgeCounter++;
1946 }
1947 }
1948
1949 // post-process in the case of a second base triangle
1950 if(triangleEdges[1][0] != -1) {
1951 if(triangleEdges[1][1] == -1) {
1952 // we need to consider the edges of the first triangle which share
1953 // a vertex with the second triangle's edge.
1954 // given the ordering, the following edge pairs are forbidden:
1955 // (opposite edges of the tetrahedron)
1956 // 0/5
1957 // 1/3
1958 // 2/4
1959 SimplexId forbiddenEdge = -1;
1960 switch(triangleEdges[1][0]) {
1961 case 0:
1962 forbiddenEdge = 5;
1963 break;
1964 case 1:
1965 forbiddenEdge = 3;
1966 break;
1967 case 2:
1968 forbiddenEdge = 4;
1969 break;
1970 case 3:
1971 forbiddenEdge = 1;
1972 break;
1973 case 4:
1974 forbiddenEdge = 2;
1975 break;
1976 case 5:
1977 forbiddenEdge = 0;
1978 break;
1979 }
1980 for(SimplexId i = 0; i < (SimplexId)triangleEdges[0].size(); i++) {
1981 if(triangleEdges[0][i] != forbiddenEdge) {
1982 if(triangleEdges[1][1] != -1) {
1983 triangleEdges[1][2] = triangleEdges[0][i];
1984 break;
1985 } else {
1986 triangleEdges[1][1] = triangleEdges[0][i];
1987 }
1988 }
1989 }
1990 }
1991 }
1992
1993 // post-process in case of exactly one vertex on the jacobi edge
1994 for(int i = 0; i < 2; i++) {
1995 if(triangleEdges[i][0] != -1) {
1996 // the jacobi vertex has to be the last one
1997 if(triangleEdges[i][2] == -1) {
1998 // add whatever edge connected to equalVertexLocalId
1999 switch(equalVertexLocalId) {
2000 case 0:
2001 case 1:
2002 triangleEdges[i][2] = 0;
2003 break;
2004 case 2:
2005 case 3:
2006 triangleEdges[i][2] = 5;
2007 break;
2008 }
2009 }
2010 }
2011 }
2012
2013 // 3. crop the resulting triangles to the saddleEdge
2014 // for each edge recorded previously, get the u,v coordinates of the
2015 // intersection point.
2016 // take its barycentric coordinates on the saddle edge.
2017 // see if it lies in it (in between 0 and 1)
2018 // figure 7 of the paper
2019 double d0, d1;
2020 std::pair<double, double> uv0, uv1;
2021 std::array<std::pair<double, double>, 3> uv{};
2022 std::array<double, 3> t{};
2023
2024 SimplexId createdVertices = 0;
2025
2026 for(int i = 0; i < 2; i++) {
2027 if(triangleEdges[i][0] != -1) {
2028 // this is a valid triangle, let's go ahead.
2029
2030 SimplexId lowerVertexNumber = 0;
2031 SimplexId upperVertexNumber = 0;
2032 SimplexId greyVertexNumber = 0;
2033 // iterate over the edges and compute the edge intersection (range)
2034 for(SimplexId j = 0; j < (SimplexId)triangleEdges[i].size(); j++) {
2035
2036 SimplexId vertexId0 = 0, vertexId1 = 0;
2037 if(triangulation) {
2038 triangulation->getCellVertex(
2039 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j]], vertexId0);
2040 triangulation->getCellVertex(
2041 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1],
2042 vertexId1);
2043 } else {
2044 vertexId0
2045 = tetList_[5 * tetId + 1
2046 + edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2047 vertexId1
2048 = tetList_[5 * tetId + 1
2049 + edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2050 }
2051
2052 if((j < (SimplexId)triangleEdges[i].size() - 1)
2053 && (triangleEdges[i][j] == triangleEdges[i][j + 1])) {
2054
2055 // special case of a jacobi edge
2056 uv[j].first = ((const dataTypeU *)uField_)[vertexId0];
2057 uv[j].second = ((const dataTypeV *)vField_)[vertexId0];
2058
2059 uv[j + 1].first = ((const dataTypeU *)uField_)[vertexId1];
2060 uv[j + 1].second = ((const dataTypeV *)vField_)[vertexId1];
2061
2062 } else if((!j)
2063 || ((j)
2064 && (triangleEdges[i][j] != triangleEdges[i][j - 1]))) {
2065
2066 // regular intersection case
2067 d0 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2068 uv0.first = ((const dataTypeU *)uField_)[vertexId0];
2069 uv0.second = ((const dataTypeV *)vField_)[vertexId0];
2070
2071 d1 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2072 uv1.first = ((const dataTypeU *)uField_)[vertexId1];
2073 uv1.second = ((const dataTypeV *)vField_)[vertexId1];
2074
2075 uv[j].first
2076 = uv0.first + (d0 / (d0 - d1)) * (uv1.first - uv0.first);
2077 uv[j].second
2078 = uv0.second + (d0 / (d0 - d1)) * (uv1.second - uv0.second);
2079 }
2080
2081 // now determine the line parameterization of this intersection
2082 // point on the saddle edge
2083 if(fabs(rangePoint1.first - rangePoint0.first)
2084 > fabs(rangePoint1.second - rangePoint0.second)) {
2085 t[j] = (uv[j].first - rangePoint0.first)
2086 / (rangePoint1.first - rangePoint0.first);
2087 } else {
2088 t[j] = (uv[j].second - rangePoint0.second)
2089 / (rangePoint1.second - rangePoint0.second);
2090 }
2091
2092 if((t[j] <= 1) && (t[j] >= 0))
2093 greyVertexNumber++;
2094 else if(t[j] < 0)
2095 lowerVertexNumber++;
2096 else
2097 upperVertexNumber++;
2098 }
2099 // at this point, we know the uv coordinates (and the edge param)
2100 // of each vertex of the base triangle.
2101 // we can proceed with the cropping
2102
2103 // 4. triangulate the result
2104 if(greyVertexNumber == 3) {
2105 createdVertices += computeCase0<dataTypeU, dataTypeV>(
2106 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2107 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2108 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2109 triangulation);
2110 } else if(lowerVertexNumber == 3 || upperVertexNumber == 3) {
2111 // well do nothing (empty triangle)
2112 } else if((lowerVertexNumber == 1) && (upperVertexNumber == 1)
2113 && (greyVertexNumber == 1)) {
2114 createdVertices += computeCase1<dataTypeU, dataTypeV>(
2115 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2116 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2117 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2118 triangulation);
2119 } else if(((lowerVertexNumber == 2) && (upperVertexNumber == 1))
2120 || ((lowerVertexNumber == 1) && (upperVertexNumber == 2))) {
2121 createdVertices += computeCase2<dataTypeU, dataTypeV>(
2122 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2123 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2124 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2125 triangulation);
2126 } else if((greyVertexNumber == 1)
2127 && ((lowerVertexNumber == 2) || (upperVertexNumber == 2))) {
2128 createdVertices += computeCase3<dataTypeU, dataTypeV>(
2129 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2130 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2131 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2132 triangulation);
2133 } else if(((greyVertexNumber == 2))
2134 && ((lowerVertexNumber == 1) || (upperVertexNumber == 1))) {
2135 createdVertices += computeCase4<dataTypeU, dataTypeV>(
2136 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2137 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2138 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2139 triangulation);
2140 }
2141 }
2142 }
2143
2144 // in case of 2 triangles, remesh locally
2145 // NOTE: here we just snap vertices together if they are colinear
2146 // the mergeFiberSurfaces() function will take care of removing any
2147 // duplicate
2148 if((triangleEdges[1][0] != -1) && (createdVertices > 3)) {
2149
2150 std::vector<SimplexId> createdVertexList(createdVertices);
2151 for(SimplexId i = 0; i < (SimplexId)createdVertices; i++) {
2152 createdVertexList[i]
2153 = polygonEdgeVertexLists_[polygonEdgeId]->size() - 1 - i;
2154 }
2155
2156 std::vector<bool> snappedVertices(createdVertices, false);
2157 std::vector<SimplexId> colinearVertices;
2158 colinearVertices.reserve(createdVertices);
2159
2160 for(SimplexId i = 0; i < createdVertices; i++) {
2161 colinearVertices.clear();
2162
2163 if(!snappedVertices[i]) {
2164 colinearVertices.push_back(i);
2165 for(SimplexId j = 0; j < createdVertices; j++) {
2166 if((i != j) && (!snappedVertices[j])) {
2167 // not the same vertex
2168 // not snapped already
2169
2170 if((*polygonEdgeVertexLists_[polygonEdgeId])[createdVertexList[i]]
2171 .t_
2172 == (*polygonEdgeVertexLists_[polygonEdgeId])
2173 [createdVertexList[j]]
2174 .t_) {
2175 colinearVertices.push_back(j);
2176 }
2177 }
2178 }
2179 }
2180
2181 if((colinearVertices.size() == 4) || (colinearVertices.size() == 3)) {
2182 // 3 co-linear vertices with 1 duplicate
2183
2184 // we just need to find the pair of duplicates and snap both of
2185 // them to another vertex
2186 std::pair<SimplexId, SimplexId> minPair;
2187 double minDistance = -1;
2188 for(SimplexId j = 0; j < (SimplexId)colinearVertices.size(); j++) {
2189 for(SimplexId k = 0; k < (SimplexId)colinearVertices.size(); k++) {
2190 if(j != k) {
2191
2192 double distance = Geometry::distance(
2193 (*polygonEdgeVertexLists_[polygonEdgeId])
2194 [createdVertexList[colinearVertices[j]]]
2195 .p_.data(),
2196 (*polygonEdgeVertexLists_[polygonEdgeId])
2197 [createdVertexList[colinearVertices[k]]]
2198 .p_.data());
2199
2200 // bool basePointSnap = true;
2201 // for(int l = 0; l < 3; l++){
2202 // if((*polygonEdgeVertexLists_[polygonEdgeId])[
2203 // createdVertexList[colinearVertices[j]]].p_[l]
2204 // !=
2205 // (*polygonEdgeVertexLists_[polygonEdgeId])[
2206 // createdVertexList[colinearVertices[k]]].p_[l]){
2207 // basePointSnap = false;
2208 // break;
2209 // }
2210 // }
2211
2212 // if(!basePointSnap){
2213 if((minDistance < 0) || (distance < minDistance)) {
2214 minDistance = distance;
2215 minPair.first = j;
2216 minPair.second = k;
2217 }
2218 // }
2219 }
2220 }
2221 }
2222 if((minDistance != -1) && (minDistance < prec_dbl)) {
2223 // if((minDistance != -1)&&(minDistance <
2224 // pointSnappingThreshold_)){
2225 // snap them to another colinear vertex
2226 for(SimplexId j = 0; j < (SimplexId)colinearVertices.size(); j++) {
2227 if((j != minPair.first) && (j != minPair.second)) {
2228 // snap minPair.first to j
2229
2230 for(int k = 0; k < 3; k++) {
2231 (*polygonEdgeVertexLists_[polygonEdgeId])
2232 [createdVertexList[colinearVertices[minPair.first]]]
2233 .p_[k]
2234 = (*polygonEdgeVertexLists_[polygonEdgeId])
2235 [createdVertexList[colinearVertices[j]]]
2236 .p_[k];
2237 }
2238 (*polygonEdgeVertexLists_[polygonEdgeId])
2239 [createdVertexList[colinearVertices[minPair.first]]]
2240 .uv_
2241 = (*polygonEdgeVertexLists_[polygonEdgeId])
2242 [createdVertexList[colinearVertices[j]]]
2243 .uv_;
2244
2245 (*polygonEdgeVertexLists_[polygonEdgeId])
2246 [createdVertexList[colinearVertices[minPair.first]]]
2247 .t_
2248 = (*polygonEdgeVertexLists_[polygonEdgeId])
2249 [createdVertexList[colinearVertices[j]]]
2250 .t_;
2251
2252 (*polygonEdgeVertexLists_[polygonEdgeId])
2253 [createdVertexList[colinearVertices[minPair.first]]]
2254 .isBasePoint_
2255 = (*polygonEdgeVertexLists_[polygonEdgeId])
2256 [createdVertexList[colinearVertices[j]]]
2257 .isBasePoint_;
2258
2259 (*polygonEdgeVertexLists_[polygonEdgeId])
2260 [createdVertexList[colinearVertices[minPair.first]]]
2261 .isIntersectionPoint_
2262 = (*polygonEdgeVertexLists_[polygonEdgeId])
2263 [createdVertexList[colinearVertices[j]]]
2264 .isIntersectionPoint_;
2265 if((*polygonEdgeVertexLists_[polygonEdgeId])
2266 [createdVertexList[colinearVertices[j]]]
2267 .meshEdge_.first
2268 != -1) {
2269 (*polygonEdgeVertexLists_[polygonEdgeId])
2270 [createdVertexList[colinearVertices[minPair.first]]]
2271 .meshEdge_
2272 = (*polygonEdgeVertexLists_[polygonEdgeId])
2273 [createdVertexList[colinearVertices[j]]]
2274 .meshEdge_;
2275 }
2276
2277 // snap minPair.second to j
2278 for(int k = 0; k < 3; k++) {
2279 (*polygonEdgeVertexLists_[polygonEdgeId])
2280 [createdVertexList[colinearVertices[minPair.second]]]
2281 .p_[k]
2282 = (*polygonEdgeVertexLists_[polygonEdgeId])
2283 [createdVertexList[colinearVertices[j]]]
2284 .p_[k];
2285 }
2286 (*polygonEdgeVertexLists_[polygonEdgeId])
2287 [createdVertexList[colinearVertices[minPair.second]]]
2288 .uv_
2289 = (*polygonEdgeVertexLists_[polygonEdgeId])
2290 [createdVertexList[colinearVertices[j]]]
2291 .uv_;
2292
2293 (*polygonEdgeVertexLists_[polygonEdgeId])
2294 [createdVertexList[colinearVertices[minPair.second]]]
2295 .t_
2296 = (*polygonEdgeVertexLists_[polygonEdgeId])
2297 [createdVertexList[colinearVertices[j]]]
2298 .t_;
2299
2300 (*polygonEdgeVertexLists_[polygonEdgeId])
2301 [createdVertexList[colinearVertices[minPair.second]]]
2302 .isBasePoint_
2303 = (*polygonEdgeVertexLists_[polygonEdgeId])
2304 [createdVertexList[colinearVertices[j]]]
2305 .isBasePoint_;
2306
2307 (*polygonEdgeVertexLists_[polygonEdgeId])
2308 [createdVertexList[colinearVertices[minPair.second]]]
2309 .isIntersectionPoint_
2310 = (*polygonEdgeVertexLists_[polygonEdgeId])
2311 [createdVertexList[colinearVertices[j]]]
2312 .isIntersectionPoint_;
2313 if((*polygonEdgeVertexLists_[polygonEdgeId])
2314 [createdVertexList[colinearVertices[j]]]
2315 .meshEdge_.first
2316 != -1) {
2317 (*polygonEdgeVertexLists_[polygonEdgeId])
2318 [createdVertexList[colinearVertices[minPair.second]]]
2319 .meshEdge_
2320 = (*polygonEdgeVertexLists_[polygonEdgeId])
2321 [createdVertexList[colinearVertices[j]]]
2322 .meshEdge_;
2323 }
2324
2325 snappedVertices[colinearVertices[minPair.first]] = true;
2326 snappedVertices[colinearVertices[minPair.second]] = true;
2327
2328 // we're done
2329 break;
2330 }
2331 }
2332 }
2333 }
2334 }
2335 }
2336
2337 return createdVertices;
2338 }
2339
2340 return 0;
2341}
2342
2343template <class dataTypeU, class dataTypeV>
2345
2346#ifndef TTK_ENABLE_KAMIKAZE
2347 if(!tetNumber_)
2348 return -1;
2349#endif
2350
2351 // Algorithm
2352 // 1) loop over each triangle and mark the containing tetrahedron
2353 // 2) for each tet, in parallel, check for pairwise triangle intersections
2354 // (based on the triangles' range projection)
2355 // Given a pair of intersecting triangles:
2356 // 3) given the point of intersection, take one of its coordinates (u or v)
2357 // and compute an iso-contour I on both triangles
2358 // 4) express the barycentric coordinates of I in both triangles, 2 cases:
2359 // a) the intersection is the entire fiber (old code)
2360 // b) the intersection is a segment of fiber
2361
2362 // NOTE:
2363 // the topological aspect of the code is OK.
2364 // if any bug, it's very likely to be geometry accuracy related.
2365
2366 std::vector<std::vector<IntersectionTriangle>> tetIntersections(tetNumber_);
2367 std::vector<std::vector<Vertex>> tetNewVertices(tetNumber_);
2368
2369 // fill the information prior to the parallel pass
2370 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
2371 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
2372 j++) {
2373
2374 SimplexId tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
2375
2376 tetIntersections[tetId].emplace_back();
2377
2378 tetIntersections[tetId].back().caseId_
2379 = (*polygonEdgeTriangleLists_[i])[j].caseId_;
2380 tetIntersections[tetId].back().polygonEdgeId_ = i;
2381 tetIntersections[tetId].back().triangleId_ = j;
2382 for(int k = 0; k < 3; k++) {
2383 tetIntersections[tetId].back().vertexIds_[k]
2384 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
2385 tetIntersections[tetId].back().uv_[k]
2386 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2387 .uv_;
2388 tetIntersections[tetId].back().t_[k]
2389 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2390 .t_;
2391 for(int l = 0; l < 3; l++) {
2392 tetIntersections[tetId].back().p_[k][l]
2393 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2394 .p_[l];
2395 }
2396 }
2397 tetIntersections[tetId].back().intersection_.first = -DBL_MAX;
2398 tetIntersections[tetId].back().intersection_.second = -DBL_MAX;
2399 }
2400 }
2401
2402 std::vector<SimplexId> tetList;
2403 for(SimplexId i = 0; i < (SimplexId)tetIntersections.size(); i++) {
2404 if(tetIntersections[i].size() > 1)
2405 tetList.push_back(i);
2406 }
2407
2408#ifdef TTK_ENABLE_OPENMP
2409#pragma omp parallel for num_threads(threadNumber_)
2410#endif
2411 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
2412 SimplexId tetId = tetList[i];
2413
2414 // pre-process by merging nearby vertices...?
2415
2416 SimplexId newTriangleNumber = 1;
2417 SimplexId newVertexNumber = 1;
2418
2419 for(SimplexId j = 0; j < (SimplexId)tetIntersections[tetId].size(); j++) {
2420
2421 if(j > 1000) {
2422 this->printWrn("Preventing an infinite loop!");
2423 this->printWrn("More than 1000 re-meshed triangles in tet #"
2424 + std::to_string(tetId));
2425 this->printWrn("Extra-thin triangles keep on intersecting?!");
2426 break;
2427 }
2428
2429 // if we re-mesh, we add new triangles at the end of the list.
2430 // there's no need to check intersections with those.
2431 SimplexId originalTriangleNumber
2432 = (SimplexId)tetIntersections[tetId].size();
2433
2434 for(SimplexId k = 0; k < originalTriangleNumber; k++) {
2435
2436 SimplexId polygonEdgeId0 = tetIntersections[tetId][j].polygonEdgeId_;
2437 SimplexId polygonEdgeId1 = tetIntersections[tetId][k].polygonEdgeId_;
2438
2439 if((j != k) && (polygonEdgeId0 != polygonEdgeId1)) {
2440 // cases 3, 4 and 6 of the fiber surface table (multiple triangle
2441 // per tet given a single edge). we don't need to re-mesh that.
2442
2443 // grab the range projection for the triangle j
2444 std::pair<double, double> edge0point0, edge0point1;
2445
2446 getTriangleRangeExtremities(
2447 tetId, j, tetIntersections, edge0point0, edge0point1);
2448
2449 // now do the same thing for the triangle k
2450 std::pair<double, double> edge1point0, edge1point1;
2451
2452 getTriangleRangeExtremities(
2453 tetId, k, tetIntersections, edge1point0, edge1point1);
2454
2455 // compute the intersection
2456 std::pair<double, double> intersection;
2457 bool hasIntersection = Geometry::computeSegmentIntersection(
2458 edge0point0.first, edge0point0.second, edge0point1.first,
2459 edge0point1.second, edge1point0.first, edge1point0.second,
2460 edge1point1.first, edge1point1.second, intersection.first,
2461 intersection.second);
2462
2463 if((hasIntersection)
2464 // check if that intersection has been registered before
2465 // in the end, only one intersection per triangle, no matter what
2466 /*&&(((fabs(tetIntersections[tetId][j].intersection_.first
2467 - intersection.first) > Geometry::powIntTen(-FLT_DIG))
2468 ||(fabs(tetIntersections[tetId][j].intersection_.second
2469 - intersection.second) > Geometry::powIntTen(-FLT_DIG)))
2470 &&((fabs(tetIntersections[tetId][k].intersection_.first
2471 - intersection.first) > Geometry::powIntTen(-FLT_DIG))
2472 ||(fabs(tetIntersections[tetId][k].intersection_.second
2473 - intersection.second) > Geometry::powIntTen(-FLT_DIG))))*/){
2474
2475 computeTriangleIntersection(tetId, j, k, polygonEdgeId0,
2476 polygonEdgeId1, intersection,
2477 newVertexNumber, newTriangleNumber,
2478 tetIntersections, tetNewVertices);
2479 }
2480 }
2481 }
2482 }
2483 }
2484
2485 // now copy the new vertices
2486 for(SimplexId i = 0; i < (SimplexId)tetNewVertices.size(); i++) {
2487 for(SimplexId j = 0; j < (SimplexId)tetNewVertices[i].size(); j++) {
2488
2489 SimplexId localId = (*globalVertexList_).size();
2490 tetNewVertices[i][j].localId_ = localId;
2491 (*globalVertexList_).push_back(tetNewVertices[i][j]);
2492 }
2493 }
2494
2495 for(SimplexId i = 0; i < (SimplexId)tetIntersections.size(); i++) {
2496 for(SimplexId j = 0; j < (SimplexId)tetIntersections[i].size(); j++) {
2497
2498 if(((tetIntersections[i][j].intersection_.first != -DBL_MAX)
2499 && (tetIntersections[i][j].intersection_.second != -DBL_MAX))
2500 || (tetIntersections[i][j].triangleId_ < 0)) {
2501
2502 SimplexId triangleId = tetIntersections[i][j].triangleId_;
2503
2504 if(triangleId < 0) {
2505
2506 // this is a new triangle
2507 triangleId = (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2508 .polygonEdgeId_])
2509 .size();
2510 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2511 .resize(triangleId + 1);
2512 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2513 .back()
2514 .tetId_
2515 = i;
2516 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2517 .back()
2518 .caseId_
2519 = tetIntersections[i][j].caseId_;
2520 }
2521
2522 for(int k = 0; k < 3; k++) {
2523
2524 SimplexId vertexId = tetIntersections[i][j].vertexIds_[k];
2525
2526 if(vertexId < 0) {
2527 // newly created vertex
2528 vertexId = tetNewVertices[i][-(vertexId + 1)].localId_;
2529 }
2530 (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2531 .polygonEdgeId_])[triangleId]
2532 .vertexIds_[k]
2533 = vertexId;
2534 }
2535 }
2536 }
2537 }
2538
2539 return 0;
2540}
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
int threadNumber_
Definition: BaseClass.h:95
Minimalist debugging class.
Definition: Debug.h:88
int debugLevel_
Definition: Debug.h:379
virtual int setDebugLevel(const int &debugLevel)
Definition: Debug.cpp:147
TTK processing package that computes fiber surfaces.
Definition: FiberSurface.h:49
std::vector< std::vector< Vertex > * > polygonEdgeVertexLists_
Definition: FiberSurface.h:497
const void * vField_
Definition: FiberSurface.h:483
int setGlobalVertexList(std::vector< Vertex > *globalList)
Definition: FiberSurface.h:126
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
std::array< SimplexId, 12 > edgeImplicitEncoding_
Definition: FiberSurface.h:487
double pointSnappingThreshold_
Definition: FiberSurface.h:491
double edgeCollapseThreshold_
Definition: FiberSurface.h:490
std::vector< std::vector< Triangle > * > polygonEdgeTriangleLists_
Definition: FiberSurface.h:498
const std::vector< std::vector< SimplexId > > * tetNeighbors_
Definition: FiberSurface.h:486
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 setInputField(const void *uField, const void *vField)
Definition: FiberSurface.h:131
int setTetList(const SimplexId *tetList)
Definition: FiberSurface.h:173
int computeCase3(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) 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 finalize(const bool &mergeDuplicatedVertices=false, const bool &removeSmallEdges=false, const bool &edgeFlips=false, const bool &intersectionRemesh=false)
const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > * polygon_
Definition: FiberSurface.h:494
int flipEdges() const
int computeCase4(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
int computeCase1(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
Definition: FiberSurface.h:801
int computeSurface(const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
int remeshIntersections() const
int computeCase0(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
Definition: FiberSurface.h:656
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
int processTetrahedron(const SimplexId &tetId, const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
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 setPolygonEdgeNumber(const SimplexId &polygonEdgeNumber)
Definition: FiberSurface.h:166
int mergeVertices(const double &distanceThreshold) const
int setTetNumber(const SimplexId &tetNumber)
Definition: FiberSurface.h:184
int computeCase2(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
Definition: FiberSurface.h:990
std::vector< Vertex > * globalVertexList_
Definition: FiberSurface.h:496
SimplexId pointNumber_
Definition: FiberSurface.h:482
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
int setPointSet(const float *pointSet)
Definition: FiberSurface.h:154
int setPointNumber(const SimplexId &number)
Definition: FiberSurface.h:149
int setPointMerging(const bool &onOff)
Definition: FiberSurface.h:139
int setPointMergingThreshold(const double &threshold)
Definition: FiberSurface.h:144
int computeBaseTriangle(const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, std::array< std::array< double, 3 >, 3 > &basePoints, std::array< std::pair< double, double >, 3 > &basePointProjections, std::array< double, 3 > &basePointParameterization, std::array< std::pair< SimplexId, SimplexId >, 3 > &baseEdges, const triangulationType *const triangulation) const
Definition: FiberSurface.h:536
int setPolygon(const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > *polygon)
Definition: FiberSurface.h:159
const SimplexId * tetList_
Definition: FiberSurface.h:485
const void * uField_
Definition: FiberSurface.h:483
SimplexId polygonEdgeNumber_
Definition: FiberSurface.h:482
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
int setVertexList(const SimplexId &polygonEdgeId, std::vector< Vertex > *vertexList)
Definition: FiberSurface.h:209
int setTetNeighbors(const std::vector< std::vector< SimplexId > > *tetNeighbors)
Definition: FiberSurface.h:179
int setTriangleList(const SimplexId &polygonEdgeId, std::vector< Triangle > *triangleList)
Definition: FiberSurface.h:189
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
void preconditionTriangulation(AbstractTriangulation *triangulation)
Definition: FiberSurface.h:202
int computeContour(const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const std::vector< SimplexId > &seedTetList, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
TTK optional package for octree based range queries in bivariate volumetric data.
double getElapsedTime()
Definition: Timer.h:15
bool computeSegmentIntersection(const T &xA, const T &yA, const T &xB, const T &yB, const T &xC, const T &yC, const T &xD, const T &yD, T &x, T &y)
Definition: Geometry.cpp:230
int computeBarycentricCoordinates(const T *p0, const T *p1, const T *p, std::array< T, 2 > &baryCentrics, const int &dimension=3)
Definition: Geometry.cpp:124
bool isTriangleColinear(const T *p0, const T *p1, const T *p2, const T *tolerance=nullptr)
Definition: Geometry.cpp:448
T powIntTen(const int n)
Compute the nth power of ten.
Definition: Geometry.h:360
T powInt(const T val, const int n)
Definition: Geometry.h:338
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< std::array< double, 3 >, 3 > p_
Definition: FiberSurface.h:231
std::pair< double, double > intersection_
Definition: FiberSurface.h:232
std::array< std::pair< double, double >, 3 > uv_
Definition: FiberSurface.h:229
std::array< SimplexId, 3 > vertexIds_
Definition: FiberSurface.h:228
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
std::pair< SimplexId, SimplexId > meshEdge_
Definition: FiberSurface.h:57
printMsg(debug::output::GREEN+"                           "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)