TTK
Loading...
Searching...
No Matches
DiscreteMorseSandwich.h
Go to the documentation of this file.
1
21
22#pragma once
23
24#include <DiscreteGradient.h>
25
26#include <algorithm>
27#include <numeric>
28
29namespace ttk {
30 class DiscreteMorseSandwich : virtual public Debug {
31 public:
33
43 int type;
44
46 : birth{b}, death{d}, type{t} {
47 }
48 };
49
52 }
53
54 inline void setInputOffsets(const SimplexId *const offsets) {
55 this->dg_.setInputOffsets(offsets);
56 }
57
58 inline void setComputeMinSad(const bool data) {
59 this->ComputeMinSad = data;
60 }
61 inline void setComputeSadSad(const bool data) {
62 this->ComputeSadSad = data;
63 }
64 inline void setComputeSadMax(const bool data) {
65 this->ComputeSadMax = data;
66 }
67
68 template <typename triangulationType>
69 inline int buildGradient(const void *const scalars,
70 const size_t scalarsMTime,
71 const SimplexId *const offsets,
72 const triangulationType &triangulation) {
73 this->dg_.setDebugLevel(this->debugLevel_);
75 this->dg_.setInputOffsets(offsets);
76 this->dg_.setInputScalarField(scalars, scalarsMTime);
77 return this->dg_.buildGradient(triangulation);
78 }
79
92 this->dg_ = std::move(dg);
93 // reset gradient pointer to local storage
94 this->dg_.setLocalGradient();
95 }
97 return std::move(this->dg_);
98 }
99
100 template <typename triangulationType>
101 inline SimplexId
103 const triangulationType &triangulation) {
104 return this->dg_.getCellGreaterVertex(c, triangulation);
105 }
106
107 inline const std::vector<std::vector<SimplexId>> &
109 return this->s2Children_;
110 }
111
126 template <typename triangulationType>
127 int computePersistencePairs(std::vector<PersistencePair> &pairs,
128 const SimplexId *const offsets,
129 const triangulationType &triangulation,
130 const bool ignoreBoundary,
131 const bool compute2SaddlesChildren = false);
132
141 std::vector<SimplexId> boundary;
146 std::array<SimplexId, 2> critVertsIds;
147 };
148
149 protected:
159 template <typename triangulationType>
160 std::vector<std::vector<SimplexId>>
161 getSaddle1ToMinima(const std::vector<SimplexId> &criticalEdges,
162 const triangulationType &triangulation) const;
163
179 template <typename triangulationType,
180 typename GFS,
181 typename GFSN,
182 typename OB>
183 std::vector<std::vector<SimplexId>>
184 getSaddle2ToMaxima(const std::vector<SimplexId> &criticalCells,
185 const GFS &getFaceStar,
186 const GFSN &getFaceStarNumber,
187 const OB &isOnBoundary,
188 const triangulationType &triangulation) const;
189
201 template <typename triangulationType>
202 void getMinSaddlePairs(std::vector<PersistencePair> &pairs,
203 std::vector<bool> &pairedMinima,
204 std::vector<bool> &paired1Saddles,
205 const std::vector<SimplexId> &criticalEdges,
206 const std::vector<SimplexId> &critEdgesOrder,
207 const SimplexId *const offsets,
208 const triangulationType &triangulation) const;
209
221 template <typename triangulationType>
222 void getMaxSaddlePairs(std::vector<PersistencePair> &pairs,
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;
229
243 template <typename triangulationType>
244 void getSaddleSaddlePairs(std::vector<PersistencePair> &pairs,
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;
253
263 template <typename triangulationType>
265 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
266 std::array<std::vector<SimplexId>, 4> &critCellsOrder,
267 const SimplexId *const offsets,
268 const triangulationType &triangulation,
269 const bool sortEdges) const;
270
282 void displayStats(
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;
289
297 using tripletType = std::array<SimplexId, 3>;
298
311 void tripletsToPersistencePairs(std::vector<PersistencePair> &pairs,
312 std::vector<bool> &pairedExtrema,
313 std::vector<bool> &pairedSaddles,
314 std::vector<SimplexId> &reps,
315 std::vector<tripletType> &triplets,
316 const SimplexId *const saddlesOrder,
317 const SimplexId *const extremaOrder,
318 const SimplexId pairDim) const;
319
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;
349
355 template <size_t n>
356 struct Simplex {
361 std::array<SimplexId, n> vertsOrder_{};
364 friend bool operator<(const Simplex<n> &lhs, const Simplex<n> &rhs) {
365 return lhs.vertsOrder_ < rhs.vertsOrder_;
366 }
367 };
368
372 struct EdgeSimplex : Simplex<2> {
373 template <typename triangulationType>
374 void fillEdge(const SimplexId id,
375 const SimplexId *const offsets,
376 const triangulationType &triangulation) {
377 this->id_ = id;
378 triangulation.getEdgeVertex(id, 0, this->vertsOrder_[0]);
379 triangulation.getEdgeVertex(id, 1, this->vertsOrder_[1]);
380 this->vertsOrder_[0] = offsets[this->vertsOrder_[0]];
381 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
382 // sort vertices in decreasing order
383 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
384 }
385 };
386
391 template <typename triangulationType>
392 void fillTriangle(const SimplexId id,
393 const SimplexId *const offsets,
394 const triangulationType &triangulation) {
395 this->id_ = id;
396 triangulation.getTriangleVertex(id, 0, this->vertsOrder_[0]);
397 triangulation.getTriangleVertex(id, 1, this->vertsOrder_[1]);
398 triangulation.getTriangleVertex(id, 2, this->vertsOrder_[2]);
399 this->vertsOrder_[0] = offsets[this->vertsOrder_[0]];
400 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
401 this->vertsOrder_[2] = offsets[this->vertsOrder_[2]];
402 // sort vertices in decreasing order
403 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
404 }
405 };
406
410 struct TetraSimplex : Simplex<4> {
411 template <typename triangulationType>
412 void fillTetra(const SimplexId id,
413 const SimplexId *const offsets,
414 const triangulationType &triangulation) {
415 this->id_ = id;
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]);
420 this->vertsOrder_[0] = offsets[this->vertsOrder_[0]];
421 this->vertsOrder_[1] = offsets[this->vertsOrder_[1]];
422 this->vertsOrder_[2] = offsets[this->vertsOrder_[2]];
423 this->vertsOrder_[3] = offsets[this->vertsOrder_[3]];
424 // sort vertices in decreasing order
425 std::sort(this->vertsOrder_.rbegin(), this->vertsOrder_.rend());
426 }
427 };
428
429 template <typename triangulationType>
430 void alloc(const triangulationType &triangulation) {
431 Timer tm{};
432 const auto dim{this->dg_.getDimensionality()};
433 if(dim > 3 || dim < 1) {
434 return;
435 }
436 this->firstRepMin_.resize(triangulation.getNumberOfVertices());
437 if(dim > 1) {
438 this->firstRepMax_.resize(triangulation.getNumberOfCells());
439 }
440 if(dim > 2) {
441 this->critEdges_.resize(triangulation.getNumberOfEdges());
442 this->edgeTrianglePartner_.resize(triangulation.getNumberOfEdges(), -1);
443 this->onBoundary_.resize(triangulation.getNumberOfEdges(), false);
444 this->s2Mapping_.resize(triangulation.getNumberOfTriangles(), -1);
445 this->s1Mapping_.resize(triangulation.getNumberOfEdges(), -1);
446 }
447 for(int i = 0; i < dim + 1; ++i) {
448 this->pairedCritCells_[i].resize(
449 this->dg_.getNumberOfCells(i, triangulation), false);
450 }
451 for(int i = 1; i < dim + 1; ++i) {
452 this->critCellsOrder_[i].resize(
453 this->dg_.getNumberOfCells(i, triangulation), -1);
454 }
455 this->printMsg("Memory allocations", 1.0, tm.getElapsedTime(), 1,
457 }
458
459 void clear() {
460 Timer tm{};
461 this->firstRepMin_ = {};
462 this->firstRepMax_ = {};
463 this->edgeTrianglePartner_ = {};
464 this->s2Mapping_ = {};
465 this->s1Mapping_ = {};
466 this->critEdges_ = {};
467 this->pairedCritCells_ = {};
468 this->onBoundary_ = {};
469 this->critCellsOrder_ = {};
470 this->printMsg("Memory cleanup", 1.0, tm.getElapsedTime(), 1,
472 }
473
475
476 // factor memory allocations outside computation loops
477 mutable std::vector<SimplexId> firstRepMin_{}, firstRepMax_{},
479 mutable std::vector<EdgeSimplex> critEdges_{};
480 mutable std::array<std::vector<bool>, 4> pairedCritCells_{};
481 mutable std::vector<bool> onBoundary_{};
482 mutable std::array<std::vector<SimplexId>, 4> critCellsOrder_{};
483 mutable std::vector<std::vector<SimplexId>> s2Children_{};
484
485 bool ComputeMinSad{true};
486 bool ComputeSadSad{true};
487 bool ComputeSadMax{true};
489 };
490} // namespace ttk
491
492template <typename triangulationType>
493std::vector<std::vector<SimplexId>>
495 const std::vector<SimplexId> &criticalEdges,
496 const triangulationType &triangulation) const {
497
498 Timer tm{};
499
500 std::vector<std::vector<SimplexId>> res(criticalEdges.size());
501
502 // follow vpaths from 1-saddles to minima
503#ifdef TTK_ENABLE_OPENMP
504#pragma omp parallel for num_threads(threadNumber_)
505#endif
506 for(size_t i = 0; i < criticalEdges.size(); ++i) {
507 auto &mins = res[i];
508
509 const auto followVPath = [this, &mins, &triangulation](const SimplexId v) {
510 std::vector<Cell> vpath{};
511 this->dg_.getDescendingPath(Cell{0, v}, vpath, triangulation);
512 Cell &lastCell = vpath.back();
513 if(lastCell.dim_ == 0 && this->dg_.isCellCritical(lastCell)) {
514 mins.emplace_back(lastCell.id_);
515 }
516 };
517
518 // critical edge vertices
519 SimplexId v0{}, v1{};
520 triangulation.getEdgeVertex(criticalEdges[i], 0, v0);
521 triangulation.getEdgeVertex(criticalEdges[i], 1, v1);
522
523 // follow vpath from each vertex of the critical edge
524 followVPath(v0);
525 followVPath(v1);
526 }
527
528 this->printMsg("Computed the descending 1-separatrices", 1.0,
529 tm.getElapsedTime(), this->threadNumber_, debug::LineMode::NEW,
531
532 return res;
533}
534
535template <typename triangulationType, typename GFS, typename GFSN, typename OB>
536std::vector<std::vector<SimplexId>>
538 const std::vector<SimplexId> &criticalCells,
539 const GFS &getFaceStar,
540 const GFSN &getFaceStarNumber,
541 const OB &isOnBoundary,
542 const triangulationType &triangulation) const {
543
544 Timer tm{};
545
546 const auto dim = this->dg_.getDimensionality();
547 std::vector<std::vector<SimplexId>> res(criticalCells.size());
548
549 // follow vpaths from 2-saddles to maxima
550#ifdef TTK_ENABLE_OPENMP
551#pragma omp parallel for num_threads(threadNumber_)
552#endif
553 for(size_t i = 0; i < criticalCells.size(); ++i) {
554 const auto sid = criticalCells[i];
555 auto &maxs = res[i];
556
557 const auto followVPath
558 = [this, dim, &maxs, &triangulation](const SimplexId v) {
559 std::vector<Cell> vpath{};
560 this->dg_.getAscendingPath(Cell{dim, v}, vpath, triangulation);
561 Cell &lastCell = vpath.back();
562 if(lastCell.dim_ == dim && this->dg_.isCellCritical(lastCell)) {
563 maxs.emplace_back(lastCell.id_);
564 } else if(lastCell.dim_ == dim - 1) {
565 maxs.emplace_back(-1);
566 }
567 };
568
569 const auto starNumber = getFaceStarNumber(sid);
570
571 for(SimplexId j = 0; j < starNumber; ++j) {
572 SimplexId cellId{};
573 getFaceStar(sid, j, cellId);
574 followVPath(cellId);
575 }
576
577 if(isOnBoundary(sid)) {
578 // critical saddle is on boundary
579 maxs.emplace_back(-1);
580 }
581 }
582
583 this->printMsg("Computed the ascending 1-separatrices", 1.0,
584 tm.getElapsedTime(), this->threadNumber_, debug::LineMode::NEW,
586
587 return res;
588}
589
590template <typename triangulationType>
592 std::vector<PersistencePair> &pairs,
593 std::vector<bool> &pairedMinima,
594 std::vector<bool> &paired1Saddles,
595 const std::vector<SimplexId> &criticalEdges,
596 const std::vector<SimplexId> &critEdgesOrder,
597 const SimplexId *const offsets,
598 const triangulationType &triangulation) const {
599
600 Timer tm{};
601
602 auto saddle1ToMinima = getSaddle1ToMinima(criticalEdges, triangulation);
603
604 Timer tmseq{};
605
606 auto &firstRep{this->firstRepMin_};
607 std::iota(firstRep.begin(), firstRep.end(), 0);
608 std::vector<tripletType> sadMinTriplets{};
609
610 for(size_t i = 0; i < saddle1ToMinima.size(); ++i) {
611 auto &mins = saddle1ToMinima[i];
612 const auto s1 = criticalEdges[i];
613 // remove duplicates
614 TTK_PSORT(this->threadNumber_, mins.begin(), mins.end());
615 const auto last = std::unique(mins.begin(), mins.end());
616 mins.erase(last, mins.end());
617 if(mins.size() != 2) {
618 continue;
619 }
620 sadMinTriplets.emplace_back(tripletType{s1, mins[0], mins[1]});
621 }
622
623 tripletsToPersistencePairs(pairs, pairedMinima, paired1Saddles, firstRep,
624 sadMinTriplets, critEdgesOrder.data(), offsets, 0);
625
626 const auto nMinSadPairs = pairs.size();
627
628 this->printMsg(
629 "Computed " + std::to_string(nMinSadPairs) + " min-saddle pairs", 1.0,
630 tm.getElapsedTime(), this->threadNumber_);
631
632 this->printMsg("min-saddle pairs sequential part", 1.0,
633 tmseq.getElapsedTime(), 1, debug::LineMode::NEW,
635}
636
637template <typename triangulationType>
639 std::vector<PersistencePair> &pairs,
640 std::vector<bool> &pairedMaxima,
641 std::vector<bool> &pairedSaddles,
642 const std::vector<SimplexId> &criticalSaddles,
643 const std::vector<SimplexId> &critSaddlesOrder,
644 const std::vector<SimplexId> &critMaxsOrder,
645 const triangulationType &triangulation) const {
646
647 Timer tm{};
648
649 const auto dim = this->dg_.getDimensionality();
650
651 auto saddle2ToMaxima
652 = dim == 3
653 ? getSaddle2ToMaxima(
654 criticalSaddles,
655 [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) {
656 return triangulation.getTriangleStar(a, i, r);
657 },
658 [&triangulation](const SimplexId a) {
659 return triangulation.getTriangleStarNumber(a);
660 },
661 [&triangulation](const SimplexId a) {
662 return triangulation.isTriangleOnBoundary(a);
663 },
664 triangulation)
665 : getSaddle2ToMaxima(
666 criticalSaddles,
667 [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) {
668 return triangulation.getEdgeStar(a, i, r);
669 },
670 [&triangulation](const SimplexId a) {
671 return triangulation.getEdgeStarNumber(a);
672 },
673 [&triangulation](const SimplexId a) {
674 return triangulation.isEdgeOnBoundary(a);
675 },
676 triangulation);
677
678 Timer tmseq{};
679
680 auto &firstRep{this->firstRepMax_};
681 std::iota(firstRep.begin(), firstRep.end(), 0);
682 std::vector<tripletType> sadMaxTriplets{};
683
684 for(size_t i = 0; i < saddle2ToMaxima.size(); ++i) {
685 auto &maxs = saddle2ToMaxima[i];
686 // remove duplicates
687 TTK_PSORT(this->threadNumber_, maxs.begin(), maxs.end(),
688 [](const SimplexId a, const SimplexId b) {
689 // positive values (actual maxima) before negative ones
690 // (boundary component id)
691 if(a * b >= 0) {
692 return a < b;
693 } else {
694 return a > b;
695 }
696 });
697 const auto last = std::unique(maxs.begin(), maxs.end());
698 maxs.erase(last, maxs.end());
699
700 // remove "doughnut" configurations: two ascending separatrices
701 // leading to the same maximum/boundary component
702 if(maxs.size() != 2) {
703 continue;
704 }
705
706 const auto s2 = criticalSaddles[i];
707 if(!pairedSaddles[s2]) {
708 sadMaxTriplets.emplace_back(tripletType{s2, maxs[0], maxs[1]});
709 }
710 }
711
712 const auto nMinSadPairs = pairs.size();
713
714 tripletsToPersistencePairs(pairs, pairedMaxima, pairedSaddles, firstRep,
715 sadMaxTriplets, critSaddlesOrder.data(),
716 critMaxsOrder.data(), dim - 1);
717
718 const auto nSadMaxPairs = pairs.size() - nMinSadPairs;
719
720 this->printMsg(
721 "Computed " + std::to_string(nSadMaxPairs) + " saddle-max pairs", 1.0,
722 tm.getElapsedTime(), this->threadNumber_);
723
724 this->printMsg("saddle-max pairs sequential part", 1.0,
725 tmseq.getElapsedTime(), 1, debug::LineMode::NEW,
727}
728
729template <typename triangulationType, typename Container>
731 const SimplexId s2,
732 std::vector<bool> &onBoundary,
733 std::vector<Container> &s2Boundaries,
734 const std::vector<SimplexId> &s2Mapping,
735 const std::vector<SimplexId> &s1Mapping,
736 std::vector<SimplexId> &partners,
737 std::vector<Lock> &s1Locks,
738 std::vector<Lock> &s2Locks,
739 const triangulationType &triangulation) const {
740
741 auto &boundaryIds{s2Boundaries[s2Mapping[s2]]};
742
743 const auto addBoundary = [&boundaryIds, &onBoundary](const SimplexId e) {
744 // add edge e to boundaryIds/onBoundary modulo 2
745 if(!onBoundary[e]) {
746 boundaryIds.emplace(e);
747 onBoundary[e] = true;
748 } else {
749 const auto it = boundaryIds.find(e);
750 boundaryIds.erase(it);
751 onBoundary[e] = false;
752 }
753 };
754
755 const auto clearOnBoundary = [&boundaryIds, &onBoundary]() {
756 // clear the onBoundary vector (set everything to false)
757 for(const auto e : boundaryIds) {
758 onBoundary[e] = false;
759 }
760 };
761
762 if(!boundaryIds.empty()) {
763 // restore previously computed s2 boundary
764 for(const auto e : boundaryIds) {
765 onBoundary[e] = true;
766 }
767 } else {
768 // init cascade with s2 triangle boundary (3 edges)
769 for(SimplexId i = 0; i < 3; ++i) {
770 SimplexId e{};
771 triangulation.getTriangleEdge(s2, i, e);
772 addBoundary(e);
773 }
774 }
775
776 // lock the 2-saddle to ensure that only one thread can perform the
777 // boundary expansion
778 s2Locks[s2Mapping[s2]].lock();
779
780 while(!boundaryIds.empty()) {
781 // tau: youngest edge on boundary
782 const auto tau{*boundaryIds.begin()};
783 // use the Discrete Gradient to find a triangle paired to tau
784 auto pTau{this->dg_.getPairedCell(Cell{1, tau}, triangulation)};
785 bool critical{false};
786 if(pTau == -1) {
787 // maybe tau is critical and paired to a critical triangle
788 do {
789#ifdef TTK_ENABLE_OPENMP
790#pragma omp atomic read
791#endif // TTK_ENABLE_OPENMP
792 pTau = partners[tau];
793 if(pTau == -1 || s2Boundaries[s2Mapping[pTau]].empty()) {
794 break;
795 }
796 } while(*s2Boundaries[s2Mapping[pTau]].begin() != tau);
797
798 critical = true;
799 }
800 if(pTau == -1) {
801 // tau is critical and not paired
802
803 // compare-and-swap from "Towards Lockfree Persistent Homology"
804 // using locks over 1-saddles instead of atomics (OpenMP compatibility)
805 s1Locks[s1Mapping[tau]].lock();
806 const auto cap = partners[tau];
807 if(partners[tau] == -1) {
808 partners[tau] = s2;
809 }
810 s1Locks[s1Mapping[tau]].unlock();
811
812 // cleanup before exiting
813 clearOnBoundary();
814 s2Locks[s2Mapping[s2]].unlock();
815 if(cap == -1) {
816 return tau;
817 } else {
818 return this->eliminateBoundariesSandwich(
819 s2, onBoundary, s2Boundaries, s2Mapping, s1Mapping, partners, s1Locks,
820 s2Locks, triangulation);
821 }
822
823 } else {
824 // expand boundary
825 if(critical && s2Mapping[pTau] != -1) {
826 if(s2Mapping[pTau] < s2Mapping[s2]) {
827 // pTau is an already-paired 2-saddle
828 // merge pTau boundary into s2 boundary
829
830 // make sure that pTau boundary is not modified by another
831 // thread while we merge the two boundaries...
832 s2Locks[s2Mapping[pTau]].lock();
833 for(const auto e : s2Boundaries[s2Mapping[pTau]]) {
834 addBoundary(e);
835 }
836 s2Locks[s2Mapping[pTau]].unlock();
837 if(this->Compute2SaddlesChildren) {
838 this->s2Children_[s2Mapping[s2]].emplace_back(s2Mapping[pTau]);
839 }
840
841 } else if(s2Mapping[pTau] > s2Mapping[s2]) {
842
843 // compare-and-swap from "Towards Lockfree Persistent
844 // Homology" using locks over 1-saddles
845 s1Locks[s1Mapping[tau]].lock();
846 const auto cap = partners[tau];
847 if(partners[tau] == pTau) {
848 partners[tau] = s2;
849 }
850 s1Locks[s1Mapping[tau]].unlock();
851
852 if(cap == pTau) {
853 // cleanup before exiting
854 clearOnBoundary();
855 s2Locks[s2Mapping[s2]].unlock();
856 return this->eliminateBoundariesSandwich(
857 pTau, onBoundary, s2Boundaries, s2Mapping, s1Mapping, partners,
858 s1Locks, s2Locks, triangulation);
859 }
860 }
861 } else { // pTau is a regular triangle
862 // add pTau triangle boundary (3 edges)
863 for(SimplexId i = 0; i < 3; ++i) {
864 SimplexId e{};
865 triangulation.getTriangleEdge(pTau, i, e);
866 addBoundary(e);
867 }
868 }
869 }
870 }
871
872 // cleanup before exiting
873 clearOnBoundary();
874 s2Locks[s2Mapping[s2]].unlock();
875 return -1;
876}
877
878template <typename triangulationType>
880 std::vector<PersistencePair> &pairs,
881 std::vector<bool> &paired1Saddles,
882 std::vector<bool> &paired2Saddles,
883 const bool exportBoundaries,
884 std::vector<GeneratorType> &boundaries,
885 const std::vector<SimplexId> &critical1Saddles,
886 const std::vector<SimplexId> &critical2Saddles,
887 const std::vector<SimplexId> &crit1SaddlesOrder,
888 const triangulationType &triangulation) const {
889
890 Timer tm2{};
891 const auto nSadExtrPairs = pairs.size();
892
893 // 1- and 2-saddles yet to be paired
894 std::vector<SimplexId> saddles1{}, saddles2{};
895 // filter out already paired 1-saddles (edge id)
896 for(const auto s1 : critical1Saddles) {
897 if(!paired1Saddles[s1]) {
898 saddles1.emplace_back(s1);
899 }
900 }
901 // filter out already paired 2-saddles (triangle id)
902 for(const auto s2 : critical2Saddles) {
903 if(!paired2Saddles[s2]) {
904 saddles2.emplace_back(s2);
905 }
906 }
907
908 Timer tmpar{};
909
910 if(this->Compute2SaddlesChildren) {
911 this->s2Children_.resize(saddles2.size());
912 }
913
914 // sort every triangulation edges by filtration order
915 const auto &edgesFiltrOrder{crit1SaddlesOrder};
916
917 auto &onBoundary{this->onBoundary_};
918 auto &edgeTrianglePartner{this->edgeTrianglePartner_};
919
920 const auto cmpEdges
921 = [&edgesFiltrOrder](const SimplexId a, const SimplexId b) {
922 return edgesFiltrOrder[a] > edgesFiltrOrder[b];
923 };
924 using Container = std::set<SimplexId, decltype(cmpEdges)>;
925 std::vector<Container> s2Boundaries(saddles2.size(), Container(cmpEdges));
926
927 // unpaired critical triangle id -> index in saddle2 vector
928 auto &s2Mapping{this->s2Mapping_};
929#ifdef TTK_ENABLE_OPENMP
930#pragma omp parallel for num_threads(threadNumber_)
931#endif // TTK_ENABLE_OPENMP
932 for(size_t i = 0; i < saddles2.size(); ++i) {
933 s2Mapping[saddles2[i]] = i;
934 }
935
936 // unpaired critical edge id -> index in saddle1 vector
937 auto &s1Mapping{this->s1Mapping_};
938#ifdef TTK_ENABLE_OPENMP
939#pragma omp parallel for num_threads(threadNumber_)
940#endif // TTK_ENABLE_OPENMP
941 for(size_t i = 0; i < saddles1.size(); ++i) {
942 s1Mapping[saddles1[i]] = i;
943 }
944
945 // one lock per 1-saddle
946 std::vector<Lock> s1Locks(saddles1.size());
947 // one lock per 2-saddle
948 std::vector<Lock> s2Locks(saddles2.size());
949
950 // compute 2-saddles boundaries in parallel
951
952#ifdef TTK_ENABLE_OPENMP
953#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
954 firstprivate(onBoundary)
955#endif // TTK_ENABLE_OPENMP
956 for(size_t i = 0; i < saddles2.size(); ++i) {
957 // 2-saddles sorted in increasing order
958 const auto s2 = saddles2[i];
959 this->eliminateBoundariesSandwich(s2, onBoundary, s2Boundaries, s2Mapping,
960 s1Mapping, edgeTrianglePartner, s1Locks,
961 s2Locks, triangulation);
962 }
963
964 Timer tmseq{};
965
966 // extract saddle-saddle pairs from computed boundaries
967 for(size_t i = 0; i < saddles2.size(); ++i) {
968 if(!s2Boundaries[i].empty()) {
969 const auto s2 = saddles2[i];
970 const auto s1 = *s2Boundaries[i].begin();
971 // we found a pair
972 pairs.emplace_back(s1, s2, 1);
973 paired1Saddles[s1] = true;
974 paired2Saddles[s2] = true;
975 }
976 }
977
978 if(exportBoundaries) {
979 boundaries.resize(s2Boundaries.size());
980 for(size_t i = 0; i < boundaries.size(); ++i) {
981 const auto &boundSet{s2Boundaries[i]};
982 if(boundSet.empty()) {
983 continue;
984 }
985 boundaries[i] = {
986 {boundSet.begin(), boundSet.end()},
987 saddles2[i],
988 std::array<SimplexId, 2>{
989 this->dg_.getCellGreaterVertex(Cell{2, saddles2[i]}, triangulation),
990 this->dg_.getCellGreaterVertex(
991 Cell{1, *boundSet.begin()}, triangulation),
992 }};
993 }
994 }
995
996 const auto nSadSadPairs = pairs.size() - nSadExtrPairs;
997
998 this->printMsg(
999 "Computed " + std::to_string(nSadSadPairs) + " saddle-saddle pairs", 1.0,
1000 tm2.getElapsedTime(), this->threadNumber_);
1001
1002 this->printMsg("saddle-saddle pairs sequential part", 1.0,
1003 tmseq.getElapsedTime(), 1, debug::LineMode::NEW,
1005}
1006
1007template <typename triangulationType>
1009 std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
1010 std::array<std::vector<SimplexId>, 4> &critCellsOrder,
1011 const SimplexId *const offsets,
1012 const triangulationType &triangulation,
1013 const bool sortEdges) const {
1014
1015 Timer tm{};
1016
1017 this->dg_.getCriticalPoints(criticalCellsByDim, triangulation);
1018
1019 this->printMsg("Extracted critical cells", 1.0, tm.getElapsedTime(),
1020 this->threadNumber_, debug::LineMode::NEW,
1022
1023 // memory allocations
1024 auto &critEdges{this->critEdges_};
1025 if(!sortEdges) {
1026 critEdges.resize(criticalCellsByDim[1].size());
1027 }
1028 std::vector<TriangleSimplex> critTriangles(criticalCellsByDim[2].size());
1029 std::vector<TetraSimplex> critTetras(criticalCellsByDim[3].size());
1030
1031#ifdef TTK_ENABLE_OPENMP
1032#pragma omp parallel num_threads(threadNumber_)
1033#endif // TTK_ENABLE_OPENMP
1034 {
1035 if(sortEdges) {
1036#ifdef TTK_ENABLE_OPENMP
1037#pragma omp for nowait
1038#endif // TTK_ENABLE_OPENMP
1039 for(size_t i = 0; i < critEdges.size(); ++i) {
1040 critEdges[i].fillEdge(i, offsets, triangulation);
1041 }
1042 } else {
1043#ifdef TTK_ENABLE_OPENMP
1044#pragma omp for nowait
1045#endif // TTK_ENABLE_OPENMP
1046 for(size_t i = 0; i < critEdges.size(); ++i) {
1047 critEdges[i].fillEdge(criticalCellsByDim[1][i], offsets, triangulation);
1048 }
1049 }
1050
1051#ifdef TTK_ENABLE_OPENMP
1052#pragma omp for nowait
1053#endif // TTK_ENABLE_OPENMP
1054 for(size_t i = 0; i < critTriangles.size(); ++i) {
1055 critTriangles[i].fillTriangle(
1056 criticalCellsByDim[2][i], offsets, triangulation);
1057 }
1058
1059#ifdef TTK_ENABLE_OPENMP
1060#pragma omp for
1061#endif // TTK_ENABLE_OPENMP
1062 for(size_t i = 0; i < critTetras.size(); ++i) {
1063 critTetras[i].fillTetra(criticalCellsByDim[3][i], offsets, triangulation);
1064 }
1065 }
1066
1067 TTK_PSORT(this->threadNumber_, critEdges.begin(), critEdges.end());
1068 TTK_PSORT(this->threadNumber_, critTriangles.begin(), critTriangles.end());
1069 TTK_PSORT(this->threadNumber_, critTetras.begin(), critTetras.end());
1070
1071#ifdef TTK_ENABLE_OPENMP
1072#pragma omp parallel num_threads(threadNumber_)
1073#endif // TTK_ENABLE_OPENMP
1074 {
1075#ifdef TTK_ENABLE_OPENMP
1076#pragma omp for nowait
1077#endif // TTK_ENABLE_OPENMP
1078 for(size_t i = 0; i < critEdges.size(); ++i) {
1079 critCellsOrder[1][critEdges[i].id_] = i;
1080 }
1081
1082#ifdef TTK_ENABLE_OPENMP
1083#pragma omp for nowait
1084#endif // TTK_ENABLE_OPENMP
1085 for(size_t i = 0; i < critTriangles.size(); ++i) {
1086 criticalCellsByDim[2][i] = critTriangles[i].id_;
1087 critCellsOrder[2][critTriangles[i].id_] = i;
1088 }
1089
1090#ifdef TTK_ENABLE_OPENMP
1091#pragma omp for
1092#endif // TTK_ENABLE_OPENMP
1093 for(size_t i = 0; i < critTetras.size(); ++i) {
1094 criticalCellsByDim[3][i] = critTetras[i].id_;
1095 critCellsOrder[3][critTetras[i].id_] = i;
1096 }
1097 }
1098
1099 if(sortEdges) {
1100 TTK_PSORT(this->threadNumber_, criticalCellsByDim[1].begin(),
1101 criticalCellsByDim[1].end(),
1102 [&critCellsOrder](const SimplexId a, const SimplexId b) {
1103 return critCellsOrder[1][a] < critCellsOrder[1][b];
1104 });
1105 } else {
1106#ifdef TTK_ENABLE_OPENMP
1107#pragma omp parallel for num_threads(threadNumber_)
1108#endif // TTK_ENABLE_OPENMP
1109 for(size_t i = 0; i < critEdges.size(); ++i) {
1110 criticalCellsByDim[1][i] = critEdges[i].id_;
1111 }
1112 }
1113
1114 this->printMsg("Extracted & sorted critical cells", 1.0, tm.getElapsedTime(),
1115 this->threadNumber_, debug::LineMode::NEW,
1117}
1118
1119template <typename triangulationType>
1121 std::vector<PersistencePair> &pairs,
1122 const SimplexId *const offsets,
1123 const triangulationType &triangulation,
1124 const bool ignoreBoundary,
1125 const bool compute2SaddlesChildren) {
1126
1127 // allocate memory
1128 this->alloc(triangulation);
1129
1130 Timer tm{};
1131 pairs.clear();
1132 const auto dim = this->dg_.getDimensionality();
1133 this->Compute2SaddlesChildren = compute2SaddlesChildren;
1134
1135 // get every critical cell sorted them by dimension
1136 std::array<std::vector<SimplexId>, 4> criticalCellsByDim{};
1137 // holds the critical cells order
1138 auto &critCellsOrder{this->critCellsOrder_};
1139
1140 this->extractCriticalCells(
1141 criticalCellsByDim, critCellsOrder, offsets, triangulation, dim == 3);
1142
1143 // if minima are paired
1144 auto &pairedMinima{this->pairedCritCells_[0]};
1145 // if 1-saddles are paired
1146 auto &paired1Saddles{this->pairedCritCells_[1]};
1147 // if 2-saddles are paired
1148 auto &paired2Saddles{this->pairedCritCells_[dim - 1]};
1149 // if maxima are paired
1150 auto &pairedMaxima{this->pairedCritCells_[dim]};
1151
1152 // connected components (global min/max pair)
1153 size_t nConnComp{};
1154
1155 if(this->ComputeMinSad) {
1156 // minima - saddle pairs
1157 this->getMinSaddlePairs(pairs, pairedMinima, paired1Saddles,
1158 criticalCellsByDim[1], critCellsOrder[1], offsets,
1159 triangulation);
1160
1161 // non-paired minima
1162 for(const auto min : criticalCellsByDim[0]) {
1163 if(!pairedMinima[min]) {
1164 pairs.emplace_back(min, -1, 0);
1165 pairedMinima[min] = true;
1166 nConnComp++;
1167 }
1168 }
1169 } else {
1170 // still extract the global pair
1171 const auto globMin{*std::min_element(
1172 criticalCellsByDim[0].begin(), criticalCellsByDim[0].end(),
1173 [offsets](const SimplexId a, const SimplexId b) {
1174 return offsets[a] < offsets[b];
1175 })};
1176 pairs.emplace_back(globMin, -1, 0);
1177 pairedMinima[globMin] = true;
1178 nConnComp++;
1179 }
1180
1181 if(dim > 1 && this->ComputeSadMax) {
1182 // saddle - maxima pairs
1183 this->getMaxSaddlePairs(
1184 pairs, pairedMaxima, paired2Saddles, criticalCellsByDim[dim - 1],
1185 critCellsOrder[dim - 1], critCellsOrder[dim], triangulation);
1186 }
1187
1188 if(ignoreBoundary) {
1189 // post-process saddle-max pairs: remove the one with the global
1190 // maximum (if it exists) to be (more) compatible with FTM
1191 const auto it
1192 = std::find_if(pairs.begin(), pairs.end(), [&](const PersistencePair &p) {
1193 if(p.type < dim - 1) {
1194 return false;
1195 }
1196 const Cell cmax{dim, p.death};
1197 const auto vmax{this->getCellGreaterVertex(cmax, triangulation)};
1198 return offsets[vmax] == triangulation.getNumberOfVertices() - 1;
1199 });
1200
1201 if(it != pairs.end()) {
1202 // remove saddle-max pair with global maximum
1203 paired2Saddles[it->birth] = false;
1204 pairedMaxima[it->death] = false;
1205 pairs.erase(it);
1206 }
1207 }
1208
1209 // saddle - saddle pairs
1210 if(dim == 3 && !criticalCellsByDim[1].empty()
1211 && !criticalCellsByDim[2].empty() && this->ComputeSadSad) {
1212 std::vector<GeneratorType> tmp{};
1213 this->getSaddleSaddlePairs(
1214 pairs, paired1Saddles, paired2Saddles, false, tmp, criticalCellsByDim[1],
1215 criticalCellsByDim[2], critCellsOrder[1], triangulation);
1216 }
1217
1218 if(std::is_same<triangulationType, ttk::ExplicitTriangulation>::value) {
1219 // create infinite pairs from non-paired 1-saddles, 2-saddles and maxima
1220 size_t nHandles{}, nCavities{}, nNonPairedMax{};
1221 if((dim == 2 && !ignoreBoundary && this->ComputeMinSad
1222 && this->ComputeSadMax)
1223 || (dim == 3 && this->ComputeMinSad && this->ComputeSadSad)) {
1224 // non-paired 1-saddles
1225 for(const auto s1 : criticalCellsByDim[1]) {
1226 if(!paired1Saddles[s1]) {
1227 paired1Saddles[s1] = true;
1228 // topological handles
1229 pairs.emplace_back(s1, -1, 1);
1230 nHandles++;
1231 }
1232 }
1233 }
1234 if(dim == 3 && !ignoreBoundary && this->ComputeSadMax
1235 && this->ComputeSadSad) {
1236 // non-paired 2-saddles
1237 for(const auto s2 : criticalCellsByDim[2]) {
1238 if(!paired2Saddles[s2]) {
1239 paired2Saddles[s2] = true;
1240 // cavities
1241 pairs.emplace_back(s2, -1, 2);
1242 nCavities++;
1243 }
1244 }
1245 }
1246 if(dim == 2 && !ignoreBoundary && this->ComputeSadMax) {
1247 // non-paired maxima
1248 for(const auto max : criticalCellsByDim[dim]) {
1249 if(!pairedMaxima[max]) {
1250 pairs.emplace_back(max, -1, 2);
1251 nNonPairedMax++;
1252 }
1253 }
1254 }
1255
1256 int nBoundComp
1257 = (dim == 3 ? nCavities : nHandles) + nConnComp - nNonPairedMax;
1258 nBoundComp = std::max(nBoundComp, 0);
1259
1260 // print Betti numbers
1261 std::vector<std::vector<std::string>> rows{
1262 {" #Connected components", std::to_string(nConnComp)},
1263 {" #Topological handles", std::to_string(nHandles)},
1264 {" #Cavities", std::to_string(nCavities)},
1265 {" #Boundary components", std::to_string(nBoundComp)},
1266 };
1267
1268 this->printMsg(rows, debug::Priority::DETAIL);
1269 }
1270
1271 this->printMsg(
1272 "Computed " + std::to_string(pairs.size()) + " persistence pairs", 1.0,
1273 tm.getElapsedTime(), this->threadNumber_);
1274
1275 this->displayStats(pairs, criticalCellsByDim, pairedMinima, paired1Saddles,
1276 paired2Saddles, pairedMaxima);
1277
1278 // free memory
1279 this->clear();
1280
1281 return 0;
1282}
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
Definition: OpenMP.h:46
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
int threadNumber_
Definition: BaseClass.h:95
virtual int setThreadNumber(const int threadNumber)
Definition: BaseClass.h:80
Minimalist debugging class.
Definition: Debug.h:88
int debugLevel_
Definition: Debug.h:379
virtual int setDebugLevel(const int &debugLevel)
Definition: Debug.cpp:147
int printMsg(const std::string &msg, const debug::Priority &priority=debug::Priority::INFO, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cout) const
Definition: Debug.h:118
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()
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)
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)
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 getDescendingPath(const Cell &cell, std::vector< Cell > &vpath, const triangulationType &triangulation) const
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
void fillEdge(const SimplexId id, const SimplexId *const offsets, const triangulationType &triangulation)
Type for exporting persistent generators.
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)
void fillTriangle(const SimplexId id, const SimplexId *const offsets, const triangulationType &triangulation)
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)