68 template <
typename triangulationType>
70 const size_t scalarsMTime,
72 const triangulationType &triangulation) {
92 this->
dg_ = std::move(dg);
97 return std::move(this->
dg_);
100 template <
typename triangulationType>
103 const triangulationType &triangulation) {
107 inline const std::vector<std::vector<SimplexId>> &
126 template <
typename triangulationType>
129 const triangulationType &triangulation,
130 const bool ignoreBoundary,
131 const bool compute2SaddlesChildren =
false);
159 template <
typename triangulationType>
160 std::vector<std::vector<SimplexId>>
162 const triangulationType &triangulation)
const;
179 template <
typename triangulationType,
183 std::vector<std::vector<SimplexId>>
185 const GFS &getFaceStar,
186 const GFSN &getFaceStarNumber,
187 const OB &isOnBoundary,
188 const triangulationType &triangulation)
const;
201 template <
typename triangulationType>
203 std::vector<bool> &pairedMinima,
204 std::vector<bool> &paired1Saddles,
205 const std::vector<SimplexId> &criticalEdges,
206 const std::vector<SimplexId> &critEdgesOrder,
208 const triangulationType &triangulation)
const;
221 template <
typename triangulationType>
223 std::vector<bool> &pairedMaxima,
224 std::vector<bool> &pairedSaddles,
225 const std::vector<SimplexId> &criticalSaddles,
226 const std::vector<SimplexId> &critSaddlesOrder,
227 const std::vector<SimplexId> &critMaxsOrder,
228 const triangulationType &triangulation)
const;
243 template <
typename triangulationType>
245 std::vector<bool> &paired1Saddles,
246 std::vector<bool> &paired2Saddles,
247 const bool exportBoundaries,
248 std::vector<GeneratorType> &boundaries,
249 const std::vector<SimplexId> &critical1Saddles,
250 const std::vector<SimplexId> &critical2Saddles,
251 const std::vector<SimplexId> &crit1SaddlesOrder,
252 const triangulationType &triangulation)
const;
263 template <
typename triangulationType>
265 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
266 std::array<std::vector<SimplexId>, 4> &critCellsOrder,
268 const triangulationType &triangulation,
269 const bool sortEdges)
const;
283 const std::vector<PersistencePair> &pairs,
284 const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
285 const std::vector<bool> &pairedMinima,
286 const std::vector<bool> &paired1Saddles,
287 const std::vector<bool> &paired2Saddles,
288 const std::vector<bool> &pairedMaxima)
const;
312 std::vector<bool> &pairedExtrema,
313 std::vector<bool> &pairedSaddles,
314 std::vector<SimplexId> &reps,
315 std::vector<tripletType> &triplets,
338 template <
typename triangulationType,
typename Container>
341 std::vector<bool> &onBoundary,
342 std::vector<Container> &s2Boundaries,
343 const std::vector<SimplexId> &s2Mapping,
344 const std::vector<SimplexId> &s1Mapping,
345 std::vector<SimplexId> &partners,
346 std::vector<Lock> &s1Locks,
347 std::vector<Lock> &s2Locks,
348 const triangulationType &triangulation)
const;
373 template <
typename triangulationType>
376 const triangulationType &triangulation) {
378 triangulation.getEdgeVertex(
id, 0, this->
vertsOrder_[0]);
379 triangulation.getEdgeVertex(
id, 1, this->
vertsOrder_[1]);
381 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
383 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
391 template <
typename triangulationType>
394 const triangulationType &triangulation) {
396 triangulation.getTriangleVertex(
id, 0, this->
vertsOrder_[0]);
397 triangulation.getTriangleVertex(
id, 1, this->
vertsOrder_[1]);
398 triangulation.getTriangleVertex(
id, 2, this->
vertsOrder_[2]);
400 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
401 this->vertsOrder_[2] = offsets[this->vertsOrder_[2]];
403 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
411 template <
typename triangulationType>
414 const triangulationType &triangulation) {
416 triangulation.getCellVertex(
id, 0, this->
vertsOrder_[0]);
417 triangulation.getCellVertex(
id, 1, this->
vertsOrder_[1]);
418 triangulation.getCellVertex(
id, 2, this->
vertsOrder_[2]);
419 triangulation.getCellVertex(
id, 3, this->
vertsOrder_[3]);
421 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
422 this->vertsOrder_[2] = offsets[this->vertsOrder_[2]];
423 this->vertsOrder_[3] = offsets[this->vertsOrder_[3]];
425 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
429 template <
typename triangulationType>
430 void alloc(
const triangulationType &triangulation) {
433 if(dim > 3 || dim < 1) {
436#ifdef TTK_ENABLE_OPENMP
437#pragma omp parallel master num_threads(threadNumber_)
440#ifdef TTK_ENABLE_OPENMP
443 this->
firstRepMin_.resize(triangulation.getNumberOfVertices());
445#ifdef TTK_ENABLE_OPENMP
448 this->
firstRepMax_.resize(triangulation.getNumberOfCells());
451#ifdef TTK_ENABLE_OPENMP
454 this->
critEdges_.resize(triangulation.getNumberOfEdges());
455#ifdef TTK_ENABLE_OPENMP
459 triangulation.getNumberOfEdges(), -1);
460#ifdef TTK_ENABLE_OPENMP
463 this->
onBoundary_.resize(triangulation.getNumberOfEdges(),
false);
464#ifdef TTK_ENABLE_OPENMP
467 this->
s2Mapping_.resize(triangulation.getNumberOfTriangles(), -1);
468#ifdef TTK_ENABLE_OPENMP
471 this->
s1Mapping_.resize(triangulation.getNumberOfEdges(), -1);
473 for(
int i = 0; i < dim + 1; ++i) {
474#ifdef TTK_ENABLE_OPENMP
480 for(
int i = 1; i < dim + 1; ++i) {
481#ifdef TTK_ENABLE_OPENMP
488 this->
printMsg(
"Memory allocations", 1.0, tm.getElapsedTime(), 1,
503 this->
printMsg(
"Memory cleanup", 1.0, tm.getElapsedTime(), 1,
525template <
typename triangulationType>
526std::vector<std::vector<SimplexId>>
528 const std::vector<SimplexId> &criticalEdges,
529 const triangulationType &triangulation)
const {
533 std::vector<std::vector<SimplexId>> res(criticalEdges.size());
536#ifdef TTK_ENABLE_OPENMP
537#pragma omp parallel for num_threads(threadNumber_)
539 for(
size_t i = 0; i < criticalEdges.size(); ++i) {
542 const auto followVPath = [
this, &mins, &triangulation](
const SimplexId v) {
543 std::vector<Cell> vpath{};
545 const Cell &lastCell = vpath.back();
546 if(lastCell.
dim_ == 0 && this->dg_.isCellCritical(lastCell)) {
547 mins.emplace_back(lastCell.
id_);
553 triangulation.getEdgeVertex(criticalEdges[i], 0, v0);
554 triangulation.getEdgeVertex(criticalEdges[i], 1, v1);
561 this->
printMsg(
"Computed the descending 1-separatrices", 1.0,
568template <
typename triangulationType,
typename GFS,
typename GFSN,
typename OB>
569std::vector<std::vector<SimplexId>>
571 const std::vector<SimplexId> &criticalCells,
572 const GFS &getFaceStar,
573 const GFSN &getFaceStarNumber,
574 const OB &isOnBoundary,
575 const triangulationType &triangulation)
const {
579 const auto dim = this->dg_.getDimensionality();
580 std::vector<std::vector<SimplexId>> res(criticalCells.size());
583#ifdef TTK_ENABLE_OPENMP
584#pragma omp parallel for num_threads(threadNumber_)
586 for(
size_t i = 0; i < criticalCells.size(); ++i) {
587 const auto sid = criticalCells[i];
590 const auto followVPath
591 = [
this, dim, &maxs, &triangulation](
const SimplexId v) {
592 std::vector<Cell> vpath{};
593 this->dg_.getAscendingPath(
Cell{dim, v}, vpath, triangulation);
594 const Cell &lastCell = vpath.back();
595 if(lastCell.
dim_ == dim && this->dg_.isCellCritical(lastCell)) {
596 maxs.emplace_back(lastCell.
id_);
597 }
else if(lastCell.
dim_ == dim - 1) {
598 maxs.emplace_back(-1);
602 const auto starNumber = getFaceStarNumber(sid);
604 for(
SimplexId j = 0; j < starNumber; ++j) {
606 getFaceStar(sid, j, cellId);
610 if(isOnBoundary(sid)) {
612 maxs.emplace_back(-1);
616 this->
printMsg(
"Computed the ascending 1-separatrices", 1.0,
623template <
typename triangulationType>
625 std::vector<PersistencePair> &pairs,
626 std::vector<bool> &pairedMinima,
627 std::vector<bool> &paired1Saddles,
628 const std::vector<SimplexId> &criticalEdges,
629 const std::vector<SimplexId> &critEdgesOrder,
631 const triangulationType &triangulation)
const {
635 auto saddle1ToMinima = getSaddle1ToMinima(criticalEdges, triangulation);
639 auto &firstRep{this->firstRepMin_};
640 std::iota(firstRep.begin(), firstRep.end(), 0);
641 std::vector<tripletType> sadMinTriplets{};
643 for(
size_t i = 0; i < saddle1ToMinima.size(); ++i) {
644 auto &mins = saddle1ToMinima[i];
645 const auto s1 = criticalEdges[i];
647 TTK_PSORT(this->threadNumber_, mins.begin(), mins.end());
648 const auto last = std::unique(mins.begin(), mins.end());
649 mins.erase(last, mins.end());
650 if(mins.size() != 2) {
653 sadMinTriplets.emplace_back(
tripletType{s1, mins[0], mins[1]});
656 tripletsToPersistencePairs(pairs, pairedMinima, paired1Saddles, firstRep,
657 sadMinTriplets, critEdgesOrder.data(), offsets, 0);
659 const auto nMinSadPairs = pairs.size();
662 "Computed " +
std::to_string(nMinSadPairs) +
" min-saddle pairs", 1.0,
663 tm.getElapsedTime(), this->threadNumber_);
665 this->
printMsg(
"min-saddle pairs sequential part", 1.0,
670template <
typename triangulationType>
672 std::vector<PersistencePair> &pairs,
673 std::vector<bool> &pairedMaxima,
674 std::vector<bool> &pairedSaddles,
675 const std::vector<SimplexId> &criticalSaddles,
676 const std::vector<SimplexId> &critSaddlesOrder,
677 const std::vector<SimplexId> &critMaxsOrder,
678 const triangulationType &triangulation)
const {
682 const auto dim = this->dg_.getDimensionality();
686 ? getSaddle2ToMaxima(
689 return triangulation.getTriangleStar(a, i, r);
692 return triangulation.getTriangleStarNumber(a);
695 return triangulation.isTriangleOnBoundary(a);
698 : getSaddle2ToMaxima(
701 return triangulation.getEdgeStar(a, i, r);
704 return triangulation.getEdgeStarNumber(a);
707 return triangulation.isEdgeOnBoundary(a);
713 auto &firstRep{this->firstRepMax_};
714 std::iota(firstRep.begin(), firstRep.end(), 0);
715 std::vector<tripletType> sadMaxTriplets{};
717 for(
size_t i = 0; i < saddle2ToMaxima.size(); ++i) {
718 auto &maxs = saddle2ToMaxima[i];
720 TTK_PSORT(this->threadNumber_, maxs.begin(), maxs.end(),
730 const auto last = std::unique(maxs.begin(), maxs.end());
731 maxs.erase(last, maxs.end());
735 if(maxs.size() != 2) {
739 const auto s2 = criticalSaddles[i];
740 if(!pairedSaddles[s2]) {
741 sadMaxTriplets.emplace_back(
tripletType{s2, maxs[0], maxs[1]});
745 const auto nMinSadPairs = pairs.size();
747 tripletsToPersistencePairs(pairs, pairedMaxima, pairedSaddles, firstRep,
748 sadMaxTriplets, critSaddlesOrder.data(),
749 critMaxsOrder.data(), dim - 1);
751 const auto nSadMaxPairs = pairs.size() - nMinSadPairs;
754 "Computed " +
std::to_string(nSadMaxPairs) +
" saddle-max pairs", 1.0,
755 tm.getElapsedTime(), this->threadNumber_);
757 this->
printMsg(
"saddle-max pairs sequential part", 1.0,
762template <
typename triangulationType,
typename Container>
765 std::vector<bool> &onBoundary,
766 std::vector<Container> &s2Boundaries,
767 const std::vector<SimplexId> &s2Mapping,
768 const std::vector<SimplexId> &s1Mapping,
769 std::vector<SimplexId> &partners,
770 std::vector<Lock> &s1Locks,
771 std::vector<Lock> &s2Locks,
772 const triangulationType &triangulation)
const {
774 auto &boundaryIds{s2Boundaries[s2Mapping[s2]]};
776 const auto addBoundary = [&boundaryIds, &onBoundary](
const SimplexId e) {
779 boundaryIds.emplace(e);
780 onBoundary[e] =
true;
782 const auto it = boundaryIds.find(e);
783 boundaryIds.erase(it);
784 onBoundary[e] =
false;
788 const auto clearOnBoundary = [&boundaryIds, &onBoundary]() {
790 for(
const auto e : boundaryIds) {
791 onBoundary[e] =
false;
795 if(!boundaryIds.empty()) {
797 for(
const auto e : boundaryIds) {
798 onBoundary[e] =
true;
804 triangulation.getTriangleEdge(s2, i, e);
811 s2Locks[s2Mapping[s2]].lock();
813 while(!boundaryIds.empty()) {
815 const auto tau{*boundaryIds.begin()};
817 auto pTau{this->dg_.getPairedCell(
Cell{1, tau}, triangulation)};
818 bool critical{
false};
822#ifdef TTK_ENABLE_OPENMP
823#pragma omp atomic read
825 pTau = partners[tau];
826 if(pTau == -1 || s2Boundaries[s2Mapping[pTau]].empty()) {
829 }
while(*s2Boundaries[s2Mapping[pTau]].
begin() != tau);
838 s1Locks[s1Mapping[tau]].lock();
839 const auto cap = partners[tau];
840 if(partners[tau] == -1) {
843 s1Locks[s1Mapping[tau]].unlock();
847 s2Locks[s2Mapping[s2]].unlock();
851 return this->eliminateBoundariesSandwich(
852 s2, onBoundary, s2Boundaries, s2Mapping, s1Mapping, partners, s1Locks,
853 s2Locks, triangulation);
858 if(critical && s2Mapping[pTau] != -1) {
859 if(s2Mapping[pTau] < s2Mapping[s2]) {
865 s2Locks[s2Mapping[pTau]].lock();
866 for(
const auto e : s2Boundaries[s2Mapping[pTau]]) {
869 s2Locks[s2Mapping[pTau]].unlock();
870 if(this->Compute2SaddlesChildren) {
871 this->s2Children_[s2Mapping[s2]].emplace_back(s2Mapping[pTau]);
874 }
else if(s2Mapping[pTau] > s2Mapping[s2]) {
878 s1Locks[s1Mapping[tau]].lock();
879 const auto cap = partners[tau];
880 if(partners[tau] == pTau) {
883 s1Locks[s1Mapping[tau]].unlock();
888 s2Locks[s2Mapping[s2]].unlock();
889 return this->eliminateBoundariesSandwich(
890 pTau, onBoundary, s2Boundaries, s2Mapping, s1Mapping, partners,
891 s1Locks, s2Locks, triangulation);
898 triangulation.getTriangleEdge(pTau, i, e);
907 s2Locks[s2Mapping[s2]].unlock();
911template <
typename triangulationType>
913 std::vector<PersistencePair> &pairs,
914 std::vector<bool> &paired1Saddles,
915 std::vector<bool> &paired2Saddles,
916 const bool exportBoundaries,
917 std::vector<GeneratorType> &boundaries,
918 const std::vector<SimplexId> &critical1Saddles,
919 const std::vector<SimplexId> &critical2Saddles,
920 const std::vector<SimplexId> &crit1SaddlesOrder,
921 const triangulationType &triangulation)
const {
924 const auto nSadExtrPairs = pairs.size();
927 std::vector<SimplexId> saddles1{}, saddles2{};
929 for(
const auto s1 : critical1Saddles) {
930 if(!paired1Saddles[s1]) {
931 saddles1.emplace_back(s1);
935 for(
const auto s2 : critical2Saddles) {
936 if(!paired2Saddles[s2]) {
937 saddles2.emplace_back(s2);
941 if(this->Compute2SaddlesChildren) {
942 this->s2Children_.resize(saddles2.size());
946 const auto &edgesFiltrOrder{crit1SaddlesOrder};
948 auto &onBoundary{this->onBoundary_};
949 auto &edgeTrianglePartner{this->edgeTrianglePartner_};
953 return edgesFiltrOrder[a] > edgesFiltrOrder[b];
955 using Container = std::set<
SimplexId,
decltype(cmpEdges)>;
956 std::vector<Container> s2Boundaries(saddles2.size(), Container(cmpEdges));
959 auto &s2Mapping{this->s2Mapping_};
960#ifdef TTK_ENABLE_OPENMP
961#pragma omp parallel for num_threads(threadNumber_)
963 for(
size_t i = 0; i < saddles2.size(); ++i) {
964 s2Mapping[saddles2[i]] = i;
968 auto &s1Mapping{this->s1Mapping_};
969#ifdef TTK_ENABLE_OPENMP
970#pragma omp parallel for num_threads(threadNumber_)
972 for(
size_t i = 0; i < saddles1.size(); ++i) {
973 s1Mapping[saddles1[i]] = i;
977 std::vector<Lock> s1Locks(saddles1.size());
979 std::vector<Lock> s2Locks(saddles2.size());
983#ifdef TTK_ENABLE_OPENMP4
984#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
985 firstprivate(onBoundary)
987 for(
size_t i = 0; i < saddles2.size(); ++i) {
989 const auto s2 = saddles2[i];
990 this->eliminateBoundariesSandwich(s2, onBoundary, s2Boundaries, s2Mapping,
991 s1Mapping, edgeTrianglePartner, s1Locks,
992 s2Locks, triangulation);
998 for(
size_t i = 0; i < saddles2.size(); ++i) {
999 if(!s2Boundaries[i].empty()) {
1000 const auto s2 = saddles2[i];
1001 const auto s1 = *s2Boundaries[i].begin();
1003 pairs.emplace_back(s1, s2, 1);
1004 paired1Saddles[s1] =
true;
1005 paired2Saddles[s2] =
true;
1009 if(exportBoundaries) {
1010 boundaries.resize(s2Boundaries.size());
1011 for(
size_t i = 0; i < boundaries.size(); ++i) {
1012 const auto &boundSet{s2Boundaries[i]};
1013 if(boundSet.empty()) {
1017 {boundSet.begin(), boundSet.end()},
1019 std::array<SimplexId, 2>{
1020 this->dg_.getCellGreaterVertex(
Cell{2, saddles2[i]}, triangulation),
1021 this->dg_.getCellGreaterVertex(
1022 Cell{1, *boundSet.begin()}, triangulation),
1027 const auto nSadSadPairs = pairs.size() - nSadExtrPairs;
1030 "Computed " +
std::to_string(nSadSadPairs) +
" saddle-saddle pairs", 1.0,
1031 tm2.getElapsedTime(), this->threadNumber_);
1033 this->
printMsg(
"saddle-saddle pairs sequential part", 1.0,
1038template <
typename triangulationType>
1040 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
1041 std::array<std::vector<SimplexId>, 4> &critCellsOrder,
1043 const triangulationType &triangulation,
1044 const bool sortEdges)
const {
1048 this->dg_.getCriticalPoints(criticalCellsByDim, triangulation);
1050 this->
printMsg(
"Extracted critical cells", 1.0, tm.getElapsedTime(),
1055 auto &critEdges{this->critEdges_};
1057 critEdges.resize(criticalCellsByDim[1].size());
1059 std::vector<TriangleSimplex> critTriangles(criticalCellsByDim[2].size());
1060 std::vector<TetraSimplex> critTetras(criticalCellsByDim[3].size());
1062#ifdef TTK_ENABLE_OPENMP
1063#pragma omp parallel num_threads(threadNumber_)
1067#ifdef TTK_ENABLE_OPENMP
1068#pragma omp for nowait
1070 for(
size_t i = 0; i < critEdges.size(); ++i) {
1071 critEdges[i].fillEdge(i, offsets, triangulation);
1074#ifdef TTK_ENABLE_OPENMP
1075#pragma omp for nowait
1077 for(
size_t i = 0; i < critEdges.size(); ++i) {
1078 critEdges[i].fillEdge(criticalCellsByDim[1][i], offsets, triangulation);
1082#ifdef TTK_ENABLE_OPENMP
1083#pragma omp for nowait
1085 for(
size_t i = 0; i < critTriangles.size(); ++i) {
1086 critTriangles[i].fillTriangle(
1087 criticalCellsByDim[2][i], offsets, triangulation);
1090#ifdef TTK_ENABLE_OPENMP
1093 for(
size_t i = 0; i < critTetras.size(); ++i) {
1094 critTetras[i].fillTetra(criticalCellsByDim[3][i], offsets, triangulation);
1098 TTK_PSORT(this->threadNumber_, critEdges.begin(), critEdges.end());
1099 TTK_PSORT(this->threadNumber_, critTriangles.begin(), critTriangles.end());
1100 TTK_PSORT(this->threadNumber_, critTetras.begin(), critTetras.end());
1102#ifdef TTK_ENABLE_OPENMP
1103#pragma omp parallel num_threads(threadNumber_)
1106#ifdef TTK_ENABLE_OPENMP
1107#pragma omp for nowait
1109 for(
size_t i = 0; i < critEdges.size(); ++i) {
1110 critCellsOrder[1][critEdges[i].id_] = i;
1113#ifdef TTK_ENABLE_OPENMP
1114#pragma omp for nowait
1116 for(
size_t i = 0; i < critTriangles.size(); ++i) {
1117 criticalCellsByDim[2][i] = critTriangles[i].id_;
1118 critCellsOrder[2][critTriangles[i].id_] = i;
1121#ifdef TTK_ENABLE_OPENMP
1124 for(
size_t i = 0; i < critTetras.size(); ++i) {
1125 criticalCellsByDim[3][i] = critTetras[i].id_;
1126 critCellsOrder[3][critTetras[i].id_] = i;
1132 criticalCellsByDim[1].
end(),
1134 return critCellsOrder[1][a] < critCellsOrder[1][b];
1137#ifdef TTK_ENABLE_OPENMP
1138#pragma omp parallel for num_threads(threadNumber_)
1140 for(
size_t i = 0; i < critEdges.size(); ++i) {
1141 criticalCellsByDim[1][i] = critEdges[i].id_;
1145 this->
printMsg(
"Extracted & sorted critical cells", 1.0, tm.getElapsedTime(),
1150template <
typename triangulationType>
1152 std::vector<PersistencePair> &pairs,
1154 const triangulationType &triangulation,
1155 const bool ignoreBoundary,
1156 const bool compute2SaddlesChildren) {
1159 this->alloc(triangulation);
1163 const auto dim = this->dg_.getDimensionality();
1164 this->Compute2SaddlesChildren = compute2SaddlesChildren;
1167 std::array<std::vector<SimplexId>, 4> criticalCellsByDim{};
1169 auto &critCellsOrder{this->critCellsOrder_};
1171 this->extractCriticalCells(
1172 criticalCellsByDim, critCellsOrder, offsets, triangulation, dim == 3);
1175 auto &pairedMinima{this->pairedCritCells_[0]};
1177 auto &paired1Saddles{this->pairedCritCells_[1]};
1179 auto &paired2Saddles{this->pairedCritCells_[dim - 1]};
1181 auto &pairedMaxima{this->pairedCritCells_[dim]};
1186 if(this->ComputeMinSad) {
1188 this->getMinSaddlePairs(pairs, pairedMinima, paired1Saddles,
1189 criticalCellsByDim[1], critCellsOrder[1], offsets,
1193 for(
const auto min : criticalCellsByDim[0]) {
1194 if(!pairedMinima[min]) {
1195 pairs.emplace_back(min, -1, 0);
1196 pairedMinima[min] =
true;
1202 const auto globMin{*std::min_element(
1203 criticalCellsByDim[0].
begin(), criticalCellsByDim[0].
end(),
1205 return offsets[a] < offsets[b];
1207 pairs.emplace_back(globMin, -1, 0);
1208 pairedMinima[globMin] =
true;
1212 if(dim > 1 && this->ComputeSadMax) {
1214 this->getMaxSaddlePairs(
1215 pairs, pairedMaxima, paired2Saddles, criticalCellsByDim[dim - 1],
1216 critCellsOrder[dim - 1], critCellsOrder[dim], triangulation);
1219 if(ignoreBoundary) {
1223 = std::find_if(pairs.begin(), pairs.end(), [&](
const PersistencePair &p) {
1224 if(p.type < dim - 1) {
1227 const Cell cmax{dim, p.death};
1228 const auto vmax{this->getCellGreaterVertex(cmax, triangulation)};
1229 return offsets[vmax] == triangulation.getNumberOfVertices() - 1;
1232 if(it != pairs.end()) {
1234 paired2Saddles[it->birth] =
false;
1235 pairedMaxima[it->death] =
false;
1241 if(dim == 3 && !criticalCellsByDim[1].empty()
1242 && !criticalCellsByDim[2].empty() && this->ComputeSadSad) {
1243 std::vector<GeneratorType> tmp{};
1244 this->getSaddleSaddlePairs(
1245 pairs, paired1Saddles, paired2Saddles,
false, tmp, criticalCellsByDim[1],
1246 criticalCellsByDim[2], critCellsOrder[1], triangulation);
1249 if(std::is_same<triangulationType, ttk::ExplicitTriangulation>::value) {
1251 size_t nHandles{}, nCavities{}, nNonPairedMax{};
1252 if((dim == 2 && !ignoreBoundary && this->ComputeMinSad
1253 && this->ComputeSadMax)
1254 || (dim == 3 && this->ComputeMinSad && this->ComputeSadSad)) {
1256 for(
const auto s1 : criticalCellsByDim[1]) {
1257 if(!paired1Saddles[s1]) {
1258 paired1Saddles[s1] =
true;
1260 pairs.emplace_back(s1, -1, 1);
1265 if(dim == 3 && !ignoreBoundary && this->ComputeSadMax
1266 && this->ComputeSadSad) {
1268 for(
const auto s2 : criticalCellsByDim[2]) {
1269 if(!paired2Saddles[s2]) {
1270 paired2Saddles[s2] =
true;
1272 pairs.emplace_back(s2, -1, 2);
1277 if(dim == 2 && !ignoreBoundary && this->ComputeSadMax) {
1279 for(
const auto max : criticalCellsByDim[dim]) {
1280 if(!pairedMaxima[max]) {
1281 pairs.emplace_back(max, -1, 2);
1288 = (dim == 3 ? nCavities : nHandles) + nConnComp - nNonPairedMax;
1289 nBoundComp = std::max(nBoundComp, 0);
1292 const std::vector<std::vector<std::string>> rows{
1303 "Computed " +
std::to_string(pairs.size()) +
" persistence pairs", 1.0,
1304 tm.getElapsedTime(), this->threadNumber_);
1306 this->displayStats(pairs, criticalCellsByDim, pairedMinima, paired1Saddles,
1307 paired2Saddles, pairedMaxima);
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual int setThreadNumber(const int threadNumber)
Minimalist debugging class.
virtual int setDebugLevel(const int &debugLevel)
TTK DiscreteMorseSandwich processing package.
void setComputeSadMax(const bool data)
std::array< SimplexId, 3 > tripletType
Triplet type for persistence pairs.
std::vector< SimplexId > s1Mapping_
ttk::dcg::DiscreteGradient && getGradient()
dcg::DiscreteGradient dg_
std::vector< SimplexId > firstRepMin_
std::vector< std::vector< SimplexId > > s2Children_
const std::vector< std::vector< SimplexId > > & get2SaddlesChildren() const
std::array< std::vector< SimplexId >, 4 > critCellsOrder_
SimplexId getCellGreaterVertex(const dcg::Cell &c, const triangulationType &triangulation)
std::vector< bool > onBoundary_
void alloc(const triangulationType &triangulation)
std::vector< EdgeSimplex > critEdges_
std::vector< SimplexId > edgeTrianglePartner_
void setGradient(ttk::dcg::DiscreteGradient &&dg)
Ugly hack to avoid a call to buildGradient()
void setComputeMinSad(const bool data)
int buildGradient(const void *const scalars, const size_t scalarsMTime, const SimplexId *const offsets, const triangulationType &triangulation)
void setInputOffsets(const SimplexId *const offsets)
bool Compute2SaddlesChildren
void preconditionTriangulation(AbstractTriangulation *const data)
std::vector< SimplexId > firstRepMax_
void setComputeSadSad(const bool data)
void getMinSaddlePairs(std::vector< PersistencePair > &pairs, std::vector< bool > &pairedMinima, std::vector< bool > &paired1Saddles, const std::vector< SimplexId > &criticalEdges, const std::vector< SimplexId > &critEdgesOrder, const SimplexId *const offsets, const triangulationType &triangulation) const
Compute the pairs of dimension 0.
void extractCriticalCells(std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, std::array< std::vector< SimplexId >, 4 > &critCellsOrder, const SimplexId *const offsets, const triangulationType &triangulation, const bool sortEdges) const
Extract & sort critical cell from the DiscreteGradient.
void displayStats(const std::vector< PersistencePair > &pairs, const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const std::vector< bool > &pairedMinima, const std::vector< bool > &paired1Saddles, const std::vector< bool > &paired2Saddles, const std::vector< bool > &pairedMaxima) const
Print number of pairs, critical cells per dimension & unpaired cells.
std::vector< SimplexId > s2Mapping_
void getSaddleSaddlePairs(std::vector< PersistencePair > &pairs, std::vector< bool > &paired1Saddles, std::vector< bool > &paired2Saddles, const bool exportBoundaries, std::vector< GeneratorType > &boundaries, const std::vector< SimplexId > &critical1Saddles, const std::vector< SimplexId > &critical2Saddles, const std::vector< SimplexId > &crit1SaddlesOrder, const triangulationType &triangulation) const
Compute the saddle-saddle pairs (in 3D)
int computePersistencePairs(std::vector< PersistencePair > &pairs, const SimplexId *const offsets, const triangulationType &triangulation, const bool ignoreBoundary, const bool compute2SaddlesChildren=false)
Compute the persistence pairs from the discrete gradient.
std::vector< std::vector< SimplexId > > getSaddle2ToMaxima(const std::vector< SimplexId > &criticalCells, const GFS &getFaceStar, const GFSN &getFaceStarNumber, const OB &isOnBoundary, const triangulationType &triangulation) const
Follow the ascending 1-separatrices to compute the saddles -> maxima association.
std::array< std::vector< bool >, 4 > pairedCritCells_
SimplexId eliminateBoundariesSandwich(const SimplexId s2, std::vector< bool > &onBoundary, std::vector< Container > &s2Boundaries, const std::vector< SimplexId > &s2Mapping, const std::vector< SimplexId > &s1Mapping, std::vector< SimplexId > &partners, std::vector< Lock > &s1Locks, std::vector< Lock > &s2Locks, const triangulationType &triangulation) const
Detect 1-saddles paired to a given 2-saddle.
void getMaxSaddlePairs(std::vector< PersistencePair > &pairs, std::vector< bool > &pairedMaxima, std::vector< bool > &pairedSaddles, const std::vector< SimplexId > &criticalSaddles, const std::vector< SimplexId > &critSaddlesOrder, const std::vector< SimplexId > &critMaxsOrder, const triangulationType &triangulation) const
Compute the pairs of dimension dim - 1.
void tripletsToPersistencePairs(std::vector< PersistencePair > &pairs, std::vector< bool > &pairedExtrema, std::vector< bool > &pairedSaddles, std::vector< SimplexId > &reps, std::vector< tripletType > &triplets, const SimplexId *const saddlesOrder, const SimplexId *const extremaOrder, const SimplexId pairDim) const
Compute persistence pairs from triplets.
std::vector< std::vector< SimplexId > > getSaddle1ToMinima(const std::vector< SimplexId > &criticalEdges, const triangulationType &triangulation) const
Follow the descending 1-separatrices to compute the saddles -> minima association.
TTK discreteGradient processing package.
void setInputOffsets(const SimplexId *const data)
void setInputScalarField(const void *const data, const size_t mTime)
SimplexId getNumberOfCells(const int dimension, const triangulationType &triangulation) const
void setLocalGradient()
Use local storage instead of cache.
SimplexId getCellGreaterVertex(const Cell c, const triangulationType &triangulation) const
void preconditionTriangulation(AbstractTriangulation *const data)
int buildGradient(const triangulationType &triangulation, bool bypassCache=false)
int getDimensionality() const
int getDescendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation) const
std::string to_string(__int128)
int SimplexId
Identifier type for simplices of any dimension.
T end(std::pair< T, T > &p)
T begin(std::pair< T, T > &p)
Simplex adaptation for edges
void fillEdge(const SimplexId id, const SimplexId *const offsets, const triangulationType &triangulation)
Type for exporting persistent generators.
std::array< SimplexId, 2 > critVertsIds
std::vector< SimplexId > boundary
Persistence pair struct as exported by DiscreteGradient.
PersistencePair(SimplexId b, SimplexId d, int t)
Ad-hoc struct for sorting simplices.
friend bool operator<(const Simplex< n > &lhs, const Simplex< n > &rhs)
std::array< SimplexId, n > vertsOrder_
Simplex adaptation for tetrahedra
void fillTetra(const SimplexId id, const SimplexId *const offsets, const triangulationType &triangulation)
Simplex adaptation for triangles
void fillTriangle(const SimplexId id, const SimplexId *const offsets, const triangulationType &triangulation)
printMsg(debug::output::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)