75 template <
class triangulationType = AbstractTriangulation>
77 const triangulationType *triangulation);
79 template <
class triangulationType = AbstractTriangulation>
81 const triangulationType *triangulation);
83 template <
class triangulationType = AbstractTriangulation>
85 const triangulationType *triangulation);
87 template <
class triangulationType>
90 template <
class triangulationType = AbstractTriangulation>
94 const triangulationType *triangulation,
95 bool &isLowerOnBoundary,
96 bool &isUpperOnBoundary,
97 std::vector<std::vector<ttk::SimplexId>> *upperComponents,
98 std::vector<std::vector<ttk::SimplexId>> *lowerComponents)
const;
100 template <
class triangulationType = AbstractTriangulation>
104 const triangulationType *triangulation,
105 std::vector<std::vector<ttk::SimplexId>> *upperComponents
107 std::vector<std::vector<ttk::SimplexId>> *lowerComponents
112 const std::vector<std::pair<SimplexId, SimplexId>>
113 &vertexLinkEdgeList)
const;
120 setOutput(std::vector<std::pair<SimplexId, char>> *criticalPoints) {
138 const std::vector<std::vector<std::pair<SimplexId, SimplexId>>>
160 const std::vector<std::vector<std::pair<SimplexId, SimplexId>>>
223 const SimplexId *
const offsets,
const triangulationType *triangulation) {
226#ifndef TTK_ENABLE_KAMIKAZE
227 if((!dimension_) && ((!triangulation) || (triangulation->isEmpty())))
229 if((!vertexNumber_) && ((!triangulation) || (triangulation->isEmpty())))
231 if((!vertexLinkEdgeLists_) && (!triangulation))
238 vertexNumber_ = triangulation->getNumberOfVertices();
239 dimension_ = triangulation->getCellVertexNumber(0) - 1;
242 printMsg(
"Extracting critical points...");
248#ifdef TTK_ENABLE_OPENMP
250 = std::max(1000, (
int)vertexNumber_ / (threadNumber_ * 100));
255#ifdef TTK_ENABLE_OPENMP
256#pragma omp parallel for schedule(dynamic, chunkSize) num_threads(threadNumber_)
260 if(!isRunningWithMPI()
261 || (isRunningWithMPI()
265 vertexTypes[i] = getCriticalType(i, offsets, triangulation);
268 }
else if(vertexLinkEdgeLists_) {
270#ifdef TTK_ENABLE_OPENMP
271#pragma omp parallel for schedule(dynamic, chunkSize) num_threads(threadNumber_)
275 vertexTypes[i] = getCriticalType(i, offsets, (*vertexLinkEdgeLists_)[i]);
279 SimplexId minimumNumber = 0, maximumNumber = 0, saddleNumber = 0,
280 oneSaddleNumber = 0, twoSaddleNumber = 0, monkeySaddleNumber = 0;
284 if(dimension_ == 3) {
285 for(
SimplexId i = 0; i < vertexNumber_; i++) {
286 switch(vertexTypes[i]) {
305 monkeySaddleNumber++;
309 }
else if(dimension_ == 2) {
310 for(
SimplexId i = 0; i < vertexNumber_; i++) {
311 switch(vertexTypes[i]) {
326 monkeySaddleNumber++;
333 std::vector<std::vector<std::string>> stats;
335 if(dimension_ == 3) {
336 stats.push_back({
" #1-saddles",
std::to_string(oneSaddleNumber)});
337 stats.push_back({
" #2-saddles",
std::to_string(twoSaddleNumber)});
339 if(dimension_ == 2) {
342 stats.push_back({
" #Multi-saddles",
std::to_string(monkeySaddleNumber)});
350 criticalPoints_->clear();
351 criticalPoints_->reserve(vertexNumber_);
352 for(
SimplexId i = 0; i < vertexNumber_; i++) {
354 criticalPoints_->emplace_back(i, vertexTypes[i]);
368 const triangulationType *triangulation,
369 bool &isLowerOnBoundary,
370 bool &isUpperOnBoundary,
371 std::vector<std::vector<ttk::SimplexId>> *upperComponents,
372 std::vector<std::vector<ttk::SimplexId>> *lowerComponents)
const {
375 = triangulation->getVertexNeighborNumber(vertexId);
376 std::vector<SimplexId> lowerNeighbors, upperNeighbors;
378 for(
SimplexId i = 0; i < neighborNumber; i++) {
380 triangulation->getVertexNeighbor(vertexId, i, neighborId);
382 if(offsets[neighborId] < offsets[vertexId]) {
383 lowerNeighbors.push_back(neighborId);
384 if(dimension_ == 3) {
385 if(triangulation->isVertexOnBoundary(neighborId))
386 isLowerOnBoundary =
true;
391 if(offsets[neighborId] > offsets[vertexId]) {
392 upperNeighbors.push_back(neighborId);
393 if(dimension_ == 3) {
394 if(triangulation->isVertexOnBoundary(neighborId))
395 isUpperOnBoundary =
true;
401 std::vector<UnionFind> lowerSeeds(lowerNeighbors.size());
402 std::vector<UnionFind *> lowerList(lowerNeighbors.size());
403 std::vector<UnionFind> upperSeeds(upperNeighbors.size());
404 std::vector<UnionFind *> upperList(upperNeighbors.size());
407 lowerList[i] = &(lowerSeeds[i]);
410 upperList[i] = &(upperSeeds[i]);
413 SimplexId const vertexStarSize = triangulation->getVertexStarNumber(vertexId);
415 for(
SimplexId i = 0; i < vertexStarSize; i++) {
417 triangulation->getVertexStar(vertexId, i, cellId);
419 SimplexId const cellSize = triangulation->getCellVertexNumber(cellId);
420 for(
SimplexId j = 0; j < cellSize; j++) {
422 triangulation->getCellVertex(cellId, j, neighborId0);
424 if(neighborId0 != vertexId) {
427 bool const lower0 = offsets[neighborId0] < offsets[vertexId];
430 for(
SimplexId k = j + 1; k < cellSize; k++) {
433 triangulation->getCellVertex(cellId, k, neighborId1);
435 if((neighborId1 != neighborId0) && (neighborId1 != vertexId)) {
437 bool const lower1 = offsets[neighborId1] < offsets[vertexId];
439 std::vector<SimplexId> *neighbors = &lowerNeighbors;
440 std::vector<UnionFind *> *seeds = &lowerList;
443 neighbors = &upperNeighbors;
447 if(lower0 == lower1) {
451 if((*neighbors)[l] == neighborId0) {
454 if((*neighbors)[l] == neighborId1) {
458 if((lowerId0 != -1) && (lowerId1 != -1)) {
460 (*seeds)[lowerId0], (*seeds)[lowerId1]);
461 (*seeds)[lowerId1] = (*seeds)[lowerId0];
474 lowerList[i] = lowerList[i]->find();
476 upperList[i] = upperList[i]->find();
478 std::unordered_map<UnionFind *, std::vector<ttk::SimplexId>>::iterator it;
479 std::unordered_map<UnionFind *, std::vector<ttk::SimplexId>>
481 std::unordered_map<UnionFind *, std::vector<ttk::SimplexId>>
486 it = upperComponentId.find(upperList[i]);
487 if(it != upperComponentId.end()) {
488 upperComponentId[upperList[i]].push_back(upperNeighbors[i]);
490 upperComponentId[upperList[i]]
491 = std::vector<ttk::SimplexId>(1, upperNeighbors[i]);
494 for(
auto elt : upperComponentId) {
495 upperComponents->push_back(std::vector<ttk::SimplexId>());
497 upperComponents->back().push_back(elt.second.at(i));
502 it = lowerComponentId.find(lowerList[i]);
503 if(it != lowerComponentId.end()) {
504 lowerComponentId[lowerList[i]].push_back(lowerNeighbors[i]);
506 lowerComponentId[lowerList[i]]
507 = std::vector<ttk::SimplexId>(1, lowerNeighbors[i]);
510 for(
auto elt : lowerComponentId) {
511 lowerComponents->push_back(std::vector<ttk::SimplexId>());
513 lowerComponents->back().push_back(elt.second.at(i));
531 const triangulationType *triangulation,
532 std::vector<std::vector<ttk::SimplexId>> *upperComponents,
533 std::vector<std::vector<ttk::SimplexId>> *lowerComponents)
const {
535 bool isLowerOnBoundary =
false, isUpperOnBoundary =
false;
536 std::vector<std::vector<ttk::SimplexId>> localUpperComponents;
537 std::vector<std::vector<ttk::SimplexId>> localLowerComponents;
538 if(upperComponents ==
nullptr) {
539 upperComponents = &localUpperComponents;
541 if(lowerComponents ==
nullptr) {
542 lowerComponents = &localLowerComponents;
544 getLowerUpperComponents(vertexId, offsets, triangulation, isLowerOnBoundary,
545 isUpperOnBoundary, upperComponents, lowerComponents);
546 ttk::SimplexId const lowerComponentNumber = lowerComponents->size();
547 ttk::SimplexId const upperComponentNumber = upperComponents->size();
549 if(dimension_ == 1) {
550 if(lowerComponentNumber == 0 && upperComponentNumber != 0) {
552 }
else if(lowerComponentNumber != 0 && upperComponentNumber == 0) {
554 }
else if(lowerComponentNumber == 1 && upperComponentNumber == 1) {
560 if(lowerComponentNumber == 0 && upperComponentNumber == 1) {
562 }
else if(lowerComponentNumber == 1 && upperComponentNumber == 0) {
564 }
else if(lowerComponentNumber == 1 && upperComponentNumber == 1) {
566 if((dimension_ == 3) && (triangulation->isVertexOnBoundary(vertexId))) {
568 if((isUpperOnBoundary) && (!isLowerOnBoundary))
570 if((!isUpperOnBoundary) && (isLowerOnBoundary))
578 if(dimension_ == 2 || dimension_ == 1) {
579 if((lowerComponentNumber == 2 && upperComponentNumber == 1)
580 || (lowerComponentNumber == 1 && upperComponentNumber == 2)
581 || (lowerComponentNumber == 2 && upperComponentNumber == 2)) {
592 }
else if(dimension_ == 3) {
593 if(lowerComponentNumber == 2 && upperComponentNumber == 1) {
595 }
else if(lowerComponentNumber == 1 && upperComponentNumber == 2) {