39#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
58 std::array<double, 3>
p_{};
60 std::pair<double, double>
uv_{};
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);
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;
82 template <
class dataTypeU,
class dataTypeV>
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;
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;
95 template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
96 inline int computeSurface(
const triangulationType *
const triangulation);
98#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
99 template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
101 computeSurfaceWithOctree(
const std::pair<double, double> &rangePoint0,
102 const std::pair<double, double> &rangePoint1,
103 const triangulationType *
const triangulation,
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);
113#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
114 inline void flushOctree() {
119 template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
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;
160 const std::vector<std::pair<std::pair<double, double>,
161 std::pair<double, double>>> *polygon) {
190 std::vector<Triangle> *triangleList) {
192#ifndef TTK_ENABLE_KAMIKAZE
193 if((polygonEdgeId >= 0)
204 if(triangulation !=
nullptr) {
210 std::vector<Vertex> *vertexList) {
212#ifndef TTK_ENABLE_KAMIKAZE
213 if((polygonEdgeId >= 0)
229 std::array<std::pair<double, double>, 3>
uv_{};
230 std::array<double, 3>
t_{};
231 std::array<std::array<double, 3>, 3>
p_{};
235 template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
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;
256 template <
class dataTypeU,
class dataTYpeV,
typename triangulationType>
271 const triangulationType *
const triangulation)
const;
273 template <
class dataTypeU,
class dataTYpeV,
typename triangulationType>
288 const triangulationType *
const triangulation)
const;
290 template <
class dataTypeU,
class dataTYpeV,
typename triangulationType>
305 const triangulationType *
const triangulation)
const;
307 template <
class dataTypeU,
class dataTYpeV,
typename triangulationType>
322 const triangulationType *
const triangulation)
const;
324 template <
class dataTypeU,
class dataTYpeV,
typename triangulationType>
339 const triangulationType *
const triangulation)
const;
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,
349 bool &edgeFiber)
const;
357 const std::pair<double, double> &intersection,
360 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
361 std::vector<std::vector<Vertex>> &tetNewVertices)
const;
367 const std::pair<double, double> &intersection,
368 const std::array<double, 3> &pA,
369 const std::array<double, 3> &pB,
373 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
374 std::vector<std::vector<Vertex>> &tetNewVertices)
const;
382 const std::vector<std::vector<Vertex>> &tetNewVertices,
384 std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
385 const std::pair<double, double> *intersection =
nullptr)
const;
390 flipEdges(std::vector<std::pair<SimplexId, SimplexId>> &triangles)
const;
396 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections)
402 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
403 std::pair<double, double> &extremity0,
404 std::pair<double, double> &extremity1)
const;
408 const double *p2)
const;
411 const std::pair<double, double> &uv0,
413 const std::array<double, 3> &p1,
414 const std::pair<double, double> &uv1,
423 const std::vector<std::pair<SimplexId, SimplexId>> &starNeighbors)
const;
433 const std::vector<std::vector<IntersectionTriangle>> &tetIntersections,
434 const std::vector<std::vector<Vertex>> &tetNewVertices,
439 std::array<std::array<double, 3>, 3> points{};
440 for(
int i = 0; i < 3; i++) {
443 vertexId = vertexId1;
445 vertexId = vertexId2;
448 for(
int j = 0; j < 3; j++) {
449 points[i][j] = tetIntersections[tetId][triangleId].p_[vertexId][j];
452 for(
int j = 0; j < 3; j++) {
453 points[i][j] = tetNewVertices[tetId][(-vertexId) - 1].p_[j];
459 points[0].data(), points[1].data(), points[2].data());
462 int mergeEdges(
const double &distanceThreshold)
const;
466 template <
class dataTypeU,
class dataTypeV>
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,
478 const std::vector<std::pair<SimplexId, SimplexId>> &triangles)
const;
488 0, 1, 0, 2, 0, 3, 3, 1, 2, 1, 2, 3};
493 const std::vector<std::pair<std::pair<double, double>,
500#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
506#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
507template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
509 ttk::FiberSurface::buildOctree(
const triangulationType *
const triangulation) {
516 if(octree_.empty()) {
528 octree_.build<dataTypeU, dataTypeV>(triangulation);
535template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
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 {
556 for(
int i = 0; i < 3; i++) {
565 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId0]];
566 vertexId1 = tetList_[5 * tetId + 1
567 + edgeImplicitEncoding_[2 * localEdgeId0 + 1]];
569 triangulation->getCellVertex(
570 tetId, edgeImplicitEncoding_[2 * localEdgeId0], vertexId0);
571 triangulation->getCellVertex(
572 tetId, edgeImplicitEncoding_[2 * localEdgeId0 + 1], vertexId1);
574 basePointProjections[i].first = u0;
575 basePointProjections[i].second = v0;
576 basePointParameterization[i] = t0;
582 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId1]];
583 vertexId1 = tetList_[5 * tetId + 1
584 + edgeImplicitEncoding_[2 * localEdgeId1 + 1]];
586 triangulation->getCellVertex(
587 tetId, edgeImplicitEncoding_[2 * localEdgeId1], vertexId0);
588 triangulation->getCellVertex(
589 tetId, edgeImplicitEncoding_[2 * localEdgeId1 + 1], vertexId1);
591 basePointProjections[i].first = u1;
592 basePointProjections[i].second = v1;
593 basePointParameterization[i] = t1;
599 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId2]];
600 vertexId1 = tetList_[5 * tetId + 1
601 + edgeImplicitEncoding_[2 * localEdgeId2 + 1]];
603 triangulation->getCellVertex(
604 tetId, edgeImplicitEncoding_[2 * localEdgeId2], vertexId0);
605 triangulation->getCellVertex(
606 tetId, edgeImplicitEncoding_[2 * localEdgeId2 + 1], vertexId1);
608 basePointProjections[i].first = u2;
609 basePointProjections[i].second = v2;
610 basePointParameterization[i] = t2;
614 std::array<double, 2> baryCentrics{};
615 std::array<double, 2> p0{}, p1{}, p{};
616 p0[0] = ((
const dataTypeU *)uField_)[vertexId0];
617 p0[1] = ((
const dataTypeV *)vField_)[vertexId0];
618 p1[0] = ((
const dataTypeU *)uField_)[vertexId1];
619 p1[1] = ((
const dataTypeV *)vField_)[vertexId1];
620 p[0] = basePointProjections[i].first;
621 p[1] = basePointProjections[i].second;
623 p0.data(), p1.data(), p.data(), baryCentrics, 2);
625 std::array<float, 3> pA{}, pB{};
627 triangulation->getVertexPoint(vertexId0, pA[0], pA[1], pA[2]);
628 triangulation->getVertexPoint(vertexId1, pB[0], pB[1], pB[2]);
631 for(
int j = 0; j < 3; j++) {
635 c0 = pointSet_[3 * vertexId0 + j];
636 c1 = pointSet_[3 * vertexId1 + j];
642 basePoints[i][j] = baryCentrics[0] * c0 + baryCentrics[1] * c1;
645 if(vertexId0 < vertexId1) {
646 baseEdges[i] = std::pair<SimplexId, SimplexId>(vertexId0, vertexId1);
648 baseEdges[i] = std::pair<SimplexId, SimplexId>(vertexId1, vertexId0);
655template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
671 const triangulationType *
const triangulation)
const {
674 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
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_
684 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[0] = vertexId;
685 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[1]
687 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().vertexIds_[2]
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_
699 for(
int i = 0; i < 3; i++) {
707 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId0]];
708 vertexId1 = tetList_[5 * tetId + 1
709 + edgeImplicitEncoding_[2 * localEdgeId0 + 1]];
711 triangulation->getCellVertex(
712 tetId, edgeImplicitEncoding_[2 * localEdgeId0], vertexId0);
713 triangulation->getCellVertex(
714 tetId, edgeImplicitEncoding_[2 * localEdgeId0 + 1], vertexId1);
716 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u0;
717 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v0;
718 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t0;
724 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId1]];
725 vertexId1 = tetList_[5 * tetId + 1
726 + edgeImplicitEncoding_[2 * localEdgeId1 + 1]];
728 triangulation->getCellVertex(
729 tetId, edgeImplicitEncoding_[2 * localEdgeId1], vertexId0);
730 triangulation->getCellVertex(
731 tetId, edgeImplicitEncoding_[2 * localEdgeId1 + 1], vertexId1);
733 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u1;
734 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v1;
735 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t1;
741 = tetList_[5 * tetId + 1 + edgeImplicitEncoding_[2 * localEdgeId2]];
742 vertexId1 = tetList_[5 * tetId + 1
743 + edgeImplicitEncoding_[2 * localEdgeId2 + 1]];
745 triangulation->getCellVertex(
746 tetId, edgeImplicitEncoding_[2 * localEdgeId2], vertexId0);
747 triangulation->getCellVertex(
748 tetId, edgeImplicitEncoding_[2 * localEdgeId2 + 1], vertexId1);
750 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first = u2;
751 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second = v2;
752 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t2;
756 std::array<double, 2> baryCentrics{};
757 std::array<double, 2> p0{}, p1{}, p{};
758 p0[0] = ((
const dataTypeU *)uField_)[vertexId0];
759 p0[1] = ((
const dataTypeV *)vField_)[vertexId0];
760 p1[0] = ((
const dataTypeU *)uField_)[vertexId1];
761 p1[1] = ((
const dataTypeV *)vField_)[vertexId1];
762 p[0] = (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.first;
763 p[1] = (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_.second;
765 p0.data(), p1.data(), p.data(), baryCentrics, 2);
769 triangulation->getVertexPoint(vertexId0, pA[0], pA[1], pA[2]);
770 triangulation->getVertexPoint(vertexId1, pB[0], pB[1], pB[2]);
773 for(
int j = 0; j < 3; j++) {
777 c0 = pointSet_[3 * vertexId0 + j];
778 c1 = pointSet_[3 * vertexId1 + j];
784 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
785 = baryCentrics[0] * c0 + baryCentrics[1] * c1;
788 if(vertexId0 < vertexId1)
789 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
790 = std::pair<SimplexId, SimplexId>(vertexId0, vertexId1);
792 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
793 = std::pair<SimplexId, SimplexId>(vertexId1, vertexId0);
800template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
816 const triangulationType *
const triangulation)
const {
818 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
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_
826 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
827 = std::pair<SimplexId, SimplexId>(-1, -1);
831 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
832 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 3);
834 for(
int i = 0; i < 3; i++) {
836 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
837 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 1;
838 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
843 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
846 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
849 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
854 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
857 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
860 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
865 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
868 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
871 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i]
879 std::array<std::array<double, 3>, 3> basePoints{};
880 std::array<std::pair<double, double>, 3> basePointProjections{};
881 std::array<double, 3> basePointParameterization{};
882 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
884 computeBaseTriangle<dataTypeU, dataTypeV>(
885 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
886 u2, v2, basePoints, basePointProjections, basePointParameterization,
887 baseEdges, triangulation);
892 if((t0 >= 0) && (t0 <= 1)) {
895 if((t1 >= 0) && (t1 <= 1)) {
898 if((t2 >= 0) && (t2 <= 1)) {
903 for(
int i = 0; i < 5; i++) {
905 SimplexId vertexId0 = -1, vertexId1 = -1;
910 for(
int j = 0; j < 3; j++) {
911 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
912 = basePoints[pivotVertexId][j];
915 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
916 = basePointParameterization[pivotVertexId];
917 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
918 = basePointProjections[pivotVertexId];
919 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
920 = baseEdges[pivotVertexId];
929 vertexId0 = pivotVertexId;
930 vertexId1 = (pivotVertexId + 2) % 3;
931 if(basePointParameterization[(pivotVertexId + 2) % 3] > 1) {
942 vertexId0 = pivotVertexId;
943 vertexId1 = (pivotVertexId + 1) % 3;
944 if(basePointParameterization[(pivotVertexId + 1) % 3] > 1) {
952 vertexId0 = (pivotVertexId + 2) % 3;
953 vertexId1 = (pivotVertexId + 1) % 3;
954 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
962 vertexId0 = (pivotVertexId + 2) % 3;
963 vertexId1 = (pivotVertexId + 1) % 3;
964 if(basePointParameterization[(pivotVertexId + 2) % 3] < 0) {
972 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
974 interpolateBasePoints(
975 basePoints[vertexId0], basePointProjections[vertexId0],
976 basePointParameterization[vertexId0], basePoints[vertexId1],
977 basePointProjections[vertexId1], basePointParameterization[vertexId1],
978 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
989template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1005 const triangulationType *
const triangulation)
const {
1007 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1010 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1011 for(
int i = 0; i < 4; i++) {
1012 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ =
true;
1013 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1015 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1016 = std::pair<SimplexId, SimplexId>(-1, -1);
1020 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1021 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1023 for(
int i = 0; i < 2; i++) {
1025 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1026 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 2;
1027 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1031 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1033 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1035 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1038 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1040 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1042 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1048 std::array<std::array<double, 3>, 3> basePoints{};
1049 std::array<std::pair<double, double>, 3> basePointProjections{};
1050 std::array<double, 3> basePointParameterization{};
1051 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1053 computeBaseTriangle<dataTypeU, dataTypeV>(
1054 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1055 u2, v2, basePoints, basePointProjections, basePointParameterization,
1056 baseEdges, triangulation);
1059 bool isPivotPositive =
false;
1062 if(((t0 < 0) && ((t1 < 0) || (t2 < 0)))
1063 || ((t1 < 0) && ((t0 < 0) || (t2 < 0)))
1064 || ((t2 < 0) && ((t1 < 0) || (t0 < 0)))) {
1065 isPivotPositive =
true;
1067 if(isPivotPositive) {
1074 if(isPivotPositive) {
1081 if(isPivotPositive) {
1090 for(
int i = 0; i < 4; i++) {
1092 SimplexId vertexId0 = -1, vertexId1 = -1;
1101 vertexId0 = pivotVertexId;
1102 vertexId1 = (pivotVertexId + 2) % 3;
1113 vertexId0 = pivotVertexId;
1114 vertexId1 = (pivotVertexId + 1) % 3;
1125 vertexId0 = pivotVertexId;
1126 vertexId1 = (pivotVertexId + 2) % 3;
1137 vertexId0 = pivotVertexId;
1138 vertexId1 = (pivotVertexId + 1) % 3;
1146 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1148 interpolateBasePoints(
1149 basePoints[vertexId0], basePointProjections[vertexId0],
1150 basePointParameterization[vertexId0], basePoints[vertexId1],
1151 basePointProjections[vertexId1], basePointParameterization[vertexId1], t,
1152 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1162template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1178 const triangulationType *
const triangulation)
const {
1180 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1183 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 3);
1184 for(
int i = 0; i < 3; i++) {
1185 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ =
true;
1186 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1188 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1189 = std::pair<SimplexId, SimplexId>(-1, -1);
1193 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1194 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 1);
1196 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].tetId_ = tetId;
1197 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].caseId_ = 3;
1198 (*polygonEdgeTriangleLists_[polygonEdgeId]).back().polygonEdgeId_
1201 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[0]
1203 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[1]
1205 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId].vertexIds_[2]
1209 std::array<std::array<double, 3>, 3> basePoints{};
1210 std::array<std::pair<double, double>, 3> basePointProjections{};
1211 std::array<double, 3> basePointParameterization{};
1212 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1214 computeBaseTriangle<dataTypeU, dataTypeV>(
1215 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1216 u2, v2, basePoints, basePointProjections, basePointParameterization,
1217 baseEdges, triangulation);
1220 bool isPivotPositive =
false;
1223 if((t0 <= 1) && (t0 >= 0)) {
1227 isPivotPositive =
true;
1229 if((t1 <= 1) && (t1 >= 0)) {
1233 isPivotPositive =
true;
1235 if((t2 <= 1) && (t2 >= 0)) {
1239 isPivotPositive =
true;
1243 for(
int i = 0; i < 3; i++) {
1245 SimplexId vertexId0 = -1, vertexId1 = -1;
1250 for(
int j = 0; j < 3; j++) {
1251 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].p_[j]
1252 = basePoints[pivotVertexId][j];
1255 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].t_
1256 = basePointParameterization[pivotVertexId];
1257 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].uv_
1258 = basePointProjections[pivotVertexId];
1259 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId].meshEdge_
1260 = baseEdges[pivotVertexId];
1266 vertexId0 = pivotVertexId;
1267 vertexId1 = (pivotVertexId + 1) % 3;
1277 vertexId0 = pivotVertexId;
1278 vertexId1 = (pivotVertexId + 2) % 3;
1285 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1287 interpolateBasePoints(
1288 basePoints[vertexId0], basePointProjections[vertexId0],
1289 basePointParameterization[vertexId0], basePoints[vertexId1],
1290 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1291 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1302template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1318 const triangulationType *
const triangulation)
const {
1320 SimplexId vertexId = (*polygonEdgeVertexLists_[polygonEdgeId]).size();
1323 (*polygonEdgeVertexLists_[polygonEdgeId]).resize(vertexId + 4);
1324 for(
int i = 0; i < 4; i++) {
1325 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isBasePoint_ =
true;
1326 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].isIntersectionPoint_
1328 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1329 = std::pair<SimplexId, SimplexId>(-1, -1);
1333 SimplexId triangleId = (*polygonEdgeTriangleLists_[polygonEdgeId]).size();
1334 (*polygonEdgeTriangleLists_[polygonEdgeId]).resize(triangleId + 2);
1336 for(
int i = 0; i < 2; i++) {
1338 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].tetId_ = tetId;
1339 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].caseId_ = 4;
1340 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].polygonEdgeId_
1344 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1346 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1348 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1351 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[0]
1353 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[1]
1355 (*polygonEdgeTriangleLists_[polygonEdgeId])[triangleId + i].vertexIds_[2]
1361 std::array<std::array<double, 3>, 3> basePoints{};
1362 std::array<std::pair<double, double>, 3> basePointProjections{};
1363 std::array<double, 3> basePointParameterization{};
1364 std::array<std::pair<SimplexId, SimplexId>, 3> baseEdges{};
1366 computeBaseTriangle<dataTypeU, dataTypeV>(
1367 tetId, localEdgeId0, t0, u0, v0, localEdgeId1, t1, u1, v1, localEdgeId2, t2,
1368 u2, v2, basePoints, basePointProjections, basePointParameterization,
1369 baseEdges, triangulation);
1372 bool isPivotPositive =
false;
1377 isPivotPositive =
true;
1380 isPivotPositive =
false;
1385 isPivotPositive =
true;
1388 isPivotPositive =
false;
1392 isPivotPositive =
true;
1395 isPivotPositive =
false;
1399 for(
int i = 0; i < 4; i++) {
1401 SimplexId vertexId0 = -1, vertexId1 = -1;
1409 vertexId0 = pivotVertexId;
1410 vertexId1 = (pivotVertexId + 2) % 3;
1419 vertexId0 = pivotVertexId;
1420 vertexId1 = (pivotVertexId + 1) % 3;
1427 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_ = t;
1429 interpolateBasePoints(
1430 basePoints[vertexId0], basePointProjections[vertexId0],
1431 basePointParameterization[vertexId0], basePoints[vertexId1],
1432 basePointProjections[vertexId1], basePointParameterization[vertexId1],
1433 t, (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i]);
1441 for(
int j = 0; j < 3; j++) {
1442 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1443 = basePoints[(pivotVertexId + 2) % 3][j];
1446 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1447 = basePointParameterization[(pivotVertexId + 2) % 3];
1448 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1449 = basePointProjections[(pivotVertexId + 2) % 3];
1450 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1451 = baseEdges[(pivotVertexId + 2) % 3];
1454 for(
int j = 0; j < 3; j++) {
1455 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].p_[j]
1456 = basePoints[(pivotVertexId + 1) % 3][j];
1458 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].t_
1459 = basePointParameterization[(pivotVertexId + 1) % 3];
1460 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].uv_
1461 = basePointProjections[(pivotVertexId + 1) % 3];
1462 (*polygonEdgeVertexLists_[polygonEdgeId])[vertexId + i].meshEdge_
1463 = baseEdges[(pivotVertexId + 1) % 3];
1472template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1474 const std::pair<double, double> &rangePoint0,
1475 const std::pair<double, double> &rangePoint1,
1476 const std::vector<SimplexId> &seedTetList,
1477 const triangulationType *
const triangulation,
1480#ifndef TTK_ENABLE_KAMIKAZE
1487 if(!polygonEdgeNumber_)
1489 if(!globalVertexList_)
1493 std::vector<bool> visitedTets(triangulation->getNumberOfCells(),
false);
1494 std::queue<SimplexId> tetQueue;
1498 tetQueue.push(seedTetList[i]);
1508 if(!visitedTets[tetId]) {
1510 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1511 tetId, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1513 if(createdVertices) {
1516 = triangulation->getCellNeighborNumber(tetId);
1518 for(
SimplexId i = 0; i < tetNeighborNumber; i++) {
1521 triangulation->getCellNeighbor(tetId, i, neighborId);
1523 if(!visitedTets[neighborId])
1524 tetQueue.push(neighborId);
1528 visitedTets[tetId] =
true;
1531 }
while(tetQueue.size());
1536template <
class dataTypeU,
class dataTypeV>
1539 std::pair<std::pair<double, double>, std::pair<double, double>>> &edgeList,
1540 const std::vector<SimplexId> &seedTetList,
1541 const std::vector<SimplexId> *edgeIdList)
const {
1543#ifndef TTK_ENABLE_KAMIKAZE
1556 if(!polygonEdgeNumber_)
1558 if(!globalVertexList_)
1562 std::vector<bool> visitedTets(tetNumber_,
false);
1563 std::queue<SimplexId> tetQueue;
1567 tetQueue.push(seedTetList[i]);
1577 if(!visitedTets[tetId]) {
1579 std::vector<std::vector<SimplexId>> threadedTetQueue(edgeList.size());
1580#ifdef TTK_ENABLE_OPENMP
1581#pragma omp parallel for num_threads(threadNumber_)
1588 polygonEdgeId = (*edgeIdList)[i];
1591 createdVertices = processTetrahedron<dataTypeU, dataTypeV>(
1592 tetId, edgeList[i].first, edgeList[i].second, polygonEdgeId);
1594 if(createdVertices) {
1598 if(!visitedTets[(*tetNeighbors_)[tetId][j]]) {
1599 threadedTetQueue[i].push_back((*tetNeighbors_)[tetId][j]);
1605 visitedTets[tetId] =
true;
1609 tetQueue.push(threadedTetQueue[i][j]);
1614 }
while(tetQueue.size());
1619template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1621 const std::pair<double, double> &rangePoint0,
1622 const std::pair<double, double> &rangePoint1,
1623 const triangulationType *
const triangulation,
1626#ifndef TTK_ENABLE_KAMIKAZE
1627 if((!tetNumber_) && (!triangulation))
1629 if((!tetList_) && (!triangulation))
1635 if((!pointSet_) && (!triangulation))
1637 if(!polygonEdgeNumber_)
1639 if(!globalVertexList_)
1646 tetNumber = triangulation->getNumberOfCells();
1649#ifdef TTK_ENABLE_OPENMP
1650#pragma omp parallel for num_threads(threadNumber_)
1652 for(
SimplexId i = 0; i < tetNumber; i++) {
1654 processTetrahedron<dataTypeU, dataTypeV>(
1655 i, rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1661template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1663 const triangulationType *
const triangulation) {
1665#ifndef TTK_ENABLE_KAMIKAZE
1666 if((!tetNumber_) && (!triangulation))
1668 if((!tetList_) && (!triangulation))
1674 if((!pointSet_) && (!triangulation))
1678 if(polygonEdgeNumber_ != (
SimplexId)polygon_->size())
1680 if(!globalVertexList_)
1686#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1687 if(!octree_.empty()) {
1689#ifdef TTK_ENABLE_OPENMP
1690#pragma omp parallel for num_threads(threadNumber_)
1692 for(
SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1694 computeSurfaceWithOctree<dataTypeU, dataTypeV>(
1695 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1699#ifdef TTK_ENABLE_OPENMP
1700#pragma omp parallel for num_threads(threadNumber_)
1702 for(
SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1703 computeSurface<dataTypeU, dataTypeV>(
1704 (*polygon_)[i].first, (*polygon_)[i].second, triangulation, i);
1709#ifdef TTK_ENABLE_OPENMP
1710#pragma omp parallel for num_threads(threadNumber_)
1712 for(
SimplexId i = 0; i < polygonEdgeNumber_; i++) {
1714 computeSurface<dataTypeU, dataTypeV>(
1715 (*polygon_)[i].first, (*polygon_)[i].second, i);
1719 finalize<dataTypeU, dataTypeV>(pointSnapping_,
false,
false,
false);
1726#ifdef TTK_ENABLE_FIBER_SURFACE_WITH_RANGE_OCTREE
1727template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1728inline int ttk::FiberSurface::computeSurfaceWithOctree(
1729 const std::pair<double, double> &rangePoint0,
1730 const std::pair<double, double> &rangePoint1,
1731 const triangulationType *
const triangulation,
1734#ifndef TTK_ENABLE_KAMIKAZE
1735 if((!tetNumber_) && (!triangulation))
1737 if((!tetList_) && (!triangulation))
1743 if((!pointSet_) && (!triangulation))
1745 if(!polygonEdgeNumber_)
1747 if(!globalVertexList_)
1751 std::vector<SimplexId> tetList;
1752 octree_.rangeSegmentQuery(rangePoint0, rangePoint1, tetList);
1754#ifdef TTK_ENABLE_OPENMP
1755#pragma omp parallel for num_threads(threadNumber_)
1758 processTetrahedron<dataTypeU, dataTypeV>(
1759 tetList[i], rangePoint0, rangePoint1, triangulation, polygonEdgeId);
1766template <
class dataTypeU,
class dataTypeV>
1768 const bool &removeSmallEdges,
1769 const bool &edgeFlips,
1770 const bool &intersectionRemesh) {
1775 fiberSurfaceVertexNumber += (*polygonEdgeVertexLists_[i]).size();
1778 (*globalVertexList_).resize(fiberSurfaceVertexNumber);
1779 fiberSurfaceVertexNumber = 0;
1783 (*polygonEdgeVertexLists_[i])[j].polygonEdgeId_ = i;
1784 (*polygonEdgeVertexLists_[i])[j].localId_ = j;
1785 (*polygonEdgeVertexLists_[i])[j].globalId_ = fiberSurfaceVertexNumber;
1786 (*globalVertexList_)[fiberSurfaceVertexNumber]
1787 = (*polygonEdgeVertexLists_[i])[j];
1788 fiberSurfaceVertexNumber++;
1794 for(
int k = 0; k < 3; k++) {
1795 (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k]
1796 = (*polygonEdgeVertexLists_[i])[(*polygonEdgeTriangleLists_[i])[j]
1808 if(intersectionRemesh) {
1809 remeshIntersections<dataTypeU, dataTypeV>();
1812 if((mergeDuplicatedVertices) || (removeSmallEdges)) {
1814 mergeVertices(pointSnappingThreshold_);
1820 if(removeSmallEdges)
1821 mergeEdges(edgeCollapseThreshold_);
1825 polygonEdgeVertexLists_[i]->clear();
1831template <
class dataTypeU,
class dataTypeV,
typename triangulationType>
1834 const std::pair<double, double> &rangePoint0,
1835 const std::pair<double, double> &rangePoint1,
1836 const triangulationType *
const triangulation,
1839 double rangeEdge[2];
1840 rangeEdge[0] = rangePoint0.first - rangePoint1.first;
1841 rangeEdge[1] = rangePoint0.second - rangePoint1.second;
1843 double rangeNormal[2];
1844 rangeNormal[0] = -rangeEdge[1];
1845 rangeNormal[1] = rangeEdge[0];
1854 for(
int i = 0; i < 4; i++) {
1857 if(!triangulation) {
1858 vertexId = tetList_[5 * tetId + 1 + i];
1860 triangulation->getCellVertex(tetId, i, vertexId);
1863 double projectedVertex[2];
1864 projectedVertex[0] = ((
const dataTypeU *)uField_)[vertexId];
1865 projectedVertex[1] = ((
const dataTypeV *)vField_)[vertexId];
1867 double vertexRangeEdge[2];
1868 vertexRangeEdge[0] = projectedVertex[0] - rangePoint0.first;
1869 vertexRangeEdge[1] = projectedVertex[1] - rangePoint0.second;
1871 d[i] = vertexRangeEdge[0] * rangeNormal[0]
1872 + vertexRangeEdge[1] * rangeNormal[1];
1874 if(fabs(d[i]) < prec_dbl)
1883 equalVertexLocalId = i;
1888 if(!((upperNumber == 0) || (lowerNumber == 0))) {
1891 std::array<SimplexId, 2> triangleEdgeNumbers{};
1892 std::array<std::array<SimplexId, 3>, 2> triangleEdges{
1893 {{-1, -1, -1}, {-1, -1, -1}}};
1903 for(
int i = 0; i < 4; i++) {
1918 for(
SimplexId j = jStart; j != jEnd; j += jStep) {
1920 if(((d[i] > 0) && (d[j] < 0)) || ((d[i] < 0) && (d[j] > 0))) {
1923 if(triangleEdgeNumbers[0] == 3) {
1924 triangleEdges[1][triangleEdgeNumbers[1]] = edgeCounter;
1925 triangleEdgeNumbers[1]++;
1927 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1928 triangleEdgeNumbers[0]++;
1932 if((d[i] == 0) && (d[j] == 0)) {
1938 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1939 triangleEdgeNumbers[0]++;
1941 triangleEdges[0][triangleEdgeNumbers[0]] = edgeCounter;
1942 triangleEdgeNumbers[0]++;
1950 if(triangleEdges[1][0] != -1) {
1951 if(triangleEdges[1][1] == -1) {
1960 switch(triangleEdges[1][0]) {
1981 if(triangleEdges[0][i] != forbiddenEdge) {
1982 if(triangleEdges[1][1] != -1) {
1983 triangleEdges[1][2] = triangleEdges[0][i];
1986 triangleEdges[1][1] = triangleEdges[0][i];
1994 for(
int i = 0; i < 2; i++) {
1995 if(triangleEdges[i][0] != -1) {
1997 if(triangleEdges[i][2] == -1) {
1999 switch(equalVertexLocalId) {
2002 triangleEdges[i][2] = 0;
2006 triangleEdges[i][2] = 5;
2020 std::pair<double, double> uv0, uv1;
2021 std::array<std::pair<double, double>, 3> uv{};
2022 std::array<double, 3> t{};
2026 for(
int i = 0; i < 2; i++) {
2027 if(triangleEdges[i][0] != -1) {
2038 triangulation->getCellVertex(
2039 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j]], vertexId0);
2040 triangulation->getCellVertex(
2041 tetId, edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1],
2045 = tetList_[5 * tetId + 1
2046 + edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2048 = tetList_[5 * tetId + 1
2049 + edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2052 if((j < (
SimplexId)triangleEdges[i].size() - 1)
2053 && (triangleEdges[i][j] == triangleEdges[i][j + 1])) {
2056 uv[j].first = ((
const dataTypeU *)uField_)[vertexId0];
2057 uv[j].second = ((
const dataTypeV *)vField_)[vertexId0];
2059 uv[j + 1].first = ((
const dataTypeU *)uField_)[vertexId1];
2060 uv[j + 1].second = ((
const dataTypeV *)vField_)[vertexId1];
2064 && (triangleEdges[i][j] != triangleEdges[i][j - 1]))) {
2067 d0 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j]]];
2068 uv0.first = ((
const dataTypeU *)uField_)[vertexId0];
2069 uv0.second = ((
const dataTypeV *)vField_)[vertexId0];
2071 d1 = d[edgeImplicitEncoding_[2 * triangleEdges[i][j] + 1]];
2072 uv1.first = ((
const dataTypeU *)uField_)[vertexId1];
2073 uv1.second = ((
const dataTypeV *)vField_)[vertexId1];
2076 = uv0.first + (d0 / (d0 - d1)) * (uv1.first - uv0.first);
2078 = uv0.second + (d0 / (d0 - d1)) * (uv1.second - uv0.second);
2083 if(fabs(rangePoint1.first - rangePoint0.first)
2084 > fabs(rangePoint1.second - rangePoint0.second)) {
2085 t[j] = (uv[j].first - rangePoint0.first)
2086 / (rangePoint1.first - rangePoint0.first);
2088 t[j] = (uv[j].second - rangePoint0.second)
2089 / (rangePoint1.second - rangePoint0.second);
2092 if((t[j] <= 1) && (t[j] >= 0))
2095 lowerVertexNumber++;
2097 upperVertexNumber++;
2104 if(greyVertexNumber == 3) {
2105 createdVertices += computeCase0<dataTypeU, dataTypeV>(
2106 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2107 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2108 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2110 }
else if(lowerVertexNumber == 3 || upperVertexNumber == 3) {
2112 }
else if((lowerVertexNumber == 1) && (upperVertexNumber == 1)
2113 && (greyVertexNumber == 1)) {
2114 createdVertices += computeCase1<dataTypeU, dataTypeV>(
2115 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2116 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2117 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2119 }
else if(((lowerVertexNumber == 2) && (upperVertexNumber == 1))
2120 || ((lowerVertexNumber == 1) && (upperVertexNumber == 2))) {
2121 createdVertices += computeCase2<dataTypeU, dataTypeV>(
2122 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2123 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2124 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2126 }
else if((greyVertexNumber == 1)
2127 && ((lowerVertexNumber == 2) || (upperVertexNumber == 2))) {
2128 createdVertices += computeCase3<dataTypeU, dataTypeV>(
2129 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2130 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2131 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2133 }
else if(((greyVertexNumber == 2))
2134 && ((lowerVertexNumber == 1) || (upperVertexNumber == 1))) {
2135 createdVertices += computeCase4<dataTypeU, dataTypeV>(
2136 polygonEdgeId, tetId, triangleEdges[i][0], t[0], uv[0].first,
2137 uv[0].second, triangleEdges[i][1], t[1], uv[1].first, uv[1].second,
2138 triangleEdges[i][2], t[2], uv[2].first, uv[2].second,
2148 if((triangleEdges[1][0] != -1) && (createdVertices > 3)) {
2150 std::vector<SimplexId> createdVertexList(createdVertices);
2152 createdVertexList[i]
2153 = polygonEdgeVertexLists_[polygonEdgeId]->size() - 1 - i;
2156 std::vector<bool> snappedVertices(createdVertices,
false);
2157 std::vector<SimplexId> colinearVertices;
2158 colinearVertices.reserve(createdVertices);
2160 for(
SimplexId i = 0; i < createdVertices; i++) {
2161 colinearVertices.clear();
2163 if(!snappedVertices[i]) {
2164 colinearVertices.push_back(i);
2165 for(
SimplexId j = 0; j < createdVertices; j++) {
2166 if((i != j) && (!snappedVertices[j])) {
2170 if((*polygonEdgeVertexLists_[polygonEdgeId])[createdVertexList[i]]
2172 == (*polygonEdgeVertexLists_[polygonEdgeId])
2173 [createdVertexList[j]]
2175 colinearVertices.push_back(j);
2181 if((colinearVertices.size() == 4) || (colinearVertices.size() == 3)) {
2186 std::pair<SimplexId, SimplexId> minPair;
2187 double minDistance = -1;
2193 (*polygonEdgeVertexLists_[polygonEdgeId])
2194 [createdVertexList[colinearVertices[j]]]
2196 (*polygonEdgeVertexLists_[polygonEdgeId])
2197 [createdVertexList[colinearVertices[k]]]
2213 if((minDistance < 0) || (distance < minDistance)) {
2214 minDistance = distance;
2222 if((minDistance != -1) && (minDistance < prec_dbl)) {
2227 if((j != minPair.first) && (j != minPair.second)) {
2230 for(
int k = 0; k < 3; k++) {
2231 (*polygonEdgeVertexLists_[polygonEdgeId])
2232 [createdVertexList[colinearVertices[minPair.first]]]
2234 = (*polygonEdgeVertexLists_[polygonEdgeId])
2235 [createdVertexList[colinearVertices[j]]]
2238 (*polygonEdgeVertexLists_[polygonEdgeId])
2239 [createdVertexList[colinearVertices[minPair.first]]]
2241 = (*polygonEdgeVertexLists_[polygonEdgeId])
2242 [createdVertexList[colinearVertices[j]]]
2245 (*polygonEdgeVertexLists_[polygonEdgeId])
2246 [createdVertexList[colinearVertices[minPair.first]]]
2248 = (*polygonEdgeVertexLists_[polygonEdgeId])
2249 [createdVertexList[colinearVertices[j]]]
2252 (*polygonEdgeVertexLists_[polygonEdgeId])
2253 [createdVertexList[colinearVertices[minPair.first]]]
2255 = (*polygonEdgeVertexLists_[polygonEdgeId])
2256 [createdVertexList[colinearVertices[j]]]
2259 (*polygonEdgeVertexLists_[polygonEdgeId])
2260 [createdVertexList[colinearVertices[minPair.first]]]
2261 .isIntersectionPoint_
2262 = (*polygonEdgeVertexLists_[polygonEdgeId])
2263 [createdVertexList[colinearVertices[j]]]
2264 .isIntersectionPoint_;
2265 if((*polygonEdgeVertexLists_[polygonEdgeId])
2266 [createdVertexList[colinearVertices[j]]]
2269 (*polygonEdgeVertexLists_[polygonEdgeId])
2270 [createdVertexList[colinearVertices[minPair.first]]]
2272 = (*polygonEdgeVertexLists_[polygonEdgeId])
2273 [createdVertexList[colinearVertices[j]]]
2278 for(
int k = 0; k < 3; k++) {
2279 (*polygonEdgeVertexLists_[polygonEdgeId])
2280 [createdVertexList[colinearVertices[minPair.second]]]
2282 = (*polygonEdgeVertexLists_[polygonEdgeId])
2283 [createdVertexList[colinearVertices[j]]]
2286 (*polygonEdgeVertexLists_[polygonEdgeId])
2287 [createdVertexList[colinearVertices[minPair.second]]]
2289 = (*polygonEdgeVertexLists_[polygonEdgeId])
2290 [createdVertexList[colinearVertices[j]]]
2293 (*polygonEdgeVertexLists_[polygonEdgeId])
2294 [createdVertexList[colinearVertices[minPair.second]]]
2296 = (*polygonEdgeVertexLists_[polygonEdgeId])
2297 [createdVertexList[colinearVertices[j]]]
2300 (*polygonEdgeVertexLists_[polygonEdgeId])
2301 [createdVertexList[colinearVertices[minPair.second]]]
2303 = (*polygonEdgeVertexLists_[polygonEdgeId])
2304 [createdVertexList[colinearVertices[j]]]
2307 (*polygonEdgeVertexLists_[polygonEdgeId])
2308 [createdVertexList[colinearVertices[minPair.second]]]
2309 .isIntersectionPoint_
2310 = (*polygonEdgeVertexLists_[polygonEdgeId])
2311 [createdVertexList[colinearVertices[j]]]
2312 .isIntersectionPoint_;
2313 if((*polygonEdgeVertexLists_[polygonEdgeId])
2314 [createdVertexList[colinearVertices[j]]]
2317 (*polygonEdgeVertexLists_[polygonEdgeId])
2318 [createdVertexList[colinearVertices[minPair.second]]]
2320 = (*polygonEdgeVertexLists_[polygonEdgeId])
2321 [createdVertexList[colinearVertices[j]]]
2325 snappedVertices[colinearVertices[minPair.first]] =
true;
2326 snappedVertices[colinearVertices[minPair.second]] =
true;
2337 return createdVertices;
2343template <
class dataTypeU,
class dataTypeV>
2346#ifndef TTK_ENABLE_KAMIKAZE
2366 std::vector<std::vector<IntersectionTriangle>> tetIntersections(tetNumber_);
2367 std::vector<std::vector<Vertex>> tetNewVertices(tetNumber_);
2374 SimplexId tetId = (*polygonEdgeTriangleLists_[i])[j].tetId_;
2376 tetIntersections[tetId].emplace_back();
2378 tetIntersections[tetId].back().caseId_
2379 = (*polygonEdgeTriangleLists_[i])[j].caseId_;
2380 tetIntersections[tetId].back().polygonEdgeId_ = i;
2381 tetIntersections[tetId].back().triangleId_ = j;
2382 for(
int k = 0; k < 3; k++) {
2383 tetIntersections[tetId].back().vertexIds_[k]
2384 = (*polygonEdgeTriangleLists_[i])[j].vertexIds_[k];
2385 tetIntersections[tetId].back().uv_[k]
2386 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2388 tetIntersections[tetId].back().t_[k]
2389 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2391 for(
int l = 0; l < 3; l++) {
2392 tetIntersections[tetId].back().p_[k][l]
2393 = (*globalVertexList_)[tetIntersections[tetId].back().vertexIds_[k]]
2397 tetIntersections[tetId].back().intersection_.first = -DBL_MAX;
2398 tetIntersections[tetId].back().intersection_.second = -DBL_MAX;
2402 std::vector<SimplexId> tetList;
2404 if(tetIntersections[i].size() > 1)
2405 tetList.push_back(i);
2408#ifdef TTK_ENABLE_OPENMP
2409#pragma omp parallel for num_threads(threadNumber_)
2422 this->printWrn(
"Preventing an infinite loop!");
2423 this->printWrn(
"More than 1000 re-meshed triangles in tet #"
2424 + std::to_string(tetId));
2425 this->printWrn(
"Extra-thin triangles keep on intersecting?!");
2432 = (
SimplexId)tetIntersections[tetId].size();
2434 for(
SimplexId k = 0; k < originalTriangleNumber; k++) {
2436 SimplexId polygonEdgeId0 = tetIntersections[tetId][j].polygonEdgeId_;
2437 SimplexId polygonEdgeId1 = tetIntersections[tetId][k].polygonEdgeId_;
2439 if((j != k) && (polygonEdgeId0 != polygonEdgeId1)) {
2444 std::pair<double, double> edge0point0, edge0point1;
2446 getTriangleRangeExtremities(
2447 tetId, j, tetIntersections, edge0point0, edge0point1);
2450 std::pair<double, double> edge1point0, edge1point1;
2452 getTriangleRangeExtremities(
2453 tetId, k, tetIntersections, edge1point0, edge1point1);
2456 std::pair<double, double> intersection;
2458 edge0point0.first, edge0point0.second, edge0point1.first,
2459 edge0point1.second, edge1point0.first, edge1point0.second,
2460 edge1point1.first, edge1point1.second, intersection.first,
2461 intersection.second);
2463 if((hasIntersection)
2475 computeTriangleIntersection(tetId, j, k, polygonEdgeId0,
2476 polygonEdgeId1, intersection,
2477 newVertexNumber, newTriangleNumber,
2478 tetIntersections, tetNewVertices);
2489 SimplexId localId = (*globalVertexList_).size();
2490 tetNewVertices[i][j].localId_ = localId;
2491 (*globalVertexList_).push_back(tetNewVertices[i][j]);
2498 if(((tetIntersections[i][j].intersection_.first != -DBL_MAX)
2499 && (tetIntersections[i][j].intersection_.second != -DBL_MAX))
2500 || (tetIntersections[i][j].triangleId_ < 0)) {
2502 SimplexId triangleId = tetIntersections[i][j].triangleId_;
2504 if(triangleId < 0) {
2507 triangleId = (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2510 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2511 .resize(triangleId + 1);
2512 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2516 (*polygonEdgeTriangleLists_[tetIntersections[i][j].polygonEdgeId_])
2519 = tetIntersections[i][j].caseId_;
2522 for(
int k = 0; k < 3; k++) {
2524 SimplexId vertexId = tetIntersections[i][j].vertexIds_[k];
2528 vertexId = tetNewVertices[i][-(vertexId + 1)].localId_;
2530 (*polygonEdgeTriangleLists_[tetIntersections[i][j]
2531 .polygonEdgeId_])[triangleId]
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual int preconditionCellNeighbors()
Minimalist debugging class.
virtual int setDebugLevel(const int &debugLevel)
TTK processing package that computes fiber surfaces.
std::vector< std::vector< Vertex > * > polygonEdgeVertexLists_
int setGlobalVertexList(std::vector< Vertex > *globalList)
int interpolateBasePoints(const std::array< double, 3 > &p0, const std::pair< double, double > &uv0, const double &t0, const std::array< double, 3 > &p1, const std::pair< double, double > &uv1, const double &t1, const double &t, Vertex &v) const
std::array< SimplexId, 12 > edgeImplicitEncoding_
double pointSnappingThreshold_
double edgeCollapseThreshold_
std::vector< std::vector< Triangle > * > polygonEdgeTriangleLists_
const std::vector< std::vector< SimplexId > > * tetNeighbors_
int getNumberOfCommonVertices(const SimplexId &tetId, const SimplexId &triangleId0, const SimplexId &triangleId1, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections) const
int mergeEdges(const double &distanceThreshold) const
int setInputField(const void *uField, const void *vField)
int setTetList(const SimplexId *tetList)
int computeCase3(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
int getTriangleRangeExtremities(const SimplexId &tetId, const SimplexId &triangleId, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::pair< double, double > &extremity0, std::pair< double, double > &extremity1) const
int finalize(const bool &mergeDuplicatedVertices=false, const bool &removeSmallEdges=false, const bool &edgeFlips=false, const bool &intersectionRemesh=false)
const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > * polygon_
int computeCase4(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
int computeCase1(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
int computeSurface(const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
int remeshIntersections() const
int computeCase0(const SimplexId &polygonEdgeId, const SimplexId &tetId, const SimplexId &localEdgeId0, const double &t0, const double &u0, const double &v0, const SimplexId &localEdgeId1, const double &t1, const double &u1, const double &v1, const SimplexId &localEdgeId2, const double &t2, const double &u2, const double &v2, const triangulationType *const triangulation) const
int snapVertexBarycentrics() const
bool isEdgeAngleCollapsible(const SimplexId &source, const SimplexId &destination, const SimplexId &pivotVertexId, const std::vector< std::pair< SimplexId, SimplexId > > &starNeighbors) const
int processTetrahedron(const SimplexId &tetId, const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
bool isIntersectionTriangleColinear(const SimplexId &tetId, const SimplexId &triangleId, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, const std::vector< std::vector< Vertex > > &tetNewVertices, const SimplexId &vertexId0, const SimplexId &vertexId1, const SimplexId &vertexId2) const
bool hasDuplicatedVertices(const double *p0, const double *p1, const double *p2) const
int setPolygonEdgeNumber(const SimplexId &polygonEdgeNumber)
int mergeVertices(const double &distanceThreshold) const
int setTetNumber(const SimplexId &tetNumber)
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
std::vector< Vertex > * globalVertexList_
int computeTriangleFiber(const SimplexId &tetId, const SimplexId &triangleId, const std::pair< double, double > &intersection, const std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::array< double, 3 > &pA, std::array< double, 3 > &pB, SimplexId &pivotVertexId, bool &edgeFiber) const
int setPointSet(const float *pointSet)
int setPointNumber(const SimplexId &number)
int setPointMerging(const bool &onOff)
int setPointMergingThreshold(const double &threshold)
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
int setPolygon(const std::vector< std::pair< std::pair< double, double >, std::pair< double, double > > > *polygon)
const SimplexId * tetList_
SimplexId polygonEdgeNumber_
int snapToBasePoint(const std::vector< std::vector< double > > &basePoints, const std::vector< std::pair< double, double > > &uv, const std::vector< double > &t, Vertex &v) const
int computeTriangleIntersection(const SimplexId &tetId, const SimplexId &triangleId0, const SimplexId &triangleId1, const SimplexId &polygonEdgeId0, const SimplexId &polygonEdgeId1, const std::pair< double, double > &intersection, SimplexId &newVertexNumber, SimplexId &newTriangleNumber, std::vector< std::vector< IntersectionTriangle > > &tetIntersections, std::vector< std::vector< Vertex > > &tetNewVertices) const
bool isEdgeFlippable(const SimplexId &edgeVertexId0, const SimplexId &edgeVertexId1, const SimplexId &otherVertexId0, const SimplexId &otherVertexId1) const
int setVertexList(const SimplexId &polygonEdgeId, std::vector< Vertex > *vertexList)
int setTetNeighbors(const std::vector< std::vector< SimplexId > > *tetNeighbors)
int setTriangleList(const SimplexId &polygonEdgeId, std::vector< Triangle > *triangleList)
int createNewIntersectionTriangle(const SimplexId &tetId, const SimplexId &triangleId, const SimplexId &vertexId0, const SimplexId &vertexId1, const SimplexId &vertexId2, const std::vector< std::vector< Vertex > > &tetNewVertices, SimplexId &newTriangleNumber, std::vector< std::vector< IntersectionTriangle > > &tetIntersections, const std::pair< double, double > *intersection=nullptr) const
void preconditionTriangulation(AbstractTriangulation *triangulation)
int computeContour(const std::pair< double, double > &rangePoint0, const std::pair< double, double > &rangePoint1, const std::vector< SimplexId > &seedTetList, const triangulationType *const triangulation, const SimplexId &polygonEdgeId=0) const
TTK optional package for octree based range queries in bivariate volumetric data.
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)
int computeBarycentricCoordinates(const T *p0, const T *p1, const T *p, std::array< T, 2 > &baryCentrics, const int &dimension=3)
bool isTriangleColinear(const T *p0, const T *p1, const T *p2, const T *tolerance=nullptr)
T powIntTen(const int n)
Compute the nth power of ten.
T powInt(const T val, const int n)
T distance(const T *p0, const T *p1, const int &dimension=3)
int SimplexId
Identifier type for simplices of any dimension.
std::array< std::array< double, 3 >, 3 > p_
std::pair< double, double > intersection_
std::array< double, 3 > t_
std::array< std::pair< double, double >, 3 > uv_
std::array< SimplexId, 3 > vertexIds_
std::array< SimplexId, 3 > vertexIds_
bool isIntersectionPoint_
std::pair< double, double > uv_
std::array< double, 3 > p_
std::pair< SimplexId, SimplexId > meshEdge_
printMsg(debug::output::GREEN+" "+debug::output::ENDCOLOR+debug::output::GREEN+"▒"+debug::output::ENDCOLOR+debug::output::GREEN+"▒▒▒▒▒▒▒▒▒▒▒▒▒░"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)