TTK
Loading...
Searching...
No Matches
LocalizedTopologicalSimplification.h
Go to the documentation of this file.
1
19
20#pragma once
21
22// for numerical perturbation
23#include <boost/math/special_functions/next.hpp>
24
25#include <Debug.h>
26#include <Triangulation.h>
27
28#include <Propagation.h>
29
30namespace ttk {
31 namespace lts {
32
33 inline std::string toFixed(const float &number, const int precision = 2) {
34 std::stringstream vFraction;
35 vFraction << std::fixed << std::setprecision(precision) << number;
36 return vFraction.str();
37 }
38
39 template <typename IT>
40 std::string
41 toFixed(const IT &number0, const IT &number1, const int precision = 2) {
42 return toFixed(((float)number0) / ((float)number1), precision);
43 }
44
46
47 public:
48 enum class PAIR_TYPE {
52 };
53
55 this->setDebugMsgPrefix("LTS");
56 }
58
60 ttk::AbstractTriangulation *triangulation) const {
61 return triangulation->preconditionVertexNeighbors();
62 }
63
66 template <typename IT>
67 int allocateMemory(std::vector<IT> &segmentation,
68 std::vector<IT> &queueMask,
69 std::vector<IT> &localOrder,
70 std::vector<Propagation<IT> *> &propagationMask,
71 std::vector<std::tuple<IT, IT, IT>> &sortedIndices,
72
73 const IT &nVertices) const {
74 ttk::Timer timer;
75
76 constexpr char msg[] = "Allocating Memory";
77 this->printMsg(
78 msg, 0, 0, this->threadNumber_, debug::LineMode::REPLACE);
79
80 segmentation.resize(nVertices);
81 queueMask.resize(nVertices);
82 localOrder.resize(nVertices);
83 propagationMask.resize(nVertices);
84 sortedIndices.resize(nVertices);
85
86 this->printMsg(msg, 1, timer.getElapsedTime(), this->threadNumber_);
87
88 return 0;
89 }
90
93 template <typename IT>
94 int initializeMemory(IT *segmentation,
95 IT *queueMask,
96 IT *localOrder,
97 Propagation<IT> **propagationMask,
98
99 const IT &nVertices) const {
100 ttk::Timer timer;
101
102 constexpr char msg[] = "Initializing Memory";
103 this->printMsg(
104 msg, 0, 0, this->threadNumber_, debug::LineMode::REPLACE);
105
106#ifdef TTK_ENABLE_OPENMP
107#pragma omp parallel for num_threads(this->threadNumber_)
108#endif // TTK_ENABLE_OPENMP
109 for(IT i = 0; i < nVertices; i++) {
110 segmentation[i] = -1;
111 queueMask[i] = -1;
112 localOrder[i] = -1;
113 propagationMask[i] = nullptr;
114 }
115
116 this->printMsg(msg, 1, timer.getElapsedTime(), this->threadNumber_);
117
118 return 0;
119 }
120
128 template <typename IT, class TT>
130 std::vector<Propagation<IT>> &propagations,
131 IT *authorizationMask, // assumes preservation mask is initialized as -1
132 IT *maximaBuffer,
133
134 const IT *authorizedExtremaIndices,
135 const IT &nAuthorizedExtremaIndices,
136 const IT *order,
137 const TT *triangulation) const {
138
139 ttk::Timer timer;
140 this->printMsg("Initializing Propagations", 0, 0, this->threadNumber_,
142
143 const IT nVertices = triangulation->getNumberOfVertices();
144#ifdef TTK_ENABLE_OPENMP
145#pragma omp parallel for num_threads( \
146 this->threadNumber_) if(nAuthorizedExtremaIndices > 1000)
147#endif // TTK_ENABLE_OPENMP
148 for(IT i = 0; i < nAuthorizedExtremaIndices; i++)
149 authorizationMask[authorizedExtremaIndices[i]] = -2;
150
151 IT writeIndex = 0;
152// find discarded maxima
153#ifdef TTK_ENABLE_OPENMP
154#pragma omp parallel for num_threads(this->threadNumber_)
155#endif // TTK_ENABLE_OPENMP
156 for(IT v = 0; v < nVertices; v++) {
157
158 // if v needs to be preserved then skip
159 if(authorizationMask[v] == -2)
160 continue;
161
162 // check if v has larger neighbors
163 bool hasLargerNeighbor = false;
164 const IT &vOrder = order[v];
165 const IT nNeighbors = triangulation->getVertexNeighborNumber(v);
166 for(IT n = 0; n < nNeighbors; n++) {
167 IT u{-1};
168 triangulation->getVertexNeighbor(v, n, u);
169 if(vOrder < order[u]) {
170 hasLargerNeighbor = true;
171 break;
172 }
173 }
174
175 // if v has larger neighbors then v can not be maximum
176 if(hasLargerNeighbor)
177 continue;
178
179 // get local write index for this thread
180 IT localWriteIndex = 0;
181#ifdef TTK_ENABLE_OPENMP
182#pragma omp atomic capture
183#endif // TTK_ENABLE_OPENMP
184 localWriteIndex = writeIndex++;
185
186// atomic write to prevent cache line conflicts
187#ifdef TTK_ENABLE_OPENMP
188#pragma omp atomic write
189#endif // TTK_ENABLE_OPENMP
190 maximaBuffer[localWriteIndex] = v;
191 }
192
193 // sort maxima in ascending order
194 std::sort(maximaBuffer, maximaBuffer + writeIndex,
195 [=](const IT &a, const IT &b) -> bool {
196 return order[a] < order[b];
197 });
198
199 // if there are no authorized maxima then always skip most persistent
200 // prop
201 if(nAuthorizedExtremaIndices < 1)
202 writeIndex--;
203
204 // init propagations
205 propagations.resize(writeIndex);
206#ifdef TTK_ENABLE_OPENMP
207#pragma omp parallel for num_threads(this->threadNumber_)
208#endif // TTK_ENABLE_OPENMP
209 for(IT p = 0; p < writeIndex; p++) {
210 propagations[p].criticalPoints.push_back(maximaBuffer[p]);
211 }
212
213 this->printMsg("Initializing Propagations ("
214 + std::to_string(writeIndex) + "|"
215 + std::to_string(nVertices) + ")",
216 1, timer.getElapsedTime(), this->threadNumber_);
217
218 return 0;
219 }
220
226 template <typename IT, typename TT>
228 Propagation<IT> &propagation,
229 Propagation<IT> **propagationMask,
230 IT *segmentation, // used here to also store registered larger vertices
231 IT *queueMask, // used to mark vertices that have already been added to
232 // the queue by this thread
233
234 const TT *triangulation,
235 const IT *order) const {
236
237 // pointer used to compare against representative
238 auto *currentPropagation = &propagation;
239
240 // add extremumIndex (stored in segmentId) to queue
241 IT segmentId = currentPropagation->criticalPoints[0];
242 auto queue = &currentPropagation->queue;
243 queue->emplace(order[segmentId], segmentId);
244
245 queueMask[segmentId] = segmentId;
246
247 // grow segment until prop reaches a saddle and then decide if prop
248 // should continue
249 IT v = -1;
250 while(!queue->empty()) {
251 v = std::get<1>(queue->top());
252 queue->pop();
253
254 // continue if this thread has already seen this vertex
255 if(propagationMask[v] != nullptr)
256 continue;
257
258 const IT &orderV = order[v];
259
260 // add neighbors to queue AND check if v is a saddle
261 bool isSaddle = false;
262 const IT nNeighbors = triangulation->getVertexNeighborNumber(v);
263
264 IT numberOfLargerNeighbors = 0;
265 IT numberOfLargerNeighborsThisThreadVisited = 0;
266 for(IT n = 0; n < nNeighbors; n++) {
267 IT u{-1};
268 triangulation->getVertexNeighbor(v, n, u);
269
270 const IT &orderU = order[u];
271
272 // if larger neighbor
273 if(orderU > orderV) {
274 numberOfLargerNeighbors++;
275
276 auto uPropagation = propagationMask[u];
277 if(uPropagation == nullptr
278 || currentPropagation != uPropagation->find())
279 isSaddle = true;
280 else
281 numberOfLargerNeighborsThisThreadVisited++;
282 } else if(queueMask[u] != segmentId) {
283 queue->emplace(orderU, u);
284 queueMask[u] = segmentId;
285 }
286 }
287
288 // if v is a saddle we have to check if the current thread is the last
289 // visitor
290 if(isSaddle) {
291 currentPropagation->criticalPoints.push_back(v);
292
293 IT numberOfRegisteredLargerVertices = 0;
294#ifdef TTK_ENABLE_OPENMP
295#pragma omp atomic capture
296#endif // TTK_ENABLE_OPENMP
297 {
298 segmentation[v] -= numberOfLargerNeighborsThisThreadVisited;
299 numberOfRegisteredLargerVertices = segmentation[v];
300 }
301
302 // if this thread did not register the last remaining larger
303 // vertices then terminate propagation
304 if(numberOfRegisteredLargerVertices != -numberOfLargerNeighbors - 1)
305 return 0;
306
307 // get most dominant propagation
308 std::vector<Propagation<IT> *> neighborPropagations(
309 nNeighbors, nullptr);
310 for(IT n = 0; n < nNeighbors; n++) {
311 IT u{-1};
312 triangulation->getVertexNeighbor(v, n, u);
313 if(propagationMask[u] != nullptr) {
314 auto &neighborPropagation = neighborPropagations[n];
315 neighborPropagation = propagationMask[u]->find();
316 if(order[neighborPropagation->criticalPoints[0]]
317 > order[currentPropagation->criticalPoints[0]])
318 currentPropagation = neighborPropagation;
319 }
320 }
321
322 // merge propagations
323 for(IT n = 0; n < nNeighbors; n++) {
325 currentPropagation, neighborPropagations[n]);
326 }
327
328 queue = &currentPropagation->queue;
329 segmentId = currentPropagation->criticalPoints[0];
330 }
331
332 // mark vertex as visited and continue
333 propagationMask[v] = currentPropagation;
334 segmentation[v] = segmentId;
335 currentPropagation->segmentSize++;
336 }
337
338 this->printErr(
339 "Simple propagations should never reach global minimum/maximum.");
340
341 return 1;
342 }
343
348 template <typename IT, typename DT, typename TT>
350 Propagation<IT> &propagation,
351 Propagation<IT> **propagationMask,
352 IT *segmentation, // used here to also store registered larger vertices
353 IT *queueMask, // used to mark vertices that have already been added to
354 // the queue by this thread
355
356 const TT *triangulation,
357 const IT *order,
358 const DT *scalars,
359 const DT persistenceThreshold) const {
360
361 // pointer used to compare against representative
362 auto *currentPropagation = &propagation;
363
364 // add extremumIndex (stored in segmentId) to queue
365 IT segmentId = currentPropagation->criticalPoints[0];
366 auto queue = &currentPropagation->queue;
367 queue->emplace(order[segmentId], segmentId);
368
369 DT s0 = scalars[segmentId];
370
371 queueMask[segmentId] = segmentId;
372
373 // grow segment until prop reaches a saddle and then decide if prop
374 // should continue
375 IT v = -1;
376 while(!queue->empty()) {
377 v = std::get<1>(queue->top());
378 queue->pop();
379
380 // continue if this thread has already seen this vertex
381 if(propagationMask[v] != nullptr)
382 continue;
383
384 // if propagation exceeds persistence threshold stop
385 const DT &s1 = scalars[v];
386 const DT &sd = s0 < s1 ? s1 - s0 : s0 - s1;
387 if(sd > persistenceThreshold) {
388 currentPropagation->aborted = true;
389 return 0;
390 }
391
392 const IT &orderV = order[v];
393
394 // add neighbors to queue AND check if v is a saddle
395 bool isSaddle = false;
396 const IT nNeighbors = triangulation->getVertexNeighborNumber(v);
397
398 IT numberOfLargerNeighbors = 0;
399 IT numberOfLargerNeighborsThisThreadVisited = 0;
400 for(IT n = 0; n < nNeighbors; n++) {
401 IT u{-1};
402 triangulation->getVertexNeighbor(v, n, u);
403
404 const IT &orderU = order[u];
405
406 // if larger neighbor
407 if(orderU > orderV) {
408 numberOfLargerNeighbors++;
409
410 auto uPropagation = propagationMask[u];
411 if(uPropagation == nullptr
412 || currentPropagation != uPropagation->find())
413 isSaddle = true;
414 else
415 numberOfLargerNeighborsThisThreadVisited++;
416 } else if(queueMask[u] != segmentId) {
417 queue->emplace(orderU, u);
418 queueMask[u] = segmentId;
419 }
420 }
421
422 // if v is a saddle we have to check if the current thread is the last
423 // visitor
424 if(isSaddle) {
425 currentPropagation->criticalPoints.push_back(v);
426
427 IT numberOfRegisteredLargerVertices = 0;
428#ifdef TTK_ENABLE_OPENMP
429#pragma omp atomic capture
430#endif // TTK_ENABLE_OPENMP
431 {
432 segmentation[v] -= numberOfLargerNeighborsThisThreadVisited;
433 numberOfRegisteredLargerVertices = segmentation[v];
434 }
435
436 // if this thread did not register the last remaining larger
437 // vertices then terminate propagation
438 if(numberOfRegisteredLargerVertices != -numberOfLargerNeighbors - 1)
439 return 0;
440
441 // get most dominant propagation
442 std::vector<Propagation<IT> *> neighborPropagations(
443 nNeighbors, nullptr);
444 for(IT n = 0; n < nNeighbors; n++) {
445 IT u{-1};
446 triangulation->getVertexNeighbor(v, n, u);
447 if(propagationMask[u] != nullptr) {
448 auto &neighborPropagation = neighborPropagations[n];
449 neighborPropagation = propagationMask[u]->find();
450 if(order[neighborPropagation->criticalPoints[0]]
451 > order[currentPropagation->criticalPoints[0]])
452 currentPropagation = neighborPropagation;
453 }
454 }
455
456 // merge propagations
457 for(IT n = 0; n < nNeighbors; n++) {
459 currentPropagation, neighborPropagations[n]);
460 }
461
462 queue = &currentPropagation->queue;
463 segmentId = currentPropagation->criticalPoints[0];
464 s0 = scalars[segmentId];
465 }
466
467 // mark vertex as visited and continue
468 propagationMask[v] = currentPropagation;
469 segmentation[v] = segmentId;
470 currentPropagation->segmentSize++;
471 }
472
473 return 1;
474 }
475
478 template <typename IT, class TT>
479 int computeSimplePropagations(std::vector<Propagation<IT>> &propagations,
480 Propagation<IT> **propagationMask,
481 IT *segmentation,
482 IT *queueMask,
483
484 const TT *triangulation,
485 const IT *inputOrder) const {
486 ttk::Timer timer;
487 const IT nPropagations = propagations.size();
488 const std::string msg
489 = "Computing Propagations (" + std::to_string(nPropagations) + ")";
490 this->printMsg(
491 msg, 0, 0, this->threadNumber_, debug::LineMode::REPLACE);
492
493 int status = 0;
494// compute propagations
495#ifdef TTK_ENABLE_OPENMP
496#pragma omp parallel for schedule(dynamic, 1) num_threads(this->threadNumber_)
497#endif // TTK_ENABLE_OPENMP
498 for(IT p = 0; p < nPropagations; p++) {
499 int localStatus = this->computeSimplePropagation<IT, TT>(
500 propagations[p], propagationMask, segmentation, queueMask,
501
502 triangulation, inputOrder);
503
504 if(localStatus)
505 status = 1;
506 }
507
508 if(status)
509 return 1;
510
511 this->printMsg(msg, 1, timer.getElapsedTime(), this->threadNumber_);
512
513 return 0;
514 }
515
518 template <typename IT, typename DT, class TT>
520 std::vector<Propagation<IT>> &propagations,
521 Propagation<IT> **propagationMask,
522 IT *segmentation,
523 IT *queueMask,
524
525 const TT *triangulation,
526 const IT *order,
527 const DT *scalars,
528 const DT persistenceThreshold) const {
529 ttk::Timer timer;
530 const IT nPropagations = propagations.size();
531 const std::string msg = "Computing Persistence-Sensitive Propagations ("
532 + std::to_string(nPropagations) + ")";
533 this->printMsg(
534 msg, 0, 0, this->threadNumber_, debug::LineMode::REPLACE);
535
536 int status = 0;
537// compute propagations
538#ifdef TTK_ENABLE_OPENMP
539#pragma omp parallel for schedule(dynamic, 1) num_threads(this->threadNumber_)
540#endif // TTK_ENABLE_OPENMP
541 for(IT p = 0; p < nPropagations; p++) {
542 int localStatus
543 = this->computePersistenceSensitivePropagation<IT, DT, TT>(
544 propagations[p], propagationMask, segmentation, queueMask,
545
546 triangulation, order, scalars, persistenceThreshold);
547
548 if(localStatus)
549 status = 1;
550 }
551
552 if(status)
553 return 1;
554
555 this->printMsg(msg, 1, timer.getElapsedTime(), this->threadNumber_);
556
557 return 0;
558 }
559
565 template <typename IT>
566 int
567 finalizePropagations(std::vector<Propagation<IT> *> &parentPropagations,
568 std::vector<Propagation<IT>> &propagations,
569
570 const IT nVertices) const {
571 ttk::Timer timer;
572
573 const IT nPropagations = propagations.size();
574
575 this->printMsg(
576 "Finalizing Propagations (" + std::to_string(nPropagations) + ")", 0,
577 timer.getElapsedTime(), this->threadNumber_,
579
580 IT nSegmentVertices = 0;
581 IT nParentPropagations = 0;
582 parentPropagations.resize(nPropagations);
583 for(IT p = 0; p < nPropagations; p++) {
584 auto *propagation = &propagations[p];
585 if(!propagation->aborted
586 && (propagation->parent == propagation
587 || propagation->parent->aborted)) {
588 nSegmentVertices = nSegmentVertices + propagation->segmentSize;
589 parentPropagations[nParentPropagations++] = propagation;
590 }
591 }
592 parentPropagations.resize(nParentPropagations);
593
594 this->printMsg("Finalizing Propagations ("
595 + std::to_string(nParentPropagations) + "|"
596 + toFixed(nParentPropagations, nPropagations) + "|"
597 + toFixed(nSegmentVertices, nVertices) + ")",
598 1, timer.getElapsedTime(), this->threadNumber_);
599
600 return 0;
601 }
602
606 template <typename IT, class TT>
607 int computeSegment(IT *segmentation,
608 Propagation<IT> *propagation,
609
610 const IT *order,
611 const TT *triangulation) const {
612
613 const IT &extremumIndex = propagation->criticalPoints[0];
614 const IT &saddleIndex = propagation->criticalPoints.back();
615 const IT &saddleOrder = order[saddleIndex];
616
617 // collect segment
618 auto &segment = propagation->segment;
619
620 segment.resize(propagation->segmentSize);
621 IT segmentIndex = 0;
622 if(propagation->segmentSize > 0) {
623 std::vector<IT> queue(propagation->segmentSize);
624 IT queueIndex = 0;
625
626 // init queue
627 {
628 queue[queueIndex++] = extremumIndex;
629 segmentation[extremumIndex] = -1000; // mark as visited
630 }
631
632 // flood fill by starting from extremum
633 while(queueIndex > 0) {
634 const IT v = queue[--queueIndex];
635
636 segment[segmentIndex++] = v;
637
638 IT nNeighbors = triangulation->getVertexNeighborNumber(v);
639 for(IT n = 0; n < nNeighbors; n++) {
640 IT u{-1};
641 triangulation->getVertexNeighbor(v, n, u);
642 auto &s = segmentation[u];
643 if(s >= 0 && order[u] > saddleOrder) {
644 s = -1000; // mark as visited
645 queue[queueIndex++] = u; // add to queue
646 }
647 }
648 }
649 }
650
651 if(segmentIndex != propagation->segmentSize) {
652 this->printErr("Segment size incorrect: "
653 + std::to_string(segmentIndex) + " "
654 + std::to_string(propagation->segmentSize));
655 return 1;
656 }
657
658 for(auto idx : propagation->segment)
659 segmentation[idx] = extremumIndex;
660
661 return 0;
662 }
663
665 template <typename IT, class TT>
666 int computeSegments(IT *segmentation,
667 std::vector<Propagation<IT> *> &propagations,
668
669 const IT *order,
670 const TT *triangulation) const {
671
672 const IT nPropagations = propagations.size();
673 const IT nVertices = triangulation->getNumberOfVertices();
674
675 ttk::Timer timer;
676 const std::string msg
677 = "Computing Segments (" + std::to_string(nPropagations) + ")";
678 this->printMsg(
679 msg, 0, 0, this->threadNumber_, debug::LineMode::REPLACE);
680
681 int status = 0;
682
683// compute segments in parallel
684#ifdef TTK_ENABLE_OPENMP
685#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
686#endif // TTK_ENABLE_OPENMP
687 for(IT p = 0; p < nPropagations; p++) {
688 int localStatus
689 = this->computeSegment<IT, TT>(segmentation, propagations[p],
690
691 order, triangulation);
692 if(localStatus)
693 status = 1;
694 }
695 if(status)
696 return 1;
697
698 // print status
699 if(this->debugLevel_ < 4 || nPropagations == 0) {
700 this->printMsg(msg, 1, timer.getElapsedTime(), this->threadNumber_);
701 } else {
702
703 IT min = propagations[0]->segmentSize;
704 IT max = min;
705 IT avg = 0;
706
707 for(IT p = 0; p < nPropagations; p++) {
708 const auto propagation = propagations[p];
709 if(min > propagation->segmentSize)
710 min = propagation->segmentSize;
711 if(max < propagation->segmentSize)
712 max = propagation->segmentSize;
713 avg += propagation->segmentSize;
714 }
715
716 avg /= nPropagations;
717
718 this->printMsg("Computing Segments (" + std::to_string(nPropagations)
719 + "|" + toFixed(min, nVertices) + "|"
720 + toFixed(avg, nVertices) + "|"
721 + toFixed(max, nVertices) + ")",
722 1, timer.getElapsedTime(), this->threadNumber_);
723 }
724
725 return 0;
726 }
727
728 template <typename IT, class TT>
730 IT *localOrder,
731 IT *localVertexSequence,
732
733 const bool &performSuperlevelSetPropagation,
734 const TT *triangulation,
735 const IT *segmentation,
736 const IT &segmentId,
737 const std::vector<IT> &boundary,
738 const std::vector<IT> &segment,
739 const IT &saddleIdx) const {
740 // init priority queue
741 std::priority_queue<std::pair<IT, IT>, std::vector<std::pair<IT, IT>>>
742 queue;
743
744 IT nSegmentVertices = segment.size();
745
746 if(performSuperlevelSetPropagation) {
747 // add saddle to queue
748 queue.emplace(std::numeric_limits<IT>::max(), saddleIdx);
749 } else {
750 // invert order
751 IT offset = -nSegmentVertices - 1; // ensures that inverted order < 0
752 for(IT i = 0; i < nSegmentVertices; i++)
753 localOrder[segment[i]] = offset - localOrder[segment[i]];
754
755 // add all boundary vertices to queue
756 for(const auto &v : boundary) {
757 queue.emplace(localOrder[v], v);
758 localOrder[v] = 0;
759 }
760
761 // also add saddle to queue
762 queue.emplace(std::numeric_limits<IT>::min(), saddleIdx);
763 }
764
765 IT q = 0;
766 while(!queue.empty()) {
767 IT v = std::get<1>(queue.top());
768 queue.pop();
769
770 localVertexSequence[q++] = v;
771
772 IT nNeighbors = triangulation->getVertexNeighborNumber(v);
773 for(IT n = 0; n < nNeighbors; n++) {
774 IT u{-1};
775 triangulation->getVertexNeighbor(v, n, u);
776
777 // if u is in segment and has not already been added to the queue
778 if(segmentation[u] == segmentId && localOrder[u] < 0) {
779 queue.emplace(localOrder[u], u);
780 localOrder[u] = 0;
781 }
782 }
783 }
784
785 if(performSuperlevelSetPropagation) {
786 // skip first idx as it corresponds to saddle
787 for(IT i = 1; i <= nSegmentVertices; i++)
788 localOrder[localVertexSequence[i]] = -i;
789 } else {
790 IT order = -nSegmentVertices;
791 // skip last idx as it corresponds to saddle
792 for(IT i = 0; i < nSegmentVertices; i++)
793 localOrder[localVertexSequence[i]] = order++;
794 }
795
796 return 0;
797 }
798
799 template <typename IT, class TT>
800 int computeLocalOrderOfSegment(IT *localOrder,
801
802 const Propagation<IT> *propagation,
803 const TT *triangulation,
804 const IT *segmentation,
805 const IT *inputOrder) const {
806 // quick escape for small segments
807 if(propagation->segmentSize == 1) {
808 localOrder[propagation->segment[0]] = -2;
809 return 0;
810 }
811
812 // init local order by input order
813 // Note: force orders to be smaller than 0 for in situ updates in
814 // propagation procedures
815 IT nVertices = triangulation->getNumberOfVertices();
816 for(const auto &v : propagation->segment)
817 localOrder[v] = inputOrder[v] - nVertices;
818
819 const IT &extremumIndex = propagation->criticalPoints[0];
820 const IT &saddleIndex = propagation->criticalPoints.back();
821
822 // vector to record boundary
823 std::vector<IT> boundary;
824
825 // make enough room for segment + saddle
826 std::vector<IT> localVertexSequence(propagation->segmentSize + 1);
827
828 int status = 0;
829 bool containsResidualExtrema = true;
830 bool performSuperlevelSetPropagation = true;
831 while(containsResidualExtrema) {
832 propagation->nIterations++;
833
834 if(this->debugLevel_ > 3 && propagation->nIterations == 20)
835 this->printWrn("Unable to converge with less than 20 iterations!");
836
837 // execute superlevel set propagation
838 status = this->computeLocalOrderOfSegmentIteration<IT, TT>(
839 localOrder, localVertexSequence.data(),
840
841 performSuperlevelSetPropagation, triangulation, segmentation,
842 extremumIndex, boundary, propagation->segment, saddleIndex);
843 if(status)
844 return 1;
845
846 performSuperlevelSetPropagation = !performSuperlevelSetPropagation;
847
848 IT boundaryWriteIdx = 0;
849 IT nResidualMaxima = 0;
850 IT nResidualMinima = 0;
851
852 for(const auto &v : propagation->segment) {
853 bool isOnSegmentBoundary = false;
854 bool hasSmallerNeighbor = false;
855 bool hasLargerNeighbor = false;
856
857 const auto &vOrder = localOrder[v];
858
859 IT nNeighbors = triangulation->getVertexNeighborNumber(v);
860 for(IT n = 0; n < nNeighbors; n++) {
861 IT u{-1};
862 triangulation->getVertexNeighbor(v, n, u);
863
864 if(u == saddleIndex) {
865 hasLargerNeighbor = true;
866 continue;
867 }
868
869 // if u is not inside segment -> v is on segment boundary
870 if(segmentation[u] != extremumIndex) {
871 isOnSegmentBoundary = true;
872 } else {
873 if(localOrder[u] > vOrder)
874 hasLargerNeighbor = true;
875 else
876 hasSmallerNeighbor = true;
877 }
878 }
879
880 if(!hasLargerNeighbor)
881 nResidualMaxima++;
882
883 if(isOnSegmentBoundary) {
884 localVertexSequence[boundaryWriteIdx++] = v;
885 } else if(!hasSmallerNeighbor) {
886 nResidualMinima++;
887 }
888 }
889
890 containsResidualExtrema = nResidualMinima > 0 || nResidualMaxima > 0;
891
892 if(containsResidualExtrema && boundary.size() == 0) {
893 boundary.resize(boundaryWriteIdx);
894 for(IT i = 0; i < boundaryWriteIdx; i++) {
895 boundary[i] = localVertexSequence[i];
896 }
897 }
898 }
899
900 return 0;
901 }
902
903 template <typename IT, class TT>
905 IT *localOrder,
906
907 const TT *triangulation,
908 const IT *segmentation,
909 const IT *inputOrder,
910 const std::vector<Propagation<IT> *> &propagations) const {
911 ttk::Timer timer;
912
913 const IT nPropagations = propagations.size();
914 this->printMsg("Computing Local Order of Segments ("
915 + std::to_string(nPropagations) + ")",
917
918 int status = 0;
919#ifdef TTK_ENABLE_OPENMP
920#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
921#endif // TTK_ENABLE_OPENMP
922 for(IT p = 0; p < nPropagations; p++) {
923 int localStatus = this->computeLocalOrderOfSegment<IT>(
924 localOrder,
925
926 propagations[p], triangulation, segmentation, inputOrder);
927 if(localStatus)
928 status = 1;
929 }
930 if(status)
931 return 1;
932
933// enforce that saddles have the highest local order
934#ifdef TTK_ENABLE_OPENMP
935#pragma omp parallel for num_threads(this->threadNumber_)
936#endif // TTK_ENABLE_OPENMP
937 for(IT p = 0; p < nPropagations; p++) {
938 localOrder[propagations[p]->criticalPoints.back()]
939 = std::numeric_limits<IT>::max();
940 }
941
942 if(this->debugLevel_ < 4 || nPropagations == 0) {
943 this->printMsg("Computing Local Order of Segments ("
944 + std::to_string(nPropagations) + ")",
945 1, timer.getElapsedTime(), this->threadNumber_);
946 } else {
947 IT min = propagations[0]->nIterations;
948 IT max = min;
949 IT avg = 0;
950
951 for(IT p = 0; p < nPropagations; p++) {
952 const auto propagation = propagations[p];
953 if(min > propagation->nIterations)
954 min = propagation->nIterations;
955 if(max < propagation->nIterations)
956 max = propagation->nIterations;
957 avg += propagation->nIterations;
958 }
959
960 avg /= nPropagations;
961
962 this->printMsg("Computing Local Order of Segments ("
963 + std::to_string(nPropagations) + "|"
964 + std::to_string(min) + "|" + std::to_string(avg)
965 + "|" + std::to_string(max) + ")",
966 1, timer.getElapsedTime(), this->threadNumber_);
967 }
968
969 return 0;
970 }
971
972 template <typename IT>
974 IT *outputOrder,
975 const std::vector<Propagation<IT> *> &parentPropagations) const {
976
977 ttk::Timer timer;
978 this->printMsg("Flattening Order Array", 0, 0, this->threadNumber_,
980
981 const IT nParentPropagations = parentPropagations.size();
982
983// flatten segemnt to order of last encountered saddles
984#ifdef TTK_ENABLE_OPENMP
985#pragma omp parallel for num_threads(this->threadNumber_)
986#endif // TTK_ENABLE_OPENMP
987 for(IT p = 0; p < nParentPropagations; p++) {
988 const auto *propagation = parentPropagations[p];
989 const auto &saddleOrder
990 = outputOrder[propagation->criticalPoints.back()];
991 for(const auto &v : propagation->segment)
992 outputOrder[v] = saddleOrder;
993 }
994
995 this->printMsg("Flattening Order Array", 1, timer.getElapsedTime(),
996 this->threadNumber_);
997
998 return 0;
999 }
1000
1001 template <typename DT, typename IT>
1002 int flattenScalars(DT *scalars,
1003 const std::vector<Propagation<IT>> &propagationsA,
1004 const std::vector<Propagation<IT>> &propagationsB
1005 = {}) const {
1006
1007 ttk::Timer timer;
1008 this->printMsg("Flattening Scalar Array", 0, 0, this->threadNumber_,
1010
1011 std::vector<const std::vector<Propagation<IT>> *> propagationsPair
1012 = {&propagationsA, &propagationsB};
1013
1014 for(const auto propagations : propagationsPair) {
1015 const IT nPropagations = propagations->size();
1016
1017// flatten scalars to saddle value
1018#ifdef TTK_ENABLE_OPENMP
1019#pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_)
1020#endif // TTK_ENABLE_OPENMP
1021 for(IT p = 0; p < nPropagations; p++) {
1022 const auto propagation = &((*propagations)[p]);
1023
1024 if(!propagation->aborted
1025 && (propagation->parent == propagation
1026 || propagation->parent->aborted)) {
1027 const auto &saddleScalar
1028 = scalars[propagation->criticalPoints.back()];
1029 for(const auto &v : propagation->segment)
1030 scalars[v] = saddleScalar;
1031 }
1032 }
1033 }
1034
1035 this->printMsg("Flattening Scalar Array", 1, timer.getElapsedTime(),
1036 this->threadNumber_);
1037
1038 return 0;
1039 }
1040
1041 template <typename IT>
1043 IT *order,
1044 const IT *localOrder,
1045 std::vector<std::tuple<IT, IT, IT>> &sortedIndices) const {
1046 ttk::Timer timer;
1047
1048 const IT nVertices = sortedIndices.size();
1049
1050// init tuples
1051#ifdef TTK_ENABLE_OPENMP
1052#pragma omp parallel for num_threads(this->threadNumber_)
1053#endif // TTK_ENABLE_OPENMP
1054 for(IT i = 0; i < nVertices; i++) {
1055 auto &t = sortedIndices[i];
1056 std::get<0>(t) = order[i];
1057 std::get<1>(t) = localOrder[i];
1058 std::get<2>(t) = i;
1059 }
1060
1061 this->printMsg("Computing Global Order", 0.2, timer.getElapsedTime(),
1062 this->threadNumber_, debug::LineMode::REPLACE);
1063
1064 TTK_PSORT(
1065 this->threadNumber_, sortedIndices.begin(), sortedIndices.end());
1066
1067 this->printMsg("Computing Global Order", 0.8, timer.getElapsedTime(),
1068 this->threadNumber_, debug::LineMode::REPLACE);
1069
1070#ifdef TTK_ENABLE_OPENMP
1071#pragma omp parallel for num_threads(this->threadNumber_)
1072#endif // TTK_ENABLE_OPENMP
1073 for(IT i = 0; i < nVertices; i++)
1074 order[std::get<2>(sortedIndices[i])] = i;
1075
1076 this->printMsg("Computing Global Order", 1, timer.getElapsedTime(),
1077 this->threadNumber_);
1078
1079 return 0;
1080 }
1081
1082 template <typename DT, typename IT>
1084 DT *scalars,
1085 const std::vector<std::tuple<IT, IT, IT>> &sortedIndices,
1086 const bool descending = false) const {
1087 ttk::Timer timer;
1088 this->printMsg("Applying numerical perturbation", 0, 0,
1090
1091 const IT nVertices = sortedIndices.size();
1092
1093 if(descending)
1094 for(IT i = 1; i < nVertices; i++) {
1095 const IT &v0 = std::get<2>(sortedIndices[i - 1]);
1096 const IT &v1 = std::get<2>(sortedIndices[i]);
1097 if(scalars[v0] >= scalars[v1])
1098 scalars[v1] = boost::math::float_next(scalars[v0]);
1099 }
1100 else
1101 for(IT i = nVertices - 1; i > 0; i--) {
1102 const IT &v0 = std::get<2>(sortedIndices[i]);
1103 const IT &v1 = std::get<2>(sortedIndices[i - 1]);
1104 if(scalars[v0] >= scalars[v1])
1105 scalars[v1] = boost::math::float_next(scalars[v0]);
1106 }
1107
1108 this->printMsg("Applying numerical perturbation", 1,
1109 timer.getElapsedTime(), this->threadNumber_);
1110
1111 return 0;
1112 }
1113
1114 template <typename IT, class TT>
1116 IT *order,
1117 IT *segmentation,
1118 IT *queueMask,
1119 IT *localOrder,
1120 Propagation<IT> **propagationMask,
1121 std::vector<Propagation<IT>> &propagations,
1122 std::vector<std::tuple<IT, IT, IT>> &sortedIndices,
1123
1124 const TT *triangulation,
1125 const IT *authorizedExtremaIndices,
1126 const IT &nAuthorizedExtremaIndices) const {
1127 const IT nVertices = triangulation->getNumberOfVertices();
1128
1129 int status = 0;
1130
1131 // init propagations
1132 status = this->initializeMemory<IT>(segmentation, queueMask, localOrder,
1133 propagationMask,
1134
1135 nVertices);
1136 if(status)
1137 return 1;
1138
1139 // init propagations
1140 status = this->initializePropagations<IT, TT>(
1141 propagations,
1142 queueMask, // use as authorization mask (will be overridden by
1143 // subsequent procedures)
1144 localOrder, // use as maxima buffer (will be overridden by subsequent
1145 // procedures)
1146
1147 authorizedExtremaIndices, nAuthorizedExtremaIndices, order,
1148 triangulation);
1149 if(status)
1150 return 1;
1151
1152 // compute propagations
1153 status = this->computeSimplePropagations<IT, TT>(
1154 propagations, propagationMask, segmentation, queueMask,
1155
1156 triangulation, order);
1157 if(status)
1158 return 1;
1159
1160 // finalize master propagations
1161 std::vector<Propagation<IT> *> parentPropagations;
1162 status
1163 = this->finalizePropagations<IT>(parentPropagations, propagations,
1164
1165 nVertices);
1166 if(status)
1167 return 1;
1168
1169 // compute segments
1170 status = this->computeSegments<IT, TT>(segmentation, parentPropagations,
1171
1172 order, triangulation);
1173 if(status)
1174 return 1;
1175
1176 // compute local order of segments
1177 status = this->computeLocalOrderOfSegments<IT, TT>(
1178 localOrder,
1179
1180 triangulation, segmentation, order, parentPropagations);
1181 if(status)
1182 return 1;
1183
1184 // flatten order
1185 status = this->flattenOrder<IT>(order, parentPropagations);
1186 if(status)
1187 return 1;
1188
1189 // compute global offsets
1190 status = this->computeGlobalOrder<IT>(order, localOrder, sortedIndices);
1191 if(status)
1192 return 1;
1193
1194 return 0;
1195 }
1196
1197 template <typename IT, typename DT, class TT>
1199 DT *scalars,
1200 IT *order,
1201 IT *segmentation,
1202 IT *queueMask,
1203 IT *localOrder,
1204 Propagation<IT> **propagationMask,
1205 std::vector<Propagation<IT>> &propagations,
1206 std::vector<std::tuple<IT, IT, IT>> &sortedIndices,
1207
1208 const TT *triangulation,
1209 const DT persistenceThreshold) const {
1210
1211 const IT nVertices = triangulation->getNumberOfVertices();
1212
1213 int status = 0;
1214
1215 // init propagations
1216 status = this->initializeMemory<IT>(segmentation, queueMask, localOrder,
1217 propagationMask,
1218
1219 nVertices);
1220 if(status)
1221 return 1;
1222
1223 // init propagations
1224 status = this->initializePropagations<IT, TT>(
1225 propagations,
1226 queueMask, // will be ignored since there are no authorized maxima
1227 localOrder, // use as maxima buffer (will be overridden by subsequent
1228 // procedures)
1229
1230 nullptr, 0, order, triangulation);
1231 if(status)
1232 return 1;
1233
1234 // compute propagations
1235 status = this->computePersistenceSensitivePropagations<IT, DT, TT>(
1236 propagations, propagationMask, segmentation, queueMask,
1237
1238 triangulation, order, scalars, persistenceThreshold);
1239 if(status)
1240 return 1;
1241
1242 // finalize master propagations
1243 std::vector<Propagation<IT> *> parentPropagations;
1244 status
1245 = this->finalizePropagations<IT>(parentPropagations, propagations,
1246
1247 nVertices);
1248 if(status)
1249 return 1;
1250
1251 // compute segments
1252 status = this->computeSegments<IT, TT>(segmentation, parentPropagations,
1253
1254 order, triangulation);
1255 if(status)
1256 return 1;
1257
1258 // compute local order of segments
1259 status = this->computeLocalOrderOfSegments<IT, TT>(
1260 localOrder,
1261
1262 triangulation, segmentation, order, parentPropagations);
1263 if(status)
1264 return 1;
1265
1266 // flatten order
1267 status = this->flattenOrder<IT>(order, parentPropagations);
1268 if(status)
1269 return 1;
1270
1271 // compute global offsets
1272 status = this->computeGlobalOrder<IT>(order, localOrder, sortedIndices);
1273 if(status)
1274 return 1;
1275
1276 status = this->flattenScalars<DT, IT>(scalars, propagations);
1277 if(status)
1278 return 1;
1279
1280 return 0;
1281 }
1282
1283 template <typename IT>
1284 int invertOrder(IT *outputOrder, const IT &nVertices) const {
1285 ttk::Timer timer;
1286 this->printMsg("Inverting Order", 0, 0, this->threadNumber_,
1288
1289 const auto nVerticesM1 = nVertices - 1;
1290
1291#ifdef TTK_ENABLE_OPENMP
1292#pragma omp parallel for num_threads(this->threadNumber_)
1293#endif // TTK_ENABLE_OPENMP
1294 for(IT v = 0; v < nVertices; v++)
1295 outputOrder[v] = nVerticesM1 - outputOrder[v];
1296
1297 this->printMsg(
1298 "Inverting Order", 1, timer.getElapsedTime(), this->threadNumber_);
1299
1300 return 0;
1301 }
1302
1303 template <typename DT, typename IT, class TT>
1305 IT *order,
1306
1307 const TT *triangulation,
1308 const IT *authorizedExtremaIndices,
1309 const IT &nAuthorizedExtremaIndices,
1310 const bool &computePerturbation) const {
1311 ttk::Timer globalTimer;
1312
1313 IT nVertices = triangulation->getNumberOfVertices();
1314
1315 // Allocating Memory
1316 int status = 0;
1317 std::vector<IT> segmentation;
1318 std::vector<IT> queueMask;
1319 std::vector<IT> localOrder;
1320 std::vector<Propagation<IT> *> propagationMask;
1321 std::vector<Propagation<IT>> propagationsMax;
1322 std::vector<Propagation<IT>> propagationsMin;
1323 std::vector<std::tuple<IT, IT, IT>> sortedIndices;
1324
1325 this->allocateMemory<IT>(segmentation, queueMask, localOrder,
1326 propagationMask, sortedIndices,
1327
1328 nVertices);
1329
1330 bool hasMinima{false};
1331 bool hasMaxima{false};
1332 for(IT i = 0; i < nAuthorizedExtremaIndices; i++) {
1333 const auto &extremum{authorizedExtremaIndices[i]};
1334 const auto nNeigh{triangulation->getVertexNeighborNumber(extremum)};
1335 if(nNeigh > 0) {
1336 // look at the first neighbor to determine if minimum or maximum
1337 SimplexId neigh{};
1338 triangulation->getVertexNeighbor(extremum, 0, neigh);
1339 if(order[extremum] > order[neigh]) {
1340 hasMaxima = true;
1341 } else {
1342 hasMinima = true;
1343 }
1344 }
1345 // don't look further if we have both minima and maxima
1346 if(hasMaxima && hasMinima) {
1347 break;
1348 }
1349 }
1350
1351 // Maxima
1352 if(hasMaxima) {
1353 this->printMsg("----------- [Removing Unauthorized Maxima]",
1355
1356 status = this->detectAndRemoveUnauthorizedMaxima<IT, TT>(
1357 order, segmentation.data(), queueMask.data(), localOrder.data(),
1358 propagationMask.data(), propagationsMax, sortedIndices,
1359
1360 triangulation, authorizedExtremaIndices, nAuthorizedExtremaIndices);
1361 if(status)
1362 return 1;
1363 }
1364
1365 // Minima
1366 if(hasMinima) {
1367 this->printMsg("----------- [Removing Unauthorized Minima]",
1369
1370 // invert order
1371 if(this->invertOrder(order, nVertices))
1372 return 1;
1373
1374 status = this->detectAndRemoveUnauthorizedMaxima<IT, TT>(
1375 order, segmentation.data(), queueMask.data(), localOrder.data(),
1376 propagationMask.data(), propagationsMin, sortedIndices,
1377
1378 triangulation, authorizedExtremaIndices, nAuthorizedExtremaIndices);
1379 if(status)
1380 return 1;
1381
1382 // revert order
1383 if(this->invertOrder(order, nVertices))
1384 return 1;
1385 }
1386
1387 // flatten scalars
1388 status = this->flattenScalars<DT, IT>(
1389 scalars, propagationsMax, propagationsMin);
1390 if(status)
1391 return 1;
1392
1393 // optionally compute perturbation
1394 if(computePerturbation) {
1396 status = this->computeNumericalPerturbation<DT, IT>(
1397 scalars, sortedIndices);
1398 if(status)
1399 return 1;
1400 }
1401
1403 this->printMsg(
1404 "Complete", 1, globalTimer.getElapsedTime(), this->threadNumber_);
1405
1407
1408 return 0;
1409 }
1410
1411 template <typename DT, typename IT, class TT>
1413 IT *order,
1414
1415 const TT *triangulation,
1416 const DT persistenceThreshold,
1417 const bool &computePerturbation,
1418 const PAIR_TYPE &pairType
1420 ttk::Timer globalTimer;
1421
1422 IT nVertices = triangulation->getNumberOfVertices();
1423
1424 // Allocating Memory
1425 int status = 0;
1426 std::vector<IT> segmentation;
1427 std::vector<IT> queueMask;
1428 std::vector<IT> localOrder;
1429 std::vector<Propagation<IT> *> propagationMask;
1430 std::vector<Propagation<IT>> propagationsMax;
1431 std::vector<Propagation<IT>> propagationsMin;
1432 std::vector<std::tuple<IT, IT, IT>> sortedIndices;
1433
1434 this->allocateMemory<IT>(segmentation, queueMask, localOrder,
1435 propagationMask, sortedIndices,
1436
1437 nVertices);
1438
1439 // Maxima
1440 if(pairType == PAIR_TYPE::EXTREMUM_SADDLE
1441 || pairType == PAIR_TYPE::MAXIMUM_SADDLE) {
1442 this->printMsg("----------- [Removing Non-Persistent Maxima]",
1444
1445 status = this->detectAndRemoveNonPersistentMaxima<IT, DT, TT>(
1446 scalars, order, segmentation.data(), queueMask.data(),
1447 localOrder.data(), propagationMask.data(), propagationsMax,
1448 sortedIndices,
1449
1450 triangulation, persistenceThreshold);
1451 if(status)
1452 return 1;
1453 }
1454
1455 // Minima
1456 if(pairType == PAIR_TYPE::EXTREMUM_SADDLE
1457 || pairType == PAIR_TYPE::MINIMUM_SADDLE) {
1458 this->printMsg("----------- [Removing Non-Persistent Minima]",
1460
1461 if(this->invertOrder(order, nVertices))
1462 return 1;
1463
1464 status = this->detectAndRemoveNonPersistentMaxima<IT, DT, TT>(
1465 scalars, order, segmentation.data(), queueMask.data(),
1466 localOrder.data(), propagationMask.data(), propagationsMin,
1467 sortedIndices,
1468
1469 triangulation, persistenceThreshold);
1470 if(status)
1471 return 1;
1472
1473 if(this->invertOrder(order, nVertices))
1474 return 1;
1475 }
1476
1477 // optionally compute perturbation
1478 if(computePerturbation) {
1480 status = this->computeNumericalPerturbation<DT, IT>(
1481 scalars, sortedIndices, pairType == PAIR_TYPE::MAXIMUM_SADDLE);
1482 if(status)
1483 return 1;
1484 }
1485
1487 this->printMsg(
1488 "Complete", 1, globalTimer.getElapsedTime(), this->threadNumber_);
1489
1491
1492 return 0;
1493 }
1494
1495 }; // class
1496 } // namespace lts
1497} // namespace ttk
#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
Minimalist debugging class.
Definition: Debug.h:88
int debugLevel_
Definition: Debug.h:379
int printWrn(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:159
void setDebugMsgPrefix(const std::string &prefix)
Definition: Debug.h:364
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
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
Definition: Debug.h:149
double getElapsedTime()
Definition: Timer.h:15
int flattenScalars(DT *scalars, const std::vector< Propagation< IT > > &propagationsA, const std::vector< Propagation< IT > > &propagationsB={}) const
int computeLocalOrderOfSegment(IT *localOrder, const Propagation< IT > *propagation, const TT *triangulation, const IT *segmentation, const IT *inputOrder) const
int computeSimplePropagations(std::vector< Propagation< IT > > &propagations, Propagation< IT > **propagationMask, IT *segmentation, IT *queueMask, const TT *triangulation, const IT *inputOrder) const
int computeSimplePropagation(Propagation< IT > &propagation, Propagation< IT > **propagationMask, IT *segmentation, IT *queueMask, const TT *triangulation, const IT *order) const
int computeSegments(IT *segmentation, std::vector< Propagation< IT > * > &propagations, const IT *order, const TT *triangulation) const
This method computes the segments of a given list of propagations.
int detectAndRemoveNonPersistentMaxima(DT *scalars, IT *order, IT *segmentation, IT *queueMask, IT *localOrder, Propagation< IT > **propagationMask, std::vector< Propagation< IT > > &propagations, std::vector< std::tuple< IT, IT, IT > > &sortedIndices, const TT *triangulation, const DT persistenceThreshold) const
int computeNumericalPerturbation(DT *scalars, const std::vector< std::tuple< IT, IT, IT > > &sortedIndices, const bool descending=false) const
int initializePropagations(std::vector< Propagation< IT > > &propagations, IT *authorizationMask, IT *maximaBuffer, const IT *authorizedExtremaIndices, const IT &nAuthorizedExtremaIndices, const IT *order, const TT *triangulation) const
int computeLocalOrderOfSegments(IT *localOrder, const TT *triangulation, const IT *segmentation, const IT *inputOrder, const std::vector< Propagation< IT > * > &propagations) const
int initializeMemory(IT *segmentation, IT *queueMask, IT *localOrder, Propagation< IT > **propagationMask, const IT &nVertices) const
int computePersistenceSensitivePropagations(std::vector< Propagation< IT > > &propagations, Propagation< IT > **propagationMask, IT *segmentation, IT *queueMask, const TT *triangulation, const IT *order, const DT *scalars, const DT persistenceThreshold) const
int finalizePropagations(std::vector< Propagation< IT > * > &parentPropagations, std::vector< Propagation< IT > > &propagations, const IT nVertices) const
int computeSegment(IT *segmentation, Propagation< IT > *propagation, const IT *order, const TT *triangulation) const
int allocateMemory(std::vector< IT > &segmentation, std::vector< IT > &queueMask, std::vector< IT > &localOrder, std::vector< Propagation< IT > * > &propagationMask, std::vector< std::tuple< IT, IT, IT > > &sortedIndices, const IT &nVertices) const
~LocalizedTopologicalSimplification() override=default
int computeGlobalOrder(IT *order, const IT *localOrder, std::vector< std::tuple< IT, IT, IT > > &sortedIndices) const
int flattenOrder(IT *outputOrder, const std::vector< Propagation< IT > * > &parentPropagations) const
int computePersistenceSensitivePropagation(Propagation< IT > &propagation, Propagation< IT > **propagationMask, IT *segmentation, IT *queueMask, const TT *triangulation, const IT *order, const DT *scalars, const DT persistenceThreshold) const
int preconditionTriangulation(ttk::AbstractTriangulation *triangulation) const
int detectAndRemoveUnauthorizedMaxima(IT *order, IT *segmentation, IT *queueMask, IT *localOrder, Propagation< IT > **propagationMask, std::vector< Propagation< IT > > &propagations, std::vector< std::tuple< IT, IT, IT > > &sortedIndices, const TT *triangulation, const IT *authorizedExtremaIndices, const IT &nAuthorizedExtremaIndices) const
int removeUnauthorizedExtrema(DT *scalars, IT *order, const TT *triangulation, const IT *authorizedExtremaIndices, const IT &nAuthorizedExtremaIndices, const bool &computePerturbation) const
int invertOrder(IT *outputOrder, const IT &nVertices) const
int computeLocalOrderOfSegmentIteration(IT *localOrder, IT *localVertexSequence, const bool &performSuperlevelSetPropagation, const TT *triangulation, const IT *segmentation, const IT &segmentId, const std::vector< IT > &boundary, const std::vector< IT > &segment, const IT &saddleIdx) const
int removeNonPersistentExtrema(DT *scalars, IT *order, const TT *triangulation, const DT persistenceThreshold, const bool &computePerturbation, const PAIR_TYPE &pairType=PAIR_TYPE::EXTREMUM_SADDLE) const
std::string toFixed(const float &number, const int precision=2)
The Topology ToolKit.
int SimplexId
Identifier type for simplices of any dimension.
Definition: DataTypes.h:22
Superlevel Set Component Propagation.
Definition: Propagation.h:11
std::vector< IT > criticalPoints
Definition: Propagation.h:20
std::vector< IT > segment
Definition: Propagation.h:23
Propagation * find()
Definition: Propagation.h:29
static int unify(Propagation< IT > *p0, Propagation< IT > *p1)
Definition: Propagation.h:42