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 {
53 bool isBasePoint_{}, isIntersectionPoint_{};
54 SimplexId localId_{}, globalId_{}, polygonEdgeId_{};
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_{};
65 SimplexId tetId_{}, caseId_{}, polygonEdgeId_{};
66 };
67
68 FiberSurface();
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
202 preconditionTriangulation(AbstractTriangulation *triangulation) {
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:
222 struct IntersectionTriangle {
223 SimplexId caseId_{};
224 // use negative values for new triangles
225 SimplexId triangleId_{};
226 SimplexId polygonEdgeId_{};
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
341 int computeTriangleFiber(
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
351 int computeTriangleIntersection(
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
363 int computeTriangleIntersection(
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
376 int createNewIntersectionTriangle(
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
392 int getNumberOfCommonVertices(
393 const SimplexId &tetId,
394 const SimplexId &triangleId0,
395 const SimplexId &triangleId1,
396 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections)
397 const;
398
399 int getTriangleRangeExtremities(
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
419 bool isEdgeAngleCollapsible(
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
430 inline bool isIntersectionTriangleColinear(
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
458 return Geometry::isTriangleColinear(
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
476 int snapVertexBarycentrics(
477 const SimplexId &tetId,
478 const std::vector<std::pair<SimplexId, SimplexId>> &triangles) const;
479
480 bool pointSnapping_{false};
481
482 SimplexId pointNumber_{}, tetNumber_{}, polygonEdgeNumber_{};
483 const void *uField_{}, *vField_{};
484 const float *pointSet_{};
485 const SimplexId *tetList_{};
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
490 double edgeCollapseThreshold_{Geometry::powIntTen(-FLT_DIG + 2)},
491 pointSnappingThreshold_{Geometry::powIntTen(-FLT_DIG + 1)};
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>
536inline int ttk::FiberSurface::computeBaseTriangle(
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;
622 Geometry::computeBarycentricCoordinates(
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>
656inline int ttk::FiberSurface::computeCase0(
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 const 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;
764 Geometry::computeBarycentricCoordinates(
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>
801inline int ttk::FiberSurface::computeCase1(
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 const 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 const triangleId
832 = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
833 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 3);
834
835 for(int i = 0; i < 3; i++) {
836
837 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
838 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 1;
839 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
840 = polygonEdgeId;
841
842 switch(i) {
843 case 0:
844 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
845 .vertexIds_[0]
846 = vertexId;
847 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
848 .vertexIds_[1]
849 = vertexId + 1;
850 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
851 .vertexIds_[2]
852 = vertexId + 2;
853 break;
854 case 1:
855 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
856 .vertexIds_[0]
857 = vertexId + 1;
858 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
859 .vertexIds_[1]
860 = vertexId + 2;
861 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
862 .vertexIds_[2]
863 = vertexId + 3;
864 break;
865 case 2:
866 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
867 .vertexIds_[0]
868 = vertexId + 2;
869 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
870 .vertexIds_[1]
871 = vertexId + 3;
872 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
873 .vertexIds_[2]
874 = vertexId + 4;
875 break;
876 }
877 }
878
879 // compute the base triangle vertices like in case 1
880 std::array<std::array<double, 3>, 3> basePoints{};
881 std::array<std::pair<double, double>, 3> basePointProjections{};
882 std::array<double, 3> basePointParameterization{};
883 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
884
885 computeBaseTriangle<dataTypeU, dataTypeV>(
886 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
887 u2, v2, basePoints, basePointProjections, basePointParameterization,
888 baseEdges, triangulation);
889
890 // find the pivot vertex for this case
891 SimplexId pivotVertexId = -1;
892
893 if((t0 >= 0) && (t0 <= 1)) {
894 pivotVertexId = 0;
895 }
896 if((t1 >= 0) && (t1 <= 1)) {
897 pivotVertexId = 1;
898 }
899 if((t2 >= 0) && (t2 <= 1)) {
900 pivotVertexId = 2;
901 }
902
903 // now get the vertex coordinates
904 for(int i = 0; i < 5; i++) {
905
906 SimplexId vertexId0 = -1, vertexId1 = -1;
907 double t{};
908
909 if(!i) {
910 // just take the pivot vertex
911 for(int j = 0; j < 3; j++) {
912 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
913 = basePoints[pivotVertexId][j];
914 }
915
916 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
917 = basePointParameterization[pivotVertexId];
918 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
919 = basePointProjections[pivotVertexId];
920 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
921 = baseEdges[pivotVertexId];
922 } else {
923
924 switch(i) {
925
926 case 1:
927 // interpolation between pivotVertexId and (pivotVertexId-1)%3
928 // if pivot is positive: interpolate at 1
929 // if not interpolate at 0
930 vertexId0 = pivotVertexId;
931 vertexId1 = (pivotVertexId + 2) % 3;
932 if(basePointParameterization[(pivotVertexId + 2) % 3] > 1) {
933 t = 1;
934 } else {
935 t = 0;
936 }
937 break;
938
939 case 2:
940 // interpolation between pivotVertexId and (pivotVertexId+1)%3
941 // if pivot is positive: interpolate at 1
942 // if not interpolate at 0
943 vertexId0 = pivotVertexId;
944 vertexId1 = (pivotVertexId + 1) % 3;
945 if(basePointParameterization[(pivotVertexId + 1) % 3] > 1) {
946 t = 1;
947 } else {
948 t = 0;
949 }
950 break;
951
952 case 3:
953 vertexId0 = (pivotVertexId + 2) % 3;
954 vertexId1 = (pivotVertexId + 1) % 3;
955 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
956 t = 0;
957 } else {
958 t = 1;
959 }
960 break;
961
962 case 4:
963 vertexId0 = (pivotVertexId + 2) % 3;
964 vertexId1 = (pivotVertexId + 1) % 3;
965 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
966 t = 1;
967 } else {
968 t = 0;
969 }
970 break;
971 }
972
973 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
974
975 interpolateBasePoints(
976 basePoints[vertexId0], basePointProjections[vertexId0],
977 basePointParameterization[vertexId0], basePoints[vertexId1],
978 basePointProjections[vertexId1], basePointParameterization[vertexId1],
979 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
980 // snapToBasePoint(
981 // basePoints, basePointProjections, basePointParameterization,
982 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
983 }
984 }
985
986 // return the number of created vertices
987 return 5;
988}
989
990template <class dataTypeU, class dataTypeV, typename triangulationType>
991inline int ttk::FiberSurface::computeCase2(
992 const SimplexId &polygonEdgeId,
993 const SimplexId &tetId,
994 const SimplexId &localEdgeId0,
995 const double &t0,
996 const double &u0,
997 const double &v0,
998 const SimplexId &localEdgeId1,
999 const double &t1,
1000 const double &u1,
1001 const double &v1,
1002 const SimplexId &localEdgeId2,
1003 const double &t2,
1004 const double &u2,
1005 const double &v2,
1006 const triangulationType *const triangulation) const {
1007
1008 SimplexId const vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1009
1010 // alloc 4 more vertices
1011 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1012 for(int i = 0; i < 4; i++) {
1013 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1014 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1015 = false;
1016 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1017 = std::pair<SimplexId, SimplexId>(-1, -1);
1018 }
1019
1020 // alloc 2 more triangles
1021 SimplexId const triangleId
1022 = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1023 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1024
1025 for(int i = 0; i < 2; i++) {
1026
1027 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1028 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 2;
1029 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1030 = polygonEdgeId;
1031
1032 if(!i) {
1033 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1034 = vertexId;
1035 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1036 = vertexId + 1;
1037 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1038 = vertexId + 2;
1039 } else {
1040 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1041 = vertexId + 1;
1042 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1043 = vertexId + 3;
1044 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1045 = vertexId + 2;
1046 }
1047 }
1048
1049 // compute the base triangle vertices like in case 1
1050 std::array<std::array<double, 3>, 3> basePoints{};
1051 std::array<std::pair<double, double>, 3> basePointProjections{};
1052 std::array<double, 3> basePointParameterization{};
1053 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1054
1055 computeBaseTriangle<dataTypeU, dataTypeV>(
1056 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1057 u2, v2, basePoints, basePointProjections, basePointParameterization,
1058 baseEdges, triangulation);
1059
1060 // find the pivot for this case
1061 bool isPivotPositive = false;
1062 SimplexId pivotVertexId = -1;
1063
1064 if(((t0 < 0) && ((t1 < 0) || (t2 < 0)))
1065 || ((t1 < 0) && ((t0 < 0) || (t2 < 0)))
1066 || ((t2 < 0) && ((t1 < 0) || (t0 < 0)))) {
1067 isPivotPositive = true;
1068 }
1069 if(isPivotPositive) {
1070 if(t0 >= 1)
1071 pivotVertexId = 0;
1072 } else {
1073 if(t0 <= 0)
1074 pivotVertexId = 0;
1075 }
1076 if(isPivotPositive) {
1077 if(t1 >= 1)
1078 pivotVertexId = 1;
1079 } else {
1080 if(t1 <= 0)
1081 pivotVertexId = 1;
1082 }
1083 if(isPivotPositive) {
1084 if(t2 >= 1)
1085 pivotVertexId = 2;
1086 } else {
1087 if(t2 <= 0)
1088 pivotVertexId = 2;
1089 }
1090
1091 // now get the vertex coordinates
1092 for(int i = 0; i < 4; i++) {
1093
1094 SimplexId vertexId0 = -1, vertexId1 = -1;
1095 double t{};
1096
1097 switch(i) {
1098
1099 case 0:
1100 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1101 // if pivot is positive: interpolate at 1
1102 // if not interpolate at 0
1103 vertexId0 = pivotVertexId;
1104 vertexId1 = (pivotVertexId + 2) % 3;
1105 if(isPivotPositive)
1106 t = 1;
1107 else
1108 t = 0;
1109 break;
1110
1111 case 1:
1112 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1113 // if pivot is positive: interpolate at 1
1114 // if not interpolate at 0
1115 vertexId0 = pivotVertexId;
1116 vertexId1 = (pivotVertexId + 1) % 3;
1117 if(isPivotPositive)
1118 t = 1;
1119 else
1120 t = 0;
1121 break;
1122
1123 case 2:
1124 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1125 // if pivot is positive: interpolate at 0
1126 // if not interpolate at 1
1127 vertexId0 = pivotVertexId;
1128 vertexId1 = (pivotVertexId + 2) % 3;
1129 if(isPivotPositive)
1130 t = 0;
1131 else
1132 t = 1;
1133 break;
1134
1135 case 3:
1136 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1137 // if pivot is positive: interpolate at 0
1138 // if not interpolate at 1
1139 vertexId0 = pivotVertexId;
1140 vertexId1 = (pivotVertexId + 1) % 3;
1141 if(isPivotPositive)
1142 t = 0;
1143 else
1144 t = 1;
1145 break;
1146 }
1147
1148 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1149
1150 interpolateBasePoints(
1151 basePoints[vertexId0], basePointProjections[vertexId0],
1152 basePointParameterization[vertexId0], basePoints[vertexId1],
1153 basePointProjections[vertexId1], basePointParameterization[vertexId1], t,
1154 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1155 // snapToBasePoint(
1156 // basePoints, basePointProjections, basePointParameterization,
1157 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1158 }
1159
1160 // return the number of created vertices
1161 return 4;
1162}
1163
1164template <class dataTypeU, class dataTypeV, typename triangulationType>
1165inline int ttk::FiberSurface::computeCase3(
1166 const SimplexId &polygonEdgeId,
1167 const SimplexId &tetId,
1168 const SimplexId &localEdgeId0,
1169 const double &t0,
1170 const double &u0,
1171 const double &v0,
1172 const SimplexId &localEdgeId1,
1173 const double &t1,
1174 const double &u1,
1175 const double &v1,
1176 const SimplexId &localEdgeId2,
1177 const double &t2,
1178 const double &u2,
1179 const double &v2,
1180 const triangulationType *const triangulation) const {
1181
1182 SimplexId const vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1183
1184 // alloc 3 more vertices
1185 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 3);
1186 for(int i = 0; i < 3; i++) {
1187 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1188 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1189 = false;
1190 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1191 = std::pair<SimplexId, SimplexId>(-1, -1);
1192 }
1193
1194 // alloc 1 more triangle
1195 SimplexId const triangleId
1196 = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1197 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 1);
1198
1199 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].tetId_ = tetId;
1200 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].caseId_ = 3;
1201 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().polygonEdgeId_
1202 = polygonEdgeId;
1203
1204 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[0]
1205 = vertexId;
1206 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[1]
1207 = vertexId + 1;
1208 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[2]
1209 = vertexId + 2;
1210
1211 // compute the base triangle vertices like in case 1
1212 std::array<std::array<double, 3>, 3> basePoints{};
1213 std::array<std::pair<double, double>, 3> basePointProjections{};
1214 std::array<double, 3> basePointParameterization{};
1215 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1216
1217 computeBaseTriangle<dataTypeU, dataTypeV>(
1218 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1219 u2, v2, basePoints, basePointProjections, basePointParameterization,
1220 baseEdges, triangulation);
1221
1222 // now find the pivot
1223 bool isPivotPositive = false;
1224 SimplexId pivotVertexId = -1;
1225
1226 if((t0 <= 1) && (t0 >= 0)) {
1227 pivotVertexId = 0;
1228 } else {
1229 if(t0 < 0)
1230 isPivotPositive = true;
1231 }
1232 if((t1 <= 1) && (t1 >= 0)) {
1233 pivotVertexId = 1;
1234 } else {
1235 if(t1 < 0)
1236 isPivotPositive = true;
1237 }
1238 if((t2 <= 1) && (t2 >= 0)) {
1239 pivotVertexId = 2;
1240 } else {
1241 if(t2 < 0)
1242 isPivotPositive = true;
1243 }
1244
1245 // now get the vertex coordinates
1246 for(int i = 0; i < 3; i++) {
1247
1248 SimplexId vertexId0 = -1, vertexId1 = -1;
1249 double t{};
1250
1251 if(!i) {
1252 // special case of the pivot vertex
1253 for(int j = 0; j < 3; j++) {
1254 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
1255 = basePoints[pivotVertexId][j];
1256 }
1257
1258 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
1259 = basePointParameterization[pivotVertexId];
1260 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
1261 = basePointProjections[pivotVertexId];
1262 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
1263 = baseEdges[pivotVertexId];
1264 } else {
1265 if(i == 1) {
1266 // interpolation between pivotVertexId and pivotVertexId+1
1267 // if pivot is positive, interpolate at value 0
1268 // else interpolate at value 1
1269 vertexId0 = pivotVertexId;
1270 vertexId1 = (pivotVertexId + 1) % 3;
1271 if(isPivotPositive)
1272 t = 0;
1273 else
1274 t = 1;
1275 }
1276 if(i == 2) {
1277 // interpolation between pivotVertexId and pivotVertexId-1
1278 // if pivot is positive, interpolate at value 0
1279 // else interpolate at value 1
1280 vertexId0 = pivotVertexId;
1281 vertexId1 = (pivotVertexId + 2) % 3;
1282 if(isPivotPositive)
1283 t = 0;
1284 else
1285 t = 1;
1286 }
1287
1288 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1289
1290 interpolateBasePoints(
1291 basePoints[vertexId0], basePointProjections[vertexId0],
1292 basePointParameterization[vertexId0], basePoints[vertexId1],
1293 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1294 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1295 // snapToBasePoint(
1296 // basePoints, basePointProjections, basePointParameterization,
1297 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1298 }
1299 }
1300
1301 // return the number of created vertices
1302 return 3;
1303}
1304
1305template <class dataTypeU, class dataTypeV, typename triangulationType>
1306inline int ttk::FiberSurface::computeCase4(
1307 const SimplexId &polygonEdgeId,
1308 const SimplexId &tetId,
1309 const SimplexId &localEdgeId0,
1310 const double &t0,
1311 const double &u0,
1312 const double &v0,
1313 const SimplexId &localEdgeId1,
1314 const double &t1,
1315 const double &u1,
1316 const double &v1,
1317 const SimplexId &localEdgeId2,
1318 const double &t2,
1319 const double &u2,
1320 const double &v2,
1321 const triangulationType *const triangulation) const {
1322
1323 SimplexId const vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1324
1325 // alloc 4 more vertices
1326 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1327 for(int i = 0; i < 4; i++) {
1328 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ = true;
1329 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1330 = false;
1331 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1332 = std::pair<SimplexId, SimplexId>(-1, -1);
1333 }
1334
1335 // alloc 2 more triangles
1336 SimplexId const triangleId
1337 = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1338 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1339
1340 for(int i = 0; i < 2; i++) {
1341
1342 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1343 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 4;
1344 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1345 = polygonEdgeId;
1346
1347 if(!i) {
1348 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1349 = vertexId;
1350 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1351 = vertexId + 1;
1352 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1353 = vertexId + 2;
1354 } else {
1355 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1356 = vertexId + 1;
1357 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1358 = vertexId + 3;
1359 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1360 = vertexId + 2;
1361 }
1362 }
1363
1364 // compute the base triangle vertices like in case 1
1365 std::array<std::array<double, 3>, 3> basePoints{};
1366 std::array<std::pair<double, double>, 3> basePointProjections{};
1367 std::array<double, 3> basePointParameterization{};
1368 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1369
1370 computeBaseTriangle<dataTypeU, dataTypeV>(
1371 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1372 u2, v2, basePoints, basePointProjections, basePointParameterization,
1373 baseEdges, triangulation);
1374
1375 // find the pivot vertex for this case
1376 bool isPivotPositive = false;
1377 SimplexId pivotVertexId = -1;
1378
1379 if(t0 > 1) {
1380 pivotVertexId = 0;
1381 isPivotPositive = true;
1382 } else if(t0 < 0) {
1383 pivotVertexId = 0;
1384 isPivotPositive = false;
1385 }
1386
1387 if(t1 > 1) {
1388 pivotVertexId = 1;
1389 isPivotPositive = true;
1390 } else if(t1 < 0) {
1391 pivotVertexId = 1;
1392 isPivotPositive = false;
1393 }
1394 if(t2 > 1) {
1395 pivotVertexId = 2;
1396 isPivotPositive = true;
1397 } else if(t2 < 0) {
1398 pivotVertexId = 2;
1399 isPivotPositive = false;
1400 }
1401
1402 // now get the vertex coordinates depending on the case
1403 for(int i = 0; i < 4; i++) {
1404
1405 SimplexId vertexId0 = -1, vertexId1 = -1;
1406 double t{};
1407
1408 if(i < 2) {
1409 if(!i) {
1410 // interpolation between pivotVertexId and (pivotVertexId-1)%3
1411 // if pivot is positive: interpolate at 1
1412 // if not interpolate at 0
1413 vertexId0 = pivotVertexId;
1414 vertexId1 = (pivotVertexId + 2) % 3;
1415 if(isPivotPositive)
1416 t = 1;
1417 else
1418 t = 0;
1419 } else {
1420 // interpolation between pivotVertexId and (pivotVertexId+1)%3
1421 // if pivot is positive: interpolate at 1
1422 // if not interpolate at 0
1423 vertexId0 = pivotVertexId;
1424 vertexId1 = (pivotVertexId + 1) % 3;
1425 if(isPivotPositive)
1426 t = 1;
1427 else
1428 t = 0;
1429 }
1430
1431 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1432
1433 interpolateBasePoints(
1434 basePoints[vertexId0], basePointProjections[vertexId0],
1435 basePointParameterization[vertexId0], basePoints[vertexId1],
1436 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1437 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1438 // snapToBasePoint(
1439 // basePoints, basePointProjections, basePointParameterization,
1440 // (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1441
1442 } else {
1443 if(i == 2) {
1444 // take (pivotVertexId-1)%3
1445 for(int j = 0; j < 3; j++) {
1446 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1447 = basePoints[(pivotVertexId + 2) % 3][j];
1448 }
1449
1450 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1451 = basePointParameterization[(pivotVertexId + 2) % 3];
1452 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1453 = basePointProjections[(pivotVertexId + 2) % 3];
1454 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1455 = baseEdges[(pivotVertexId + 2) % 3];
1456 } else {
1457 // take (pivtoVertexId+1)%3
1458 for(int j = 0; j < 3; j++) {
1459 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1460 = basePoints[(pivotVertexId + 1) % 3][j];
1461 }
1462 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1463 = basePointParameterization[(pivotVertexId + 1) % 3];
1464 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1465 = basePointProjections[(pivotVertexId + 1) % 3];
1466 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1467 = baseEdges[(pivotVertexId + 1) % 3];
1468 }
1469 }
1470 }
1471
1472 // return the number of created vertices
1473 return 4;
1474}
1475
1476template <class dataTypeU, class dataTypeV, typename triangulationType>
1477inline int ttk::FiberSurface::computeContour(
1478 const std::pair<double, double> &rangePoint0,
1479 const std::pair<double, double> &rangePoint1,
1480 const std::vector<SimplexId> &seedTetList,
1481 const triangulationType *const triangulation,
1482 const SimplexId &polygonEdgeId) const {
1483
1484#ifndef TTK_ENABLE_KAMIKAZE
1485 if(!uField_)
1486 return -4;
1487 if(!vField_)
1488 return -5;
1489 if(!triangulation)
1490 return -6;
1491 if(!polygonEdgeNumber_)
1492 return -7;
1493 if(!globalVertexList_)
1494 return -8;
1495#endif
1496
1497 std::vector<bool> visitedTets(triangulation->getNumberOfCells(), false);
1498 std::queue<SimplexId> tetQueue;
1499
1500 // init the queue
1501 for(SimplexId i = 0; i < (SimplexId)seedTetList.size(); i++) {
1502 tetQueue.push(seedTetList[i]);
1503 }
1504
1505 SimplexId createdVertices = 0;
1506
1507 do {
1508
1509 SimplexId const tetId = tetQueue.front();
1510 tetQueue.pop();
1511
1512 if(!visitedTets[tetId]) {
1513
1514 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1515 tetId, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1516
1517 if(createdVertices) {
1518 // only propagate if we created a triangle
1519 SimplexId const tetNeighborNumber
1520 = triangulation->getCellNeighborNumber(tetId);
1521
1522 for(SimplexId i = 0; i < tetNeighborNumber; i++) {
1523
1524 SimplexId neighborId = -1;
1525 triangulation->getCellNeighbor(tetId, i, neighborId);
1526
1527 if(!visitedTets[neighborId])
1528 tetQueue.push(neighborId);
1529 }
1530 }
1531
1532 visitedTets[tetId] = true;
1533 }
1534
1535 } while(tetQueue.size());
1536
1537 return 0;
1538}
1539
1540template <class dataTypeU, class dataTypeV>
1541int ttk::FiberSurface::computeContour(
1542 const std::vector<
1543 std::pair<std::pair<double, double>, std::pair<double, double>>> &edgeList,
1544 const std::vector<SimplexId> &seedTetList,
1545 const std::vector<SimplexId> *edgeIdList) const {
1546
1547#ifndef TTK_ENABLE_KAMIKAZE
1548 if(!tetNeighbors_)
1549 return -1;
1550 if(!tetNumber_)
1551 return -2;
1552 if(!tetList_)
1553 return -3;
1554 if(!uField_)
1555 return -4;
1556 if(!vField_)
1557 return -5;
1558 if(!pointSet_)
1559 return -6;
1560 if(!polygonEdgeNumber_)
1561 return -7;
1562 if(!globalVertexList_)
1563 return -8;
1564#endif
1565
1566 std::vector<bool> visitedTets(tetNumber_, false);
1567 std::queue<SimplexId> tetQueue;
1568
1569 // init the queue
1570 for(SimplexId i = 0; i < (SimplexId)seedTetList.size(); i++) {
1571 tetQueue.push(seedTetList[i]);
1572 }
1573
1574 SimplexId createdVertices = 0;
1575
1576 do {
1577
1578 SimplexId tetId = tetQueue.front();
1579 tetQueue.pop();
1580
1581 if(!visitedTets[tetId]) {
1582
1583 std::vector<std::vector<SimplexId>> threadedTetQueue(edgeList.size());
1584#ifdef TTK_ENABLE_OPENMP
1585#pragma omp parallel for num_threads(threadNumber_)
1586#endif
1587 for(SimplexId i = 0; i < (SimplexId)edgeList.size(); i++) {
1588
1589 SimplexId polygonEdgeId = 0;
1590
1591 if(edgeIdList) {
1592 polygonEdgeId = (*edgeIdList)[i];
1593 }
1594
1595 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1596 tetId, edgeList[i].first, edgeList[i].second, polygonEdgeId);
1597
1598 if(createdVertices) {
1599 // only propagate if we created a triangle
1600 for(SimplexId j = 0; j < (SimplexId)(*tetNeighbors_)[tetId].size();
1601 j++) {
1602 if(!visitedTets[(*tetNeighbors_)[tetId][j]]) {
1603 threadedTetQueue[i].push_back((*tetNeighbors_)[tetId][j]);
1604 }
1605 }
1606 }
1607 }
1608
1609 visitedTets[tetId] = true;
1610
1611 for(SimplexId i = 0; i < (SimplexId)threadedTetQueue.size(); i++) {
1612 for(SimplexId j = 0; j < (SimplexId)threadedTetQueue[i].size(); j++) {
1613 tetQueue.push(threadedTetQueue[i][j]);
1614 }
1615 }
1616 }
1617
1618 } while(tetQueue.size());
1619
1620 return 0;
1621}
1622
1623template <class dataTypeU, class dataTypeV, typename triangulationType>
1624inline int ttk::FiberSurface::computeSurface(
1625 const std::pair<double, double> &rangePoint0,
1626 const std::pair<double, double> &rangePoint1,
1627 const triangulationType *const triangulation,
1628 const SimplexId &polygonEdgeId) const {
1629
1630#ifndef TTK_ENABLE_KAMIKAZE
1631 if((!tetNumber_) && (!triangulation))
1632 return -1;
1633 if((!tetList_) && (!triangulation))
1634 return -2;
1635 if(!uField_)
1636 return -3;
1637 if(!vField_)
1638 return -4;
1639 if((!pointSet_) && (!triangulation))
1640 return -5;
1641 if(!polygonEdgeNumber_)
1642 return -6;
1643 if(!globalVertexList_)
1644 return -7;
1645#endif
1646
1647 SimplexId tetNumber = tetNumber_;
1648
1649 if(triangulation) {
1650 tetNumber = triangulation->getNumberOfCells();
1651 }
1652
1653#ifdef TTK_ENABLE_OPENMP
1654#pragma omp parallel for num_threads(threadNumber_)
1655#endif
1656 for(SimplexId i = 0; i < tetNumber; i++) {
1657
1658 processTetrahedron<dataTypeU, dataTypeV>(
1659 i, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1660 }
1661
1662 return 0;
1663}
1664
1665template <class dataTypeU, class dataTypeV, typename triangulationType>
1666inline int ttk::FiberSurface::computeSurface(
1667 const triangulationType *const triangulation) {
1668
1669#ifndef TTK_ENABLE_KAMIKAZE
1670 if((!tetNumber_) && (!triangulation))
1671 return -1;
1672 if((!tetList_) && (!triangulation))
1673 return -2;
1674 if(!uField_)
1675 return -3;
1676 if(!vField_)
1677 return -4;
1678 if((!pointSet_) && (!triangulation))
1679 return -5;
1680 if(!polygon_)
1681 return -6;
1682 if(polygonEdgeNumber_ != (SimplexId)polygon_->size())
1683 return -7;
1684 if(!globalVertexList_)
1685 return -8;
1686#endif
1687
1688 Timer t;
1689
1690#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1691 if(!octree_.empty()) {
1692
1693#ifdef TTK_ENABLE_OPENMP
1694#pragma omp parallel for num_threads(threadNumber_)
1695#endif
1696 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1697
1698 computeSurfaceWithOctree<dataTypeU, dataTypeV>(
1699 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1700 }
1701 } else {
1702 // regular extraction (the octree has not been computed)
1703#ifdef TTK_ENABLE_OPENMP
1704#pragma omp parallel for num_threads(threadNumber_)
1705#endif
1706 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1707 computeSurface<dataTypeU, dataTypeV>(
1708 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1709 }
1710 }
1711
1712#else
1713#ifdef TTK_ENABLE_OPENMP
1714#pragma omp parallel for num_threads(threadNumber_)
1715#endif
1716 for(SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1717
1718 computeSurface<dataTypeU, dataTypeV>(
1719 (*polygon_)[i].first, (*polygon_)[i].second, i);
1720 }
1721#endif
1722
1723 finalize<dataTypeU, dataTypeV>(pointSnapping_, false, false, false);
1724
1725 this->printMsg("Extracted", 1.0, t.getElapsedTime(), this->threadNumber_);
1726
1727 return 0;
1728}
1729
1730#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1731template <class dataTypeU, class dataTypeV, typename triangulationType>
1732inline int ttk::FiberSurface::computeSurfaceWithOctree(
1733 const std::pair<double, double> &rangePoint0,
1734 const std::pair<double, double> &rangePoint1,
1735 const triangulationType *const triangulation,
1736 const SimplexId &polygonEdgeId) const {
1737
1738#ifndef TTK_ENABLE_KAMIKAZE
1739 if((!tetNumber_) && (!triangulation))
1740 return -1;
1741 if((!tetList_) && (!triangulation))
1742 return -2;
1743 if(!uField_)
1744 return -3;
1745 if(!vField_)
1746 return -4;
1747 if((!pointSet_) && (!triangulation))
1748 return -5;
1749 if(!polygonEdgeNumber_)
1750 return -6;
1751 if(!globalVertexList_)
1752 return -7;
1753#endif
1754
1755 std::vector<SimplexId> tetList;
1756 octree_.rangeSegmentQuery(rangePoint0, rangePoint1, tetList);
1757
1758#ifdef TTK_ENABLE_OPENMP
1759#pragma omp parallel for num_threads(threadNumber_)
1760#endif
1761 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
1762 processTetrahedron<dataTypeU, dataTypeV>(
1763 tetList[i], rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1764 }
1765
1766 return 0;
1767}
1768#endif
1769
1770template <class dataTypeU, class dataTypeV>
1771int ttk::FiberSurface::finalize(const bool &mergeDuplicatedVertices,
1772 const bool &removeSmallEdges,
1773 const bool &edgeFlips,
1774 const bool &intersectionRemesh) {
1775
1776 // make only one vertex list
1777 SimplexId fiberSurfaceVertexNumber = 0;
1778 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1779 fiberSurfaceVertexNumber += (*polygonEdgeVertexLists_[i]).size();
1780 }
1781
1782 (*globalVertexList_).resize(fiberSurfaceVertexNumber);
1783 fiberSurfaceVertexNumber = 0;
1784 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1785 for(SimplexId j = 0; j < (SimplexId)polygonEdgeVertexLists_[i]->size();
1786 j++) {
1787 (*polygonEdgeVertexLists_[i])[j].polygonEdgeId_ = i;
1788 (*polygonEdgeVertexLists_[i])[j].localId_ = j;
1789 (*polygonEdgeVertexLists_[i])[j].globalId_ = fiberSurfaceVertexNumber;
1790 (*globalVertexList_)[fiberSurfaceVertexNumber]
1791 = (*polygonEdgeVertexLists_[i])[j];
1792 fiberSurfaceVertexNumber++;
1793 }
1794 }
1795 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
1796 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
1797 j++) {
1798 for(int k = 0; k < 3; k++) {
1799 (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k]
1800 = (*polygonEdgeVertexLists_[i])[(*polygonEdgeTriangleLists_[i])[j]
1801 .vertexIds_[k]]
1802 .globalId_;
1803 }
1804 }
1805 }
1806
1807 // NOTE:
1808 // 1) in a first step, make everything work in a POST process.
1809 // this is what really matters for tet gen.
1810 // 2) probably come back and apply in a pre-process (more degenerate
1811 // intersection cases).
1812 if(intersectionRemesh) {
1813 remeshIntersections<dataTypeU, dataTypeV>();
1814 }
1815
1816 if((mergeDuplicatedVertices) || (removeSmallEdges)) {
1817 // we need to have a complex to perform the edge collapses
1818 mergeVertices(pointSnappingThreshold_);
1819 }
1820
1821 if(edgeFlips)
1822 flipEdges();
1823
1824 if(removeSmallEdges)
1825 mergeEdges(edgeCollapseThreshold_);
1826
1827 // now we can release the memory for the threaded vertices
1828 for(SimplexId i = 0; i < (SimplexId)polygonEdgeVertexLists_.size(); i++) {
1829 polygonEdgeVertexLists_[i]->clear();
1830 }
1831
1832 return 0;
1833}
1834
1835template <class dataTypeU, class dataTypeV, typename triangulationType>
1836inline int ttk::FiberSurface::processTetrahedron(
1837 const SimplexId &tetId,
1838 const std::pair<double, double> &rangePoint0,
1839 const std::pair<double, double> &rangePoint1,
1840 const triangulationType *const triangulation,
1841 const SimplexId &polygonEdgeId) const {
1842
1843 double rangeEdge[2];
1844 rangeEdge[0] = rangePoint0.first - rangePoint1.first;
1845 rangeEdge[1] = rangePoint0.second - rangePoint1.second;
1846
1847 double rangeNormal[2];
1848 rangeNormal[0] = -rangeEdge[1];
1849 rangeNormal[1] = rangeEdge[0];
1850
1851 const double prec_dbl = Geometry::powInt(10.0, -DBL_DIG);
1852
1853 // 1. compute the distance to the range line carrying the saddleEdge
1854 SimplexId upperNumber = 0;
1855 SimplexId lowerNumber = 0;
1856 SimplexId equalVertexLocalId = -1;
1857 double d[4];
1858 for(int i = 0; i < 4; i++) {
1859
1860 SimplexId vertexId = 0;
1861 if(!triangulation) {
1862 vertexId = tetList_[5 * tetId + 1 + i];
1863 } else {
1864 triangulation->getCellVertex(tetId, i, vertexId);
1865 }
1866
1867 double projectedVertex[2];
1868 projectedVertex[0] = ((const dataTypeU *)uField_)[vertexId];
1869 projectedVertex[1] = ((const dataTypeV *)vField_)[vertexId];
1870
1871 double vertexRangeEdge[2];
1872 vertexRangeEdge[0] = projectedVertex[0] - rangePoint0.first;
1873 vertexRangeEdge[1] = projectedVertex[1] - rangePoint0.second;
1874
1875 d[i] = vertexRangeEdge[0] * rangeNormal[0]
1876 + vertexRangeEdge[1] * rangeNormal[1];
1877
1878 if(fabs(d[i]) < prec_dbl)
1879 d[i] = 0;
1880
1881 if(d[i] > 0)
1882 upperNumber++;
1883 if(d[i] < 0)
1884 lowerNumber++;
1885
1886 if(d[i] == 0) {
1887 equalVertexLocalId = i;
1888 }
1889 }
1890
1891 // 2. compute the base triangle(s)
1892 if(!((upperNumber == 0) || (lowerNumber == 0))) {
1893
1894 // the fiber surface is passing through this tetrahedron.
1895 std::array<SimplexId, 2> triangleEdgeNumbers{};
1896 std::array<std::array<SimplexId, 3>, 2> triangleEdges{
1897 {{-1, -1, -1}, {-1, -1, -1}}};
1898
1899 // implicit edge encoding
1900 // 0: O-1
1901 // 1: 0-2
1902 // 2: 0-3
1903 // 3: 3-1 [order!]
1904 // 4: 2-1 [order!]
1905 // 5: 2-3
1906 SimplexId edgeCounter = 0;
1907 for(int i = 0; i < 4; i++) {
1908
1909 SimplexId jStart = i + 1;
1910 SimplexId jEnd = 4;
1911 SimplexId jStep = 1;
1912
1913 if(i == 1) {
1914 // special ordering here
1915 // (to facilitate the creation of valid base triangles)
1916 // any two consecutive edges shall share a vertex
1917 jStart = 3;
1918 jEnd = i;
1919 jStep = -1;
1920 }
1921
1922 for(SimplexId j = jStart; j != jEnd; j += jStep) {
1923
1924 if(((d[i] > 0) && (d[j] < 0)) || ((d[i] < 0) && (d[j] > 0))) {
1925
1926 // the edge is crossed by a base triangle
1927 if(triangleEdgeNumbers[0] == 3) {
1928 triangleEdges[1][triangleEdgeNumbers[1]] = edgeCounter;
1929 triangleEdgeNumbers[1]++;
1930 } else {
1931 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1932 triangleEdgeNumbers[0]++;
1933 }
1934 }
1935
1936 if((d[i] == 0) && (d[j] == 0)) {
1937 // special case of a seed tet containing the jacobi edge.
1938 // the entire edge is on the fiber surface.
1939 // let's put this edge twice.
1940 // NOTE: in such a case, we're producing only one triangle.
1941 // so we just need to add that to the first triangle.
1942 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1943 triangleEdgeNumbers[0]++;
1944
1945 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1946 triangleEdgeNumbers[0]++;
1947 }
1948
1949 edgeCounter++;
1950 }
1951 }
1952
1953 // post-process in the case of a second base triangle
1954 if(triangleEdges[1][0] != -1) {
1955 if(triangleEdges[1][1] == -1) {
1956 // we need to consider the edges of the first triangle which share
1957 // a vertex with the second triangle's edge.
1958 // given the ordering, the following edge pairs are forbidden:
1959 // (opposite edges of the tetrahedron)
1960 // 0/5
1961 // 1/3
1962 // 2/4
1963 SimplexId forbiddenEdge = -1;
1964 switch(triangleEdges[1][0]) {
1965 case 0:
1966 forbiddenEdge = 5;
1967 break;
1968 case 1:
1969 forbiddenEdge = 3;
1970 break;
1971 case 2:
1972 forbiddenEdge = 4;
1973 break;
1974 case 3:
1975 forbiddenEdge = 1;
1976 break;
1977 case 4:
1978 forbiddenEdge = 2;
1979 break;
1980 case 5:
1981 forbiddenEdge = 0;
1982 break;
1983 }
1984 for(SimplexId i = 0; i < (SimplexId)triangleEdges[0].size(); i++) {
1985 if(triangleEdges[0][i] != forbiddenEdge) {
1986 if(triangleEdges[1][1] != -1) {
1987 triangleEdges[1][2] = triangleEdges[0][i];
1988 break;
1989 } else {
1990 triangleEdges[1][1] = triangleEdges[0][i];
1991 }
1992 }
1993 }
1994 }
1995 }
1996
1997 // post-process in case of exactly one vertex on the jacobi edge
1998 for(int i = 0; i < 2; i++) {
1999 if(triangleEdges[i][0] != -1) {
2000 // the jacobi vertex has to be the last one
2001 if(triangleEdges[i][2] == -1) {
2002 // add whatever edge connected to equalVertexLocalId
2003 switch(equalVertexLocalId) {
2004 case 0:
2005 case 1:
2006 triangleEdges[i][2] = 0;
2007 break;
2008 case 2:
2009 case 3:
2010 triangleEdges[i][2] = 5;
2011 break;
2012 }
2013 }
2014 }
2015 }
2016
2017 // 3. crop the resulting triangles to the saddleEdge
2018 // for each edge recorded previously, get the u,v coordinates of the
2019 // intersection point.
2020 // take its barycentric coordinates on the saddle edge.
2021 // see if it lies in it (in between 0 and 1)
2022 // figure 7 of the paper
2023 double d0, d1;
2024 std::pair<double, double> uv0, uv1;
2025 std::array<std::pair<double, double>, 3> uv{};
2026 std::array<double, 3> t{};
2027
2028 SimplexId createdVertices = 0;
2029
2030 for(int i = 0; i < 2; i++) {
2031 if(triangleEdges[i][0] != -1) {
2032 // this is a valid triangle, let's go ahead.
2033
2034 SimplexId lowerVertexNumber = 0;
2035 SimplexId upperVertexNumber = 0;
2036 SimplexId greyVertexNumber = 0;
2037 // iterate over the edges and compute the edge intersection (range)
2038 for(SimplexId j = 0; j < (SimplexId)triangleEdges[i].size(); j++) {
2039
2040 SimplexId vertexId0 = 0, vertexId1 = 0;
2041 if(triangulation) {
2042 triangulation->getCellVertex(
2043 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j]], vertexId0);
2044 triangulation->getCellVertex(
2045 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1],
2046 vertexId1);
2047 } else {
2048 vertexId0
2049 = tetList_[5 * tetId + 1
2050 + edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2051 vertexId1
2052 = tetList_[5 * tetId + 1
2053 + edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2054 }
2055
2056 if((j < (SimplexId)triangleEdges[i].size() - 1)
2057 && (triangleEdges[i][j] == triangleEdges[i][j + 1])) {
2058
2059 // special case of a jacobi edge
2060 uv[j].first = ((const dataTypeU *)uField_)[vertexId0];
2061 uv[j].second = ((const dataTypeV *)vField_)[vertexId0];
2062
2063 uv[j + 1].first = ((const dataTypeU *)uField_)[vertexId1];
2064 uv[j + 1].second = ((const dataTypeV *)vField_)[vertexId1];
2065
2066 } else if((!j)
2067 || ((j)
2068 && (triangleEdges[i][j] != triangleEdges[i][j - 1]))) {
2069
2070 // regular intersection case
2071 d0 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2072 uv0.first = ((const dataTypeU *)uField_)[vertexId0];
2073 uv0.second = ((const dataTypeV *)vField_)[vertexId0];
2074
2075 d1 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2076 uv1.first = ((const dataTypeU *)uField_)[vertexId1];
2077 uv1.second = ((const dataTypeV *)vField_)[vertexId1];
2078
2079 uv[j].first
2080 = uv0.first + (d0 / (d0 - d1)) * (uv1.first - uv0.first);
2081 uv[j].second
2082 = uv0.second + (d0 / (d0 - d1)) * (uv1.second - uv0.second);
2083 }
2084
2085 // now determine the line parameterization of this intersection
2086 // point on the saddle edge
2087 if(fabs(rangePoint1.first - rangePoint0.first)
2088 > fabs(rangePoint1.second - rangePoint0.second)) {
2089 t[j] = (uv[j].first - rangePoint0.first)
2090 / (rangePoint1.first - rangePoint0.first);
2091 } else {
2092 t[j] = (uv[j].second - rangePoint0.second)
2093 / (rangePoint1.second - rangePoint0.second);
2094 }
2095
2096 if((t[j] <= 1) && (t[j] >= 0))
2097 greyVertexNumber++;
2098 else if(t[j] < 0)
2099 lowerVertexNumber++;
2100 else
2101 upperVertexNumber++;
2102 }
2103 // at this point, we know the uv coordinates (and the edge param)
2104 // of each vertex of the base triangle.
2105 // we can proceed with the cropping
2106
2107 // 4. triangulate the result
2108 if(greyVertexNumber == 3) {
2109 createdVertices += computeCase0<dataTypeU, dataTypeV>(
2110 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2111 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2112 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2113 triangulation);
2114 } else if(lowerVertexNumber == 3 || upperVertexNumber == 3) {
2115 // well do nothing (empty triangle)
2116 } else if((lowerVertexNumber == 1) && (upperVertexNumber == 1)
2117 && (greyVertexNumber == 1)) {
2118 createdVertices += computeCase1<dataTypeU, dataTypeV>(
2119 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2120 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2121 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2122 triangulation);
2123 } else if(((lowerVertexNumber == 2) && (upperVertexNumber == 1))
2124 || ((lowerVertexNumber == 1) && (upperVertexNumber == 2))) {
2125 createdVertices += computeCase2<dataTypeU, dataTypeV>(
2126 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2127 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2128 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2129 triangulation);
2130 } else if((greyVertexNumber == 1)
2131 && ((lowerVertexNumber == 2) || (upperVertexNumber == 2))) {
2132 createdVertices += computeCase3<dataTypeU, dataTypeV>(
2133 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2134 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2135 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2136 triangulation);
2137 } else if(((greyVertexNumber == 2))
2138 && ((lowerVertexNumber == 1) || (upperVertexNumber == 1))) {
2139 createdVertices += computeCase4<dataTypeU, dataTypeV>(
2140 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2141 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2142 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2143 triangulation);
2144 }
2145 }
2146 }
2147
2148 // in case of 2 triangles, remesh locally
2149 // NOTE: here we just snap vertices together if they are colinear
2150 // the mergeFiberSurfaces() function will take care of removing any
2151 // duplicate
2152 if((triangleEdges[1][0] != -1) && (createdVertices > 3)) {
2153
2154 std::vector<SimplexId> createdVertexList(createdVertices);
2155 for(SimplexId i = 0; i < (SimplexId)createdVertices; i++) {
2156 createdVertexList[i]
2157 = polygonEdgeVertexLists_[polygonEdgeId]->size() - 1 - i;
2158 }
2159
2160 std::vector<bool> snappedVertices(createdVertices, false);
2161 std::vector<SimplexId> colinearVertices;
2162 colinearVertices.reserve(createdVertices);
2163
2164 for(SimplexId i = 0; i < createdVertices; i++) {
2165 colinearVertices.clear();
2166
2167 if(!snappedVertices[i]) {
2168 colinearVertices.push_back(i);
2169 for(SimplexId j = 0; j < createdVertices; j++) {
2170 if((i != j) && (!snappedVertices[j])) {
2171 // not the same vertex
2172 // not snapped already
2173
2174 if((*polygonEdgeVertexLists_[polygonEdgeId])[createdVertexList[i]]
2175 .t_
2176 == (*polygonEdgeVertexLists_[polygonEdgeId])
2177 [createdVertexList[j]]
2178 .t_) {
2179 colinearVertices.push_back(j);
2180 }
2181 }
2182 }
2183 }
2184
2185 if((colinearVertices.size() == 4) || (colinearVertices.size() == 3)) {
2186 // 3 co-linear vertices with 1 duplicate
2187
2188 // we just need to find the pair of duplicates and snap both of
2189 // them to another vertex
2190 std::pair<SimplexId, SimplexId> minPair;
2191 double minDistance = -1;
2192 for(SimplexId j = 0; j < (SimplexId)colinearVertices.size(); j++) {
2193 for(SimplexId k = 0; k < (SimplexId)colinearVertices.size(); k++) {
2194 if(j != k) {
2195
2196 double const distance = Geometry::distance(
2197 (*polygonEdgeVertexLists_[polygonEdgeId])
2198 [createdVertexList[colinearVertices[j]]]
2199 .p_.data(),
2200 (*polygonEdgeVertexLists_[polygonEdgeId])
2201 [createdVertexList[colinearVertices[k]]]
2202 .p_.data());
2203
2204 // bool basePointSnap = true;
2205 // for(int l = 0; l < 3; l++){
2206 // if((*polygonEdgeVertexLists_[polygonEdgeId])[
2207 // createdVertexList[colinearVertices[j]]].p_[l]
2208 // !=
2209 // (*polygonEdgeVertexLists_[polygonEdgeId])[
2210 // createdVertexList[colinearVertices[k]]].p_[l]){
2211 // basePointSnap = false;
2212 // break;
2213 // }
2214 // }
2215
2216 // if(!basePointSnap){
2217 if((minDistance < 0) || (distance < minDistance)) {
2218 minDistance = distance;
2219 minPair.first = j;
2220 minPair.second = k;
2221 }
2222 // }
2223 }
2224 }
2225 }
2226 if((minDistance != -1) && (minDistance < prec_dbl)) {
2227 // if((minDistance != -1)&&(minDistance <
2228 // pointSnappingThreshold_)){
2229 // snap them to another colinear vertex
2230 for(SimplexId j = 0; j < (SimplexId)colinearVertices.size(); j++) {
2231 if((j != minPair.first) && (j != minPair.second)) {
2232 // snap minPair.first to j
2233
2234 for(int k = 0; k < 3; k++) {
2235 (*polygonEdgeVertexLists_[polygonEdgeId])
2236 [createdVertexList[colinearVertices[minPair.first]]]
2237 .p_[k]
2238 = (*polygonEdgeVertexLists_[polygonEdgeId])
2239 [createdVertexList[colinearVertices[j]]]
2240 .p_[k];
2241 }
2242 (*polygonEdgeVertexLists_[polygonEdgeId])
2243 [createdVertexList[colinearVertices[minPair.first]]]
2244 .uv_
2245 = (*polygonEdgeVertexLists_[polygonEdgeId])
2246 [createdVertexList[colinearVertices[j]]]
2247 .uv_;
2248
2249 (*polygonEdgeVertexLists_[polygonEdgeId])
2250 [createdVertexList[colinearVertices[minPair.first]]]
2251 .t_
2252 = (*polygonEdgeVertexLists_[polygonEdgeId])
2253 [createdVertexList[colinearVertices[j]]]
2254 .t_;
2255
2256 (*polygonEdgeVertexLists_[polygonEdgeId])
2257 [createdVertexList[colinearVertices[minPair.first]]]
2258 .isBasePoint_
2259 = (*polygonEdgeVertexLists_[polygonEdgeId])
2260 [createdVertexList[colinearVertices[j]]]
2261 .isBasePoint_;
2262
2263 (*polygonEdgeVertexLists_[polygonEdgeId])
2264 [createdVertexList[colinearVertices[minPair.first]]]
2265 .isIntersectionPoint_
2266 = (*polygonEdgeVertexLists_[polygonEdgeId])
2267 [createdVertexList[colinearVertices[j]]]
2268 .isIntersectionPoint_;
2269 if((*polygonEdgeVertexLists_[polygonEdgeId])
2270 [createdVertexList[colinearVertices[j]]]
2271 .meshEdge_.first
2272 != -1) {
2273 (*polygonEdgeVertexLists_[polygonEdgeId])
2274 [createdVertexList[colinearVertices[minPair.first]]]
2275 .meshEdge_
2276 = (*polygonEdgeVertexLists_[polygonEdgeId])
2277 [createdVertexList[colinearVertices[j]]]
2278 .meshEdge_;
2279 }
2280
2281 // snap minPair.second to j
2282 for(int k = 0; k < 3; k++) {
2283 (*polygonEdgeVertexLists_[polygonEdgeId])
2284 [createdVertexList[colinearVertices[minPair.second]]]
2285 .p_[k]
2286 = (*polygonEdgeVertexLists_[polygonEdgeId])
2287 [createdVertexList[colinearVertices[j]]]
2288 .p_[k];
2289 }
2290 (*polygonEdgeVertexLists_[polygonEdgeId])
2291 [createdVertexList[colinearVertices[minPair.second]]]
2292 .uv_
2293 = (*polygonEdgeVertexLists_[polygonEdgeId])
2294 [createdVertexList[colinearVertices[j]]]
2295 .uv_;
2296
2297 (*polygonEdgeVertexLists_[polygonEdgeId])
2298 [createdVertexList[colinearVertices[minPair.second]]]
2299 .t_
2300 = (*polygonEdgeVertexLists_[polygonEdgeId])
2301 [createdVertexList[colinearVertices[j]]]
2302 .t_;
2303
2304 (*polygonEdgeVertexLists_[polygonEdgeId])
2305 [createdVertexList[colinearVertices[minPair.second]]]
2306 .isBasePoint_
2307 = (*polygonEdgeVertexLists_[polygonEdgeId])
2308 [createdVertexList[colinearVertices[j]]]
2309 .isBasePoint_;
2310
2311 (*polygonEdgeVertexLists_[polygonEdgeId])
2312 [createdVertexList[colinearVertices[minPair.second]]]
2313 .isIntersectionPoint_
2314 = (*polygonEdgeVertexLists_[polygonEdgeId])
2315 [createdVertexList[colinearVertices[j]]]
2316 .isIntersectionPoint_;
2317 if((*polygonEdgeVertexLists_[polygonEdgeId])
2318 [createdVertexList[colinearVertices[j]]]
2319 .meshEdge_.first
2320 != -1) {
2321 (*polygonEdgeVertexLists_[polygonEdgeId])
2322 [createdVertexList[colinearVertices[minPair.second]]]
2323 .meshEdge_
2324 = (*polygonEdgeVertexLists_[polygonEdgeId])
2325 [createdVertexList[colinearVertices[j]]]
2326 .meshEdge_;
2327 }
2328
2329 snappedVertices[colinearVertices[minPair.first]] = true;
2330 snappedVertices[colinearVertices[minPair.second]] = true;
2331
2332 // we're done
2333 break;
2334 }
2335 }
2336 }
2337 }
2338 }
2339 }
2340
2341 return createdVertices;
2342 }
2343
2344 return 0;
2345}
2346
2347template <class dataTypeU, class dataTypeV>
2348inline int ttk::FiberSurface::remeshIntersections() const {
2349
2350#ifndef TTK_ENABLE_KAMIKAZE
2351 if(!tetNumber_)
2352 return -1;
2353#endif
2354
2355 // Algorithm
2356 // 1) loop over each triangle and mark the containing tetrahedron
2357 // 2) for each tet, in parallel, check for pairwise triangle intersections
2358 // (based on the triangles' range projection)
2359 // Given a pair of intersecting triangles:
2360 // 3) given the point of intersection, take one of its coordinates (u or v)
2361 // and compute an iso-contour I on both triangles
2362 // 4) express the barycentric coordinates of I in both triangles, 2 cases:
2363 // a) the intersection is the entire fiber (old code)
2364 // b) the intersection is a segment of fiber
2365
2366 // NOTE:
2367 // the topological aspect of the code is OK.
2368 // if any bug, it's very likely to be geometry accuracy related.
2369
2370 std::vector<std::vector<IntersectionTriangle>> tetIntersections(tetNumber_);
2371 std::vector<std::vector<Vertex>> tetNewVertices(tetNumber_);
2372
2373 // fill the information prior to the parallel pass
2374 for(SimplexId i = 0; i < (SimplexId)polygonEdgeTriangleLists_.size(); i++) {
2375 for(SimplexId j = 0; j < (SimplexId)polygonEdgeTriangleLists_[i]->size();
2376 j++) {
2377
2378 SimplexId const tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
2379
2380 tetIntersections[tetId].emplace_back();
2381
2382 tetIntersections[tetId].back().caseId_
2383 = (*polygonEdgeTriangleLists_[i])[j].caseId_;
2384 tetIntersections[tetId].back().polygonEdgeId_ = i;
2385 tetIntersections[tetId].back().triangleId_ = j;
2386 for(int k = 0; k < 3; k++) {
2387 tetIntersections[tetId].back().vertexIds_[k]
2388 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
2389 tetIntersections[tetId].back().uv_[k]
2390 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2391 .uv_;
2392 tetIntersections[tetId].back().t_[k]
2393 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2394 .t_;
2395 for(int l = 0; l < 3; l++) {
2396 tetIntersections[tetId].back().p_[k][l]
2397 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2398 .p_[l];
2399 }
2400 }
2401 tetIntersections[tetId].back().intersection_.first = -DBL_MAX;
2402 tetIntersections[tetId].back().intersection_.second = -DBL_MAX;
2403 }
2404 }
2405
2406 std::vector<SimplexId> tetList;
2407 for(SimplexId i = 0; i < (SimplexId)tetIntersections.size(); i++) {
2408 if(tetIntersections[i].size() > 1)
2409 tetList.push_back(i);
2410 }
2411
2412#ifdef TTK_ENABLE_OPENMP
2413#pragma omp parallel for num_threads(threadNumber_)
2414#endif
2415 for(SimplexId i = 0; i < (SimplexId)tetList.size(); i++) {
2416 SimplexId const tetId = tetList[i];
2417
2418 // pre-process by merging nearby vertices...?
2419
2420 SimplexId newTriangleNumber = 1;
2421 SimplexId newVertexNumber = 1;
2422
2423 for(SimplexId j = 0; j < (SimplexId)tetIntersections[tetId].size(); j++) {
2424
2425 if(j > 1000) {
2426 this->printWrn("Preventing an infinite loop!");
2427 this->printWrn("More than 1000 re-meshed triangles in tet #"
2428 + std::to_string(tetId));
2429 this->printWrn("Extra-thin triangles keep on intersecting?!");
2430 break;
2431 }
2432
2433 // if we re-mesh, we add new triangles at the end of the list.
2434 // there's no need to check intersections with those.
2435 SimplexId const originalTriangleNumber
2436 = (SimplexId)tetIntersections[tetId].size();
2437
2438 for(SimplexId k = 0; k < originalTriangleNumber; k++) {
2439
2440 SimplexId const polygonEdgeId0
2441 = tetIntersections[tetId][j].polygonEdgeId_;
2442 SimplexId const polygonEdgeId1
2443 = tetIntersections[tetId][k].polygonEdgeId_;
2444
2445 if((j != k) && (polygonEdgeId0 != polygonEdgeId1)) {
2446 // cases 3, 4 and 6 of the fiber surface table (multiple triangle
2447 // per tet given a single edge). we don't need to re-mesh that.
2448
2449 // grab the range projection for the triangle j
2450 std::pair<double, double> edge0point0, edge0point1;
2451
2452 getTriangleRangeExtremities(
2453 tetId, j, tetIntersections, edge0point0, edge0point1);
2454
2455 // now do the same thing for the triangle k
2456 std::pair<double, double> edge1point0, edge1point1;
2457
2458 getTriangleRangeExtremities(
2459 tetId, k, tetIntersections, edge1point0, edge1point1);
2460
2461 // compute the intersection
2462 std::pair<double, double> intersection;
2463 bool const hasIntersection = Geometry::computeSegmentIntersection(
2464 edge0point0.first, edge0point0.second, edge0point1.first,
2465 edge0point1.second, edge1point0.first, edge1point0.second,
2466 edge1point1.first, edge1point1.second, intersection.first,
2467 intersection.second);
2468
2469 if((hasIntersection)
2470 // check if that intersection has been registered before
2471 // in the end, only one intersection per triangle, no matter what
2472 /*&&(((fabs(tetIntersections[tetId][j].intersection_.first
2473 - intersection.first) > Geometry::powIntTen(-FLT_DIG))
2474 ||(fabs(tetIntersections[tetId][j].intersection_.second
2475 - intersection.second) > Geometry::powIntTen(-FLT_DIG)))
2476 &&((fabs(tetIntersections[tetId][k].intersection_.first
2477 - intersection.first) > Geometry::powIntTen(-FLT_DIG))
2478 ||(fabs(tetIntersections[tetId][k].intersection_.second
2479 - intersection.second) > Geometry::powIntTen(-FLT_DIG))))*/){
2480
2481 computeTriangleIntersection(tetId, j, k, polygonEdgeId0,
2482 polygonEdgeId1, intersection,
2483 newVertexNumber, newTriangleNumber,
2484 tetIntersections, tetNewVertices);
2485 }
2486 }
2487 }
2488 }
2489 }
2490
2491 // now copy the new vertices
2492 for(SimplexId i = 0; i < (SimplexId)tetNewVertices.size(); i++) {
2493 for(SimplexId j = 0; j < (SimplexId)tetNewVertices[i].size(); j++) {
2494
2495 SimplexId const localId = (*globalVertexList_).size();
2496 tetNewVertices[i][j].localId_ = localId;
2497 (*globalVertexList_).push_back(tetNewVertices[i][j]);
2498 }
2499 }
2500
2501 for(SimplexId i = 0; i < (SimplexId)tetIntersections.size(); i++) {
2502 for(SimplexId j = 0; j < (SimplexId)tetIntersections[i].size(); j++) {
2503
2504 if(((tetIntersections[i][j].intersection_.first != -DBL_MAX)
2505 && (tetIntersections[i][j].intersection_.second != -DBL_MAX))
2506 || (tetIntersections[i][j].triangleId_ < 0)) {
2507
2508 SimplexId triangleId = tetIntersections[i][j].triangleId_;
2509
2510 if(triangleId < 0) {
2511
2512 // this is a new triangle
2513 triangleId = (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2514 .polygonEdgeId_])
2515 .size();
2516 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2517 .resize(triangleId + 1);
2518 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2519 .back()
2520 .tetId_
2521 = i;
2522 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2523 .back()
2524 .caseId_
2525 = tetIntersections[i][j].caseId_;
2526 }
2527
2528 for(int k = 0; k < 3; k++) {
2529
2530 SimplexId vertexId = tetIntersections[i][j].vertexIds_[k];
2531
2532 if(vertexId < 0) {
2533 // newly created vertex
2534 vertexId = tetNewVertices[i][-(vertexId + 1)].localId_;
2535 }
2536 (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2537 .polygonEdgeId_])[triangleId]
2538 .vertexIds_[k]
2539 = vertexId;
2540 }
2541 }
2542 }
2543 }
2544
2545 return 0;
2546}
ttk::AbstractTriangulation
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
Definition AbstractTriangulation.h:60
ttk::Debug
Minimalist debugging class.
Definition Debug.h:88
ttk::Debug::debugLevel_
int debugLevel_
Definition Debug.h:379
ttk::Debug::setDebugLevel
virtual int setDebugLevel(const int &debugLevel)
Definition Debug.cpp:147
ttk::FiberSurface
TTK processing package that computes fiber surfaces.
Definition FiberSurface.h:49
ttk::FiberSurface::polygonEdgeVertexLists_
std::vector< std::vector< Vertex > * > polygonEdgeVertexLists_
Definition FiberSurface.h:497
ttk::FiberSurface::vField_
const void * vField_
Definition FiberSurface.h:483
ttk::FiberSurface::setGlobalVertexList
int setGlobalVertexList(std::vector< Vertex > *globalList)
Definition FiberSurface.h:126
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::edgeImplicitEncoding_
std::array< SimplexId, 12 > edgeImplicitEncoding_
Definition FiberSurface.h:487
ttk::FiberSurface::pointSnappingThreshold_
double pointSnappingThreshold_
Definition FiberSurface.h:491
ttk::FiberSurface::edgeCollapseThreshold_
double edgeCollapseThreshold_
Definition FiberSurface.h:490
ttk::FiberSurface::polygonEdgeTriangleLists_
std::vector< std::vector< Triangle > * > polygonEdgeTriangleLists_
Definition FiberSurface.h:498
ttk::FiberSurface::tetNeighbors_
const std::vector< std::vector< SimplexId > > * tetNeighbors_
Definition FiberSurface.h:486
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::setInputField
int setInputField(const void *uField, const void *vField)
Definition FiberSurface.h:131
ttk::FiberSurface::setTetList
int setTetList(const SimplexId *tetList)
Definition FiberSurface.h:173
ttk::FiberSurface::computeCase3
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
Definition FiberSurface.h:1165
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::finalize
int finalize(const bool &mergeDuplicatedVertices=false, const bool &removeSmallEdges=false, const bool &edgeFlips=false, const bool &intersectionRemesh=false)
Definition FiberSurface.h:1771
ttk::FiberSurface::polygon_
const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > * polygon_
Definition FiberSurface.h:494
ttk::FiberSurface::flipEdges
int flipEdges() const
Definition FiberSurface.cpp:600
ttk::FiberSurface::computeCase4
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
Definition FiberSurface.h:1306
ttk::FiberSurface::computeCase1
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
ttk::FiberSurface::computeSurface
int computeSurface(const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
Definition FiberSurface.h:1624
ttk::FiberSurface::remeshIntersections
int remeshIntersections() const
Definition FiberSurface.h:2348
ttk::FiberSurface::computeCase0
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
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::processTetrahedron
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
Definition FiberSurface.h:1836
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::setPolygonEdgeNumber
int setPolygonEdgeNumber(const SimplexId &polygonEdgeNumber)
Definition FiberSurface.h:166
ttk::FiberSurface::mergeVertices
int mergeVertices(const double &distanceThreshold) const
Definition FiberSurface.cpp:1276
ttk::FiberSurface::setTetNumber
int setTetNumber(const SimplexId &tetNumber)
Definition FiberSurface.h:184
ttk::FiberSurface::computeCase2
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:991
ttk::FiberSurface::globalVertexList_
std::vector< Vertex > * globalVertexList_
Definition FiberSurface.h:496
ttk::FiberSurface::pointNumber_
SimplexId pointNumber_
Definition FiberSurface.h:482
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::setPointSet
int setPointSet(const float *pointSet)
Definition FiberSurface.h:154
ttk::FiberSurface::setPointNumber
int setPointNumber(const SimplexId &number)
Definition FiberSurface.h:149
ttk::FiberSurface::setPointMerging
int setPointMerging(const bool &onOff)
Definition FiberSurface.h:139
ttk::FiberSurface::setPointMergingThreshold
int setPointMergingThreshold(const double &threshold)
Definition FiberSurface.h:144
ttk::FiberSurface::computeBaseTriangle
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
ttk::FiberSurface::setPolygon
int setPolygon(const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > *polygon)
Definition FiberSurface.h:159
ttk::FiberSurface::tetList_
const SimplexId * tetList_
Definition FiberSurface.h:485
ttk::FiberSurface::uField_
const void * uField_
Definition FiberSurface.h:483
ttk::FiberSurface::polygonEdgeNumber_
SimplexId polygonEdgeNumber_
Definition FiberSurface.h:482
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::setVertexList
int setVertexList(const SimplexId &polygonEdgeId, std::vector< Vertex > *vertexList)
Definition FiberSurface.h:209
ttk::FiberSurface::setTetNeighbors
int setTetNeighbors(const std::vector< std::vector< SimplexId > > *tetNeighbors)
Definition FiberSurface.h:179
ttk::FiberSurface::setTriangleList
int setTriangleList(const SimplexId &polygonEdgeId, std::vector< Triangle > *triangleList)
Definition FiberSurface.h:189
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::FiberSurface::preconditionTriangulation
void preconditionTriangulation(AbstractTriangulation *triangulation)
Definition FiberSurface.h:202
ttk::FiberSurface::computeContour
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
Definition FiberSurface.h:1477
ttk::RangeDrivenOctree
TTK optional package for octree based range queries in bivariate volumetric data.
Definition RangeDrivenOctree.h:27
ttk::Timer::getElapsedTime
double getElapsedTime()
Definition Timer.h:15
ttk::Geometry::computeSegmentIntersection
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:248
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::isTriangleColinear
bool isTriangleColinear(const T *p0, const T *p1, const T *p2, const T *tolerance=nullptr)
Definition Geometry.cpp:466
ttk::Geometry::powIntTen
T powIntTen(const int n)
Compute the nth power of ten.
Definition Geometry.h:405
ttk::Geometry::powInt
T powInt(const T val, const int n)
Definition Geometry.h:383
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::IntersectionTriangle::p_
std::array< std::array< double, 3 >, 3 > p_
Definition FiberSurface.h:231
ttk::FiberSurface::IntersectionTriangle::intersection_
std::pair< double, double > intersection_
Definition FiberSurface.h:232
ttk::FiberSurface::IntersectionTriangle::uv_
std::array< std::pair< double, double >, 3 > uv_
Definition FiberSurface.h:229
ttk::FiberSurface::IntersectionTriangle::vertexIds_
std::array< SimplexId, 3 > vertexIds_
Definition FiberSurface.h:228
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
ttk::FiberSurface::Vertex::meshEdge_
std::pair< SimplexId, SimplexId > meshEdge_
Definition FiberSurface.h:57
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