93 std::vector<std::array<float, 3>>
points_{};
164 template <
typename dataType,
typename triangulationType>
169 const dataType *
const scalars,
170 const size_t scalarsMTime,
172 const triangulationType &triangulation);
186 const bool doDescending,
187 const bool doSaddleConnectors) {
197 const bool doDescending) {
206 const bool doDescending,
207 const bool doMorseSmale) {
264 template <
typename triangulationType>
267 std::vector<Separatrix> &separatrices,
268 const triangulationType &triangulation)
const;
274 template <
typename triangulationType>
276 std::vector<Separatrix> &separatrices,
277 const triangulationType &triangulation)
const;
283 template <
typename triangulationType>
285 std::vector<Separatrix> &separatrices,
286 const triangulationType &triangulation)
const;
291 template <
typename triangulationType>
293 const std::vector<Separatrix> &separatrices,
295 const triangulationType &triangulation)
const;
301 template <
typename triangulationType>
303 const std::vector<SimplexId> &saddles2,
304 std::vector<Separatrix> &separatrices,
305 std::vector<std::vector<SimplexId>> &separatricesSaddles,
306 const triangulationType &triangulation)
const;
312 template <
typename triangulationType>
315 const std::vector<Separatrix> &separatrices,
316 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
318 const triangulationType &triangulation)
const;
325 template <
typename triangulationType>
328 const size_t polSize,
329 const triangulationType &triangulation)
const;
334 template <
typename triangulationType>
336 const size_t polSize,
337 const triangulationType &triangulation)
const;
343 template <
typename triangulationType>
345 const std::vector<SimplexId> &saddles1,
346 std::vector<Separatrix> &separatrices,
347 std::vector<std::vector<SimplexId>> &separatricesSaddles,
348 const triangulationType &triangulation)
const;
354 template <
typename triangulationType>
357 const std::vector<Separatrix> &separatrices,
358 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
360 const triangulationType &triangulation)
const;
366 std::vector<std::vector<Separatrix>> &separatrices)
const;
371 template <
typename triangulationType>
374 const triangulationType &triangulation)
const;
379 template <
typename triangulationType>
382 const triangulationType &triangulation)
const;
388 template <
typename triangulationType>
390 const SimplexId *
const ascendingManifold,
391 const SimplexId *
const descendingManifold,
393 const triangulationType &triangulation)
const;
395 template <
typename dataType,
typename triangulationType>
397 const dataType *
const scalars,
399 const triangulationType &triangulation);
423template <
typename dataType,
typename triangulationType>
428 const dataType *
const scalars,
429 const size_t scalarsMTime,
431 const triangulationType &triangulation) {
432#ifndef TTK_ENABLE_KAMIKAZE
433 if(scalars ==
nullptr) {
434 this->
printErr(
"Input scalar field pointer is null.");
438 if(offsets ==
nullptr) {
439 this->
printErr(
"Input offset field pointer is null.");
448 const auto dim = triangulation.getDimensionality();
460 const auto nVerts{triangulation.getNumberOfVertices()};
463 const auto pair{std::minmax_element(offsets, offsets + nVerts)};
465 const auto globmin = std::distance(offsets, pair.first);
466 const auto globmax = std::distance(offsets, pair.second);
467 persistenceThreshold *= (scalars[globmax] - scalars[globmin]);
468 this->
printMsg(
"Absolute saddle connectors persistence threshold is "
469 + std::to_string(persistenceThreshold),
474 persistenceThreshold, scalars, offsets, triangulation);
477 std::array<std::vector<SimplexId>, 4> criticalPoints{};
481 this->
printMsg(
" Critical points extracted", 1.0, tm.getElapsedTime(),
486 std::vector<std::vector<Separatrix>> separatrices1{};
493 separatrices1.emplace_back();
496 criticalPoints[1], separatrices1.back(), triangulation);
498 this->
printMsg(
" Descending 1-separatrices computed", 1.0,
505 separatrices1.emplace_back();
508 criticalPoints[dim - 1], separatrices1.back(), triangulation);
510 this->
printMsg(
" Ascending 1-separatrices computed", 1.0,
518 separatrices1.emplace_back();
535 this->
printMsg(
" 1-separatrices set", 1.0, tmp.getElapsedTime(),
539 this->
printMsg(
"1-separatrices computed", 1.0, tm1sep.getElapsedTime(),
540 this->threadNumber_);
548 std::vector<Separatrix> separatrices;
549 std::vector<std::vector<SimplexId>> separatricesSaddles;
551 criticalPoints[2], separatrices, separatricesSaddles, triangulation);
553 outSeps2, separatrices, separatricesSaddles, offsets, triangulation);
555 this->
printMsg(
" Descending 2-separatrices computed", 1.0,
562 std::vector<Separatrix> separatrices;
563 std::vector<std::vector<SimplexId>> separatricesSaddles;
565 criticalPoints[1], separatrices, separatricesSaddles, triangulation);
567 outSeps2, separatrices, separatricesSaddles, offsets, triangulation);
569 this->
printMsg(
" Ascending 2-separatrices computed", 1.0,
576 this->
printMsg(
"2-separatrices computed", 1.0, tm2sep.getElapsedTime(),
577 this->threadNumber_);
585 criticalPoints[dim], outManifold.
ascending_, triangulation);
589 criticalPoints[0], outManifold.
descending_, triangulation);
599 "Segmentation computed", 1.0, tmp.
getElapsedTime(), this->threadNumber_);
615 + std::to_string(triangulation.getNumberOfVertices())
616 +
" points) processed",
626template <
typename triangulationType>
628 const std::vector<SimplexId> &saddles,
629 std::vector<Separatrix> &separatrices,
630 const triangulationType &triangulation)
const {
632 const SimplexId numberOfSaddles = saddles.size();
635 const SimplexId numberOfSeparatrices = 2 * numberOfSaddles;
636 separatrices.resize(numberOfSeparatrices);
639#ifdef TTK_ENABLE_OPENMP
640#pragma omp parallel for num_threads(threadNumber_)
642 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
643 const Cell saddle{1, saddles[i]};
647 const Cell &saddle1 = saddle;
649 for(
int j = 0; j < 2; ++j) {
651 triangulation.getEdgeVertex(saddle1.
id_, j, vertexId);
653 std::vector<Cell> vpath;
654 vpath.push_back(saddle1);
655 discreteGradient_.getDescendingPath(
656 Cell(0, vertexId), vpath, triangulation);
658 const Cell &lastCell = vpath.back();
659 if(lastCell.
dim_ == 0 and discreteGradient_.isCellCritical(lastCell)) {
660 separatrices[2 * i + j].source_ = saddle;
661 separatrices[2 * i + j].destination_ = lastCell;
662 separatrices[2 * i + j].geometry_ = std::move(vpath);
671template <
typename triangulationType>
673 const std::vector<SimplexId> &saddles,
674 std::vector<Separatrix> &separatrices,
675 const triangulationType &triangulation)
const {
677 const auto dim{triangulation.getDimensionality()};
680 auto getFaceStarNumber = &triangulationType::getTriangleStarNumber;
681 auto getFaceStar = &triangulationType::getTriangleStar;
684 getFaceStarNumber = &triangulationType::getEdgeStarNumber;
685 getFaceStar = &triangulationType::getEdgeStar;
688 const SimplexId numberOfSaddles = saddles.size();
690 std::vector<std::vector<Separatrix>> sepsPerSaddle(numberOfSaddles);
693#ifdef TTK_ENABLE_OPENMP
694#pragma omp parallel for num_threads(threadNumber_)
696 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
697 const Cell saddle{dim - 1, saddles[i]};
700 const auto starNumber{(triangulation.*getFaceStarNumber)(saddle.id_)};
701 for(
SimplexId j = 0; j < starNumber; ++j) {
704 (triangulation.*getFaceStar)(saddle.id_, j, sId);
706 std::vector<Cell> vpath{saddle};
707 discreteGradient_.getAscendingPath(
Cell(dim, sId), vpath, triangulation);
709 const Cell &lastCell = vpath.back();
710 if(lastCell.
dim_ == dim and discreteGradient_.isCellCritical(lastCell)) {
711 sepsPerSaddle[i].emplace_back();
712 sepsPerSaddle[i].back().source_ = saddle;
713 sepsPerSaddle[i].back().destination_ = lastCell;
714 sepsPerSaddle[i].back().geometry_ = std::move(vpath);
719 this->flattenSeparatricesVectors(sepsPerSaddle);
721 separatrices = std::move(sepsPerSaddle[0]);
726template <
typename triangulationType>
728 const std::vector<SimplexId> &saddles2,
729 std::vector<Separatrix> &separatrices,
730 const triangulationType &triangulation)
const {
732 const auto nTriangles = triangulation.getNumberOfTriangles();
734 std::vector<bool> isVisited(nTriangles,
false);
735 std::vector<SimplexId> visitedTriangles{};
737 using Vpath = std::vector<Cell>;
739 const auto dim{triangulation.getDimensionality()};
741 std::vector<std::vector<Separatrix>> sepsByThread(saddles2.size());
742 std::vector<SimplexId> saddles1{};
744#ifdef TTK_ENABLE_OPENMP
745#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
746 firstprivate(isVisited, visitedTriangles, saddles1)
748 for(
size_t i = 0; i < saddles2.size(); ++i) {
749 const Cell s2{dim - 1, saddles2[i]};
752 discreteGradient_.getDescendingWall(
753 s2, mask, triangulation,
nullptr, &saddles1);
755 for(
const auto saddle1Id : saddles1) {
756 const Cell s1{1, saddle1Id};
759 const bool isMultiConnected
760 = discreteGradient_.getAscendingPathThroughWall(
761 s1, s2, isVisited, &vpath, triangulation);
767 const auto &last = vpath.back();
769 if(!isMultiConnected && last.dim_ == s2.dim_ && last.id_ == s2.id_) {
770 sepsByThread[i].emplace_back();
771 sepsByThread[i].back().source_ = s1;
772 sepsByThread[i].back().destination_ = s2;
773 sepsByThread[i].back().geometry_ = std::move(vpath);
778 this->flattenSeparatricesVectors(sepsByThread);
780 separatrices = std::move(sepsByThread[0]);
785template <
typename triangulationType>
788 const std::vector<Separatrix> &separatrices,
790 const triangulationType &triangulation)
const {
808 std::vector<size_t> geomPointsBegId{npoints};
810 std::vector<size_t> geomCellsBegId{ncells};
813 for(
size_t i = 0; i < separatrices.size(); ++i) {
814 const auto &sep = separatrices[i];
815 const auto sepSize = sep.geometry_.size();
817 ncells += sepSize - 1;
818 geomPointsBegId.emplace_back(npoints);
819 geomCellsBegId.emplace_back(ncells);
822 const int dimensionality = triangulation.getCellVertexNumber(0) - 1;
836 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
837 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
840#ifdef TTK_ENABLE_OPENMP
841#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic)
843 for(
size_t i = 0; i < separatrices.size(); ++i) {
844 const auto &sep = separatrices[i];
845 const auto &sepGeom = sep.geometry_;
846 const auto sepId = separatrixId + i;
853 const auto saddleConnector
854 = dimensionality == 3 && src.
dim_ == 1 && dst.
dim_ == 2;
856 = saddleConnector ? 1 : std::min(dst.
dim_, dimensionality - 1);
860 return offsets[a] < offsets[b];
862 const std::array<SimplexId, 2> gVerts{
863 discreteGradient_.getCellGreaterVertex(src, triangulation),
864 discreteGradient_.getCellGreaterVertex(dst, triangulation)};
865 const auto sepFuncMax
866 = *std::max_element(gVerts.begin(), gVerts.end(), vertsOrder);
867 const std::array<SimplexId, 2> lVerts{
868 discreteGradient_.getCellLowerVertex(src, triangulation),
869 discreteGradient_.getCellLowerVertex(dst, triangulation)};
870 const auto sepFuncMin
871 = *std::min_element(lVerts.begin(), lVerts.end(), vertsOrder);
872 separatrixFunctionMaxima[sepId] = sepFuncMax;
873 separatrixFunctionMinima[sepId] = sepFuncMin;
876 const auto onBoundary
877 =
static_cast<char>(discreteGradient_.isBoundary(src, triangulation))
878 +
static_cast<char>(discreteGradient_.isBoundary(dst, triangulation));
880 for(
size_t j = 0; j < sepGeom.size(); ++j) {
881 const auto &cell = sepGeom[j];
882 std::array<float, 3> pt{};
883 triangulation.getCellIncenter(cell.id_, cell.dim_, pt.data());
886 const auto k = geomPointsBegId[i] + j;
888 points[3 * k + 0] = pt[0];
889 points[3 * k + 1] = pt[1];
890 points[3 * k + 2] = pt[2];
893 = (j == 0 || j == sepGeom.size() - 1) ? 0 : 1;
902 const auto l = geomCellsBegId[i] + j - 1;
904 cellsConn[2 * l + 0] = k - 1;
905 cellsConn[2 * l + 1] = k;
926template <
typename triangulationType>
928 const std::vector<SimplexId> &saddles1,
929 std::vector<Separatrix> &separatrices,
930 std::vector<std::vector<SimplexId>> &separatricesSaddles,
931 const triangulationType &triangulation)
const {
932 const Cell emptyCell;
934 const SimplexId numberOfSaddles = saddles1.size();
938 const SimplexId numberOfSeparatrices = numberOfSaddles;
939 separatrices.resize(numberOfSeparatrices);
940 separatricesSaddles.resize(numberOfSeparatrices);
942 const auto nEdges = triangulation.getNumberOfEdges();
943 std::vector<bool> isVisited(nEdges,
false);
944 std::vector<SimplexId> visitedEdges{};
947#ifdef TTK_ENABLE_OPENMP
948#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
949 firstprivate(isVisited, visitedEdges)
951 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
952 const Cell saddle1{1, saddles1[i]};
954 std::vector<Cell> wall;
956 discreteGradient_.getAscendingWall(
957 saddle1, mask, triangulation, &wall, &separatricesSaddles[i]);
959 separatrices[i].source_ = saddle1;
960 separatrices[i].destination_ = emptyCell;
961 separatrices[i].geometry_ = std::move(wall);
967template <
typename triangulationType>
969 const std::vector<SimplexId> &saddles2,
970 std::vector<Separatrix> &separatrices,
971 std::vector<std::vector<SimplexId>> &separatricesSaddles,
972 const triangulationType &triangulation)
const {
973 const Cell emptyCell;
975 const SimplexId numberOfSaddles = saddles2.size();
979 const SimplexId numberOfSeparatrices = numberOfSaddles;
980 separatrices.resize(numberOfSeparatrices);
981 separatricesSaddles.resize(numberOfSeparatrices);
983 const auto nTriangles = triangulation.getNumberOfTriangles();
984 std::vector<bool> isVisited(nTriangles,
false);
985 std::vector<SimplexId> visitedTriangles{};
987 const auto dim{triangulation.getDimensionality()};
990#ifdef TTK_ENABLE_OPENMP
991#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
992 firstprivate(isVisited, visitedTriangles)
994 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
995 const Cell saddle2{dim - 1, saddles2[i]};
997 std::vector<Cell> wall;
999 discreteGradient_.getDescendingWall(
1000 saddle2, mask, triangulation, &wall, &separatricesSaddles[i]);
1002 separatrices[i].source_ = saddle2;
1003 separatrices[i].destination_ = emptyCell;
1004 separatrices[i].geometry_ = std::move(wall);
1010template <
typename triangulationType>
1014 const size_t polSize,
1015 const triangulationType &triangulation)
const {
1017 for(
size_t i = 0; i < polSize; ++i) {
1019 triangulation.getEdgeStar(edgeId, i, starId);
1020 polygon[i] = starId;
1026template <
typename triangulationType>
1029 const size_t polSize,
1030 const triangulationType &triangulation)
const {
1032 for(
size_t i = 1; i < polSize; ++i) {
1035 bool isFound =
false;
1037 for(; j < polSize; ++j) {
1040 k < triangulation.getCellNeighborNumber(polygon[i - 1]); ++k) {
1042 triangulation.getCellNeighbor(polygon[i - 1], k, neighborId);
1043 if(neighborId == polygon[j]) {
1054 std::swap(polygon[j], polygon[i]);
1061template <
typename triangulationType>
1064 const std::vector<Separatrix> &separatrices,
1065 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
1067 const triangulationType &triangulation)
const {
1085 const auto noldcells{ncells};
1089 std::vector<size_t> geomCellsBegId{ncells};
1092 for(
size_t i = 0; i < separatrices.size(); ++i) {
1093 const auto &sep = separatrices[i];
1094 ncells += sep.geometry_.size();
1095 geomCellsBegId.emplace_back(ncells);
1099 std::vector<SimplexId> sepSourceIds(separatrices.size());
1100 std::vector<SimplexId> sepIds(separatrices.size());
1101 std::vector<char> sepOnBoundary(separatrices.size());
1102 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
1103 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
1105 std::vector<SimplexId> polygonNTetras(ncells - noldcells);
1106 std::vector<SimplexId> polygonEdgeIds(ncells - noldcells);
1107 std::vector<SimplexId> polygonSepInfosIds(ncells - noldcells);
1109#ifdef TTK_ENABLE_OPENMP
1110#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic)
1112 for(
size_t i = 0; i < separatrices.size(); ++i) {
1113 const auto &sep = separatrices[i];
1114 const auto &sepGeom = sep.geometry_;
1115 const auto &sepSaddles = separatricesSaddles[i];
1116 const auto sepId = separatrixId + i;
1120 const auto sepFuncMin
1121 = discreteGradient_.getCellLowerVertex(src, triangulation);
1123 if(!sepSaddles.empty()) {
1125 const auto maxId = *std::max_element(
1126 sepSaddles.begin(), sepSaddles.end(),
1128 return offsets[discreteGradient_.getCellGreaterVertex(
1129 Cell{2, a}, triangulation)]
1130 < offsets[discreteGradient_.getCellGreaterVertex(
1131 Cell{2, b}, triangulation)];
1134 = discreteGradient_.getCellGreaterVertex(
Cell{2, maxId}, triangulation);
1138 const auto maxId = *std::max_element(
1139 sepGeom.begin(), sepGeom.end(),
1140 [&triangulation, offsets,
this](
const Cell &a,
const Cell &b) {
1141 return offsets[discreteGradient_.getCellGreaterVertex(
1143 < offsets[discreteGradient_.getCellGreaterVertex(
1146 sepFuncMax = discreteGradient_.getCellGreaterVertex(maxId, triangulation);
1150 const char onBoundary
1151 = (sepSaddles.empty()
1153 : std::count_if(sepSaddles.begin(), sepSaddles.end(),
1155 return triangulation.isTriangleOnBoundary(a);
1157 + triangulation.isEdgeOnBoundary(src.id_);
1160 sepSourceIds[i] = src.id_;
1161 separatrixFunctionMaxima[sepId] = sepFuncMax;
1162 separatrixFunctionMinima[sepId] = sepFuncMin;
1163 sepOnBoundary[i] = onBoundary;
1165 for(
size_t j = 0; j < sepGeom.size(); ++j) {
1166 const auto &cell = sepGeom[j];
1168 const auto k = geomCellsBegId[i] + j - noldcells;
1170 polygonNTetras[k] = triangulation.getEdgeStarNumber(cell.id_);
1172 if(polygonNTetras[k] > 2) {
1173 polygonEdgeIds[k] = cell.id_;
1174 polygonSepInfosIds[k] = i;
1180 std::vector<SimplexId> validTetraIds{};
1181 validTetraIds.reserve(polygonNTetras.size());
1183 for(
size_t i = 0; i < polygonNTetras.size(); ++i) {
1184 if(polygonNTetras[i] > 2) {
1185 validTetraIds.emplace_back(i);
1190 size_t nnewpoints{};
1191 std::vector<SimplexId> pointsPerCell(validTetraIds.size() + 1);
1192 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1193 nnewpoints += polygonNTetras[validTetraIds[i]];
1194 pointsPerCell[i + 1] = nnewpoints;
1198 outSeps2.cl.connectivity_.resize(firstCellId + nnewpoints);
1199 auto cellsConn = &outSeps2.cl.connectivity_[firstCellId];
1201 std::vector<SimplexId> cellVertsIds(nnewpoints);
1203#ifdef TTK_ENABLE_OPENMP
1204#pragma omp parallel for num_threads(threadNumber_)
1206 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1207 const auto k = validTetraIds[i];
1210 getDualPolygon(polygonEdgeIds[k], &cellVertsIds[pointsPerCell[i]],
1211 polygonNTetras[k], triangulation);
1213 sortDualPolygonVertices(
1214 &cellVertsIds[pointsPerCell[i]], polygonNTetras[k], triangulation);
1216 for(
SimplexId j = 0; j < polygonNTetras[k]; ++j) {
1217 cellsConn[pointsPerCell[i] + j] = cellVertsIds[pointsPerCell[i] + j];
1221 TTK_PSORT(this->threadNumber_, cellVertsIds.begin(), cellVertsIds.end());
1222 const auto last = std::unique(cellVertsIds.begin(), cellVertsIds.end());
1223 cellVertsIds.erase(last, cellVertsIds.end());
1226 std::vector<SimplexId> vertId2PointsId(triangulation.getNumberOfCells());
1228 const auto noldpoints{npoints};
1229 npoints += cellVertsIds.size();
1230 ncells = noldcells + validTetraIds.size();
1233 outSeps2.pt.points_.resize(3 * npoints);
1234 auto points = &outSeps2.pt.points_[3 * noldpoints];
1235 outSeps2.cl.offsets_.resize(ncells + 1);
1236 outSeps2.cl.offsets_[0] = 0;
1237 auto cellsOff = &outSeps2.cl.offsets_[noldcells];
1238 outSeps2.cl.sourceIds_.resize(ncells);
1239 outSeps2.cl.separatrixIds_.resize(ncells);
1240 outSeps2.cl.separatrixTypes_.resize(ncells);
1241 outSeps2.cl.isOnBoundary_.resize(ncells);
1243#ifdef TTK_ENABLE_OPENMP
1244#pragma omp parallel for num_threads(threadNumber_)
1246 for(
size_t i = 0; i < cellVertsIds.size(); ++i) {
1248 triangulation.getTetraIncenter(cellVertsIds[i], &points[3 * i]);
1250 vertId2PointsId[cellVertsIds[i]] = i + noldpoints;
1253#ifdef TTK_ENABLE_OPENMP
1254#pragma omp parallel for num_threads(threadNumber_)
1256 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1257 const auto m = validTetraIds[i];
1258 const auto k = pointsPerCell[i];
1259 for(
SimplexId j = 0; j < polygonNTetras[m]; ++j) {
1260 cellsConn[k + j] = vertId2PointsId[cellsConn[k + j]];
1262 const auto l = i + noldcells;
1263 const auto n = polygonSepInfosIds[m];
1264 outSeps2.cl.sourceIds_[l] = sepSourceIds[n];
1265 outSeps2.cl.separatrixIds_[l] = sepIds[n];
1266 outSeps2.cl.separatrixTypes_[l] = 1;
1267 outSeps2.cl.isOnBoundary_[l] = sepOnBoundary[n];
1270 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1272 cellsOff[i + 1] = cellsOff[i] + polygonNTetras[validTetraIds[i]];
1275 outSeps2.pt.numberOfPoints_ = npoints;
1276 outSeps2.cl.numberOfCells_ = ncells;
1281template <
typename triangulationType>
1284 const std::vector<Separatrix> &separatrices,
1285 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
1287 const triangulationType &triangulation)
const {
1305 const auto noldcells{ncells};
1310 std::vector<size_t> geomCellsBegId{ncells};
1313 for(
size_t i = 0; i < separatrices.size(); ++i) {
1314 const auto &sep = separatrices[i];
1315 ncells += sep.geometry_.size();
1316 geomCellsBegId.emplace_back(ncells);
1323 firstCellId + 3 * (ncells - noldcells));
1324 auto cellsOff = &outSeps2.
cl.
offsets_[noldcells];
1329 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
1330 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
1334 std::vector<SimplexId> cellVertsIds(3 * (ncells - noldcells));
1336#ifdef TTK_ENABLE_OPENMP
1337#pragma omp parallel for num_threads(threadNumber_)
1339 for(
size_t i = 0; i < separatrices.size(); ++i) {
1340 const auto &sep = separatrices[i];
1341 const auto &sepGeom = sep.geometry_;
1342 const auto &sepSaddles = separatricesSaddles[i];
1343 const auto sepId = separatrixId + i;
1345 const char sepType = 2;
1348 const auto sepFuncMax
1349 = discreteGradient_.getCellGreaterVertex(src, triangulation);
1351 if(!sepSaddles.empty()) {
1353 const auto minId = *std::min_element(
1354 sepSaddles.begin(), sepSaddles.end(),
1356 return offsets[discreteGradient_.getCellLowerVertex(
1357 Cell{1, a}, triangulation)]
1358 < offsets[discreteGradient_.getCellLowerVertex(
1359 Cell{1, b}, triangulation)];
1362 = discreteGradient_.getCellLowerVertex(
Cell{1, minId}, triangulation);
1366 const auto minId = *std::min_element(
1367 sepGeom.begin(), sepGeom.end(),
1368 [&triangulation, offsets,
this](
const Cell &a,
const Cell &b) {
1369 return offsets[discreteGradient_.getCellLowerVertex(a, triangulation)]
1370 < offsets[discreteGradient_.getCellLowerVertex(
1373 sepFuncMin = discreteGradient_.getCellLowerVertex(minId, triangulation);
1375 separatrixFunctionMaxima[sepId] = sepFuncMax;
1376 separatrixFunctionMinima[sepId] = sepFuncMin;
1379 const char onBoundary
1380 = (sepSaddles.empty()
1382 : std::count_if(sepSaddles.begin(), sepSaddles.end(),
1384 return triangulation.isEdgeOnBoundary(a);
1386 + triangulation.isTriangleOnBoundary(src.id_);
1388 for(
size_t j = 0; j < sepGeom.size(); ++j) {
1389 const auto &cell = sepGeom[j];
1393 triangulation.getTriangleVertex(cell.id_, 0, v0);
1394 triangulation.getTriangleVertex(cell.id_, 1, v1);
1395 triangulation.getTriangleVertex(cell.id_, 2, v2);
1398 const auto l = geomCellsBegId[i] + j;
1400 const auto m = l - noldcells;
1402 cellsConn[3 * m + 0] = v0;
1403 cellsConn[3 * m + 1] = v1;
1404 cellsConn[3 * m + 2] = v2;
1405 cellVertsIds[3 * m + 0] = v0;
1406 cellVertsIds[3 * m + 1] = v1;
1407 cellVertsIds[3 * m + 2] = v2;
1418 TTK_PSORT(this->threadNumber_, cellVertsIds.begin(), cellVertsIds.end());
1419 const auto last = std::unique(cellVertsIds.begin(), cellVertsIds.end());
1420 cellVertsIds.erase(last, cellVertsIds.end());
1423 std::vector<size_t> vertId2PointsId(triangulation.getNumberOfVertices());
1425 const auto noldpoints{npoints};
1426 npoints += cellVertsIds.size();
1427 outSeps2.pt.points_.resize(3 * npoints);
1428 auto points = &outSeps2.pt.points_[3 * noldpoints];
1430#ifdef TTK_ENABLE_OPENMP
1431#pragma omp parallel for num_threads(threadNumber_)
1433 for(
size_t i = 0; i < cellVertsIds.size(); ++i) {
1435 triangulation.getVertexPoint(
1436 cellVertsIds[i], points[3 * i + 0], points[3 * i + 1], points[3 * i + 2]);
1438 vertId2PointsId[cellVertsIds[i]] = i + noldpoints;
1441 const auto lastOffset = noldcells == 0 ? 0 : cellsOff[-1];
1443#ifdef TTK_ENABLE_OPENMP
1444#pragma omp parallel for num_threads(threadNumber_)
1446 for(
size_t i = 0; i < ncells - noldcells; ++i) {
1447 cellsOff[i] = 3 * i + lastOffset;
1448 cellsConn[3 * i + 0] = vertId2PointsId[cellsConn[3 * i + 0]];
1449 cellsConn[3 * i + 1] = vertId2PointsId[cellsConn[3 * i + 1]];
1450 cellsConn[3 * i + 2] = vertId2PointsId[cellsConn[3 * i + 2]];
1453 cellsOff[ncells - noldcells] = cellsOff[ncells - noldcells - 1] + 3;
1455 outSeps2.pt.numberOfPoints_ = npoints;
1456 outSeps2.cl.numberOfCells_ = ncells;
1465template <
typename triangulationType>
1467 const std::vector<SimplexId> &maxima,
1469 const triangulationType &triangulation)
const {
1471 if(morseSmaleManifold ==
nullptr) {
1472 this->printErr(
"Could not compute ascending segmentation");
1478 const auto nVerts{triangulation.getNumberOfVertices()};
1479 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, -1);
1480 if(maxima.empty()) {
1485 const auto nCells{triangulation.getNumberOfCells()};
1486 std::vector<SimplexId> morseSmaleManifoldOnCells(nCells, -1);
1489 for(
const auto &
id : maxima) {
1491 morseSmaleManifoldOnCells[id] = nMax++;
1494 const auto dim{triangulation.getDimensionality()};
1497 auto getFaceStarNumber = &triangulationType::getTriangleStarNumber;
1498 auto getFaceStar = &triangulationType::getTriangleStar;
1501 getFaceStarNumber = &triangulationType::getEdgeStarNumber;
1502 getFaceStar = &triangulationType::getEdgeStar;
1503 }
else if(dim == 1) {
1505 getFaceStarNumber = &triangulationType::getVertexStarNumber;
1506 getFaceStar = &triangulationType::getVertexStar;
1510 std::vector<SimplexId> visited{};
1512 std::vector<uint8_t> isMarked(nCells, 0);
1514#ifdef TTK_ENABLE_OPENMP
1515#pragma omp parallel for num_threads(threadNumber_) firstprivate(visited)
1518 if(isMarked[i] == 1) {
1523 while(morseSmaleManifoldOnCells[curr] == -1) {
1524 if(isMarked[curr] == 1) {
1528 const auto paired{this->discreteGradient_.getPairedCell(
1529 Cell{dim, curr}, triangulation,
true)};
1531 const auto nStars{(triangulation.*getFaceStarNumber)(paired)};
1533 (triangulation.*getFaceStar)(paired, j, next);
1539 visited.emplace_back(curr);
1546 for(
const auto el : visited) {
1547 morseSmaleManifoldOnCells[el] = morseSmaleManifoldOnCells[curr];
1552#ifdef TTK_ENABLE_OPENMP
1553#pragma omp parallel for num_threads(threadNumber_)
1556 if(triangulation.getVertexStarNumber(i) < 1) {
1561 triangulation.getVertexStar(i, 0, starId);
1563 morseSmaleManifold[i] = morseSmaleManifoldOnCells[starId];
1566 this->
printMsg(
" Ascending segmentation computed", 1.0, tm.getElapsedTime(),
1573template <
typename triangulationType>
1575 const std::vector<SimplexId> &minima,
1577 const triangulationType &triangulation)
const {
1579 if(morseSmaleManifold ==
nullptr) {
1580 this->printErr(
"Could not compute descending segmentation");
1586 const auto nVerts{triangulation.getNumberOfVertices()};
1588 if(minima.size() == 1) {
1590 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, 0);
1594 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, -1);
1597 for(
const auto &cp : minima) {
1599 morseSmaleManifold[cp] = nMin++;
1602 std::vector<SimplexId> visited{};
1604#ifdef TTK_ENABLE_OPENMP
1605#pragma omp parallel for num_threads(threadNumber_) firstprivate(visited)
1608 if(morseSmaleManifold[i] != -1) {
1613 while(morseSmaleManifold[curr] == -1) {
1615 const auto pairedEdge{
1616 this->discreteGradient_.getPairedCell(
Cell{0, curr}, triangulation)};
1618 triangulation.getEdgeVertex(pairedEdge, 0, next);
1620 triangulation.getEdgeVertex(pairedEdge, 1, next);
1622 visited.emplace_back(curr);
1625 for(
const auto el : visited) {
1626 morseSmaleManifold[el] = morseSmaleManifold[curr];
1630 this->
printMsg(
" Descending segmentation computed", 1.0, tm.getElapsedTime(),
1637template <
typename triangulationType>
1640 const SimplexId *
const ascendingManifold,
1641 const SimplexId *
const descendingManifold,
1643 const triangulationType &triangulation)
const {
1645 if(ascendingManifold ==
nullptr || descendingManifold ==
nullptr
1646 || morseSmaleManifold ==
nullptr) {
1647 this->printErr(
"Could not compute final segmentation");
1653 const size_t nVerts = triangulation.getNumberOfVertices();
1657#ifdef TTK_ENABLE_OPENMP
1658#pragma omp parallel for num_threads(threadNumber_)
1660 for(
size_t i = 0; i < nVerts; ++i) {
1661 const auto d = ascendingManifold[i];
1662 const auto a = descendingManifold[i];
1663 if(a == -1 || d == -1) {
1664 morseSmaleManifold[i] = -1;
1666 morseSmaleManifold[i] = a * numberOfMaxima + d;
1671 std::vector<SimplexId> sparseRegionIds(
1672 morseSmaleManifold, morseSmaleManifold + nVerts);
1676 this->threadNumber_, sparseRegionIds.begin(), sparseRegionIds.end());
1677 const auto last = std::unique(sparseRegionIds.begin(), sparseRegionIds.end());
1678 sparseRegionIds.erase(last, sparseRegionIds.end());
1681 std::map<SimplexId, size_t> sparseToDenseRegionId{};
1683 for(
size_t i = 0; i < sparseRegionIds.size(); ++i) {
1684 sparseToDenseRegionId[sparseRegionIds[i]] = i;
1689#ifdef TTK_ENABLE_OPENMP
1690#pragma omp parallel for num_threads(threadNumber_)
1692 for(
size_t i = 0; i < nVerts; ++i) {
1693 morseSmaleManifold[i] = sparseToDenseRegionId[morseSmaleManifold[i]];
1696 this->
printMsg(
" Final segmentation computed", 1.0, tm.getElapsedTime(),
1703template <
typename dataType,
typename triangulationType>
1705 const double persistenceThreshold,
1706 const dataType *
const scalars,
1708 const triangulationType &triangulation) {
1712 const auto dim{triangulation.getDimensionality()};
1714 this->printWrn(
"Can't return saddle connectors without a 3D dataset");
1721 dms.setDebugLevel(this->debugLevel_);
1722 dms.setGradient(std::move(this->discreteGradient_));
1726 std::vector<PersPairType> dms_pairs{};
1727 dms.computePersistencePairs(dms_pairs, offsets, triangulation,
false,
true);
1728 this->discreteGradient_ = dms.getGradient();
1730 this->discreteGradient_.setLocalGradient();
1732 const auto getPersistence
1733 = [
this, &triangulation, scalars](
const PersPairType &p) {
1734 return this->discreteGradient_.getPersistence(
1735 Cell{2, p.death},
Cell{1, p.birth}, scalars, triangulation);
1739 const auto firstSadSadPair{std::distance(
1741 std::find_if(dms_pairs.begin(), dms_pairs.end(),
1742 [](
const auto &pair) { return pair.type == 1; }))};
1746 std::vector<bool> isVisited(triangulation.getNumberOfTriangles(),
false);
1747 std::vector<SimplexId> visitedTriangles{};
1752 const auto &s2Children{dms.get2SaddlesChildren()};
1753 std::vector<bool> simplifyS2(s2Children.size(),
true);
1755 for(
size_t i = firstSadSadPair; i < dms_pairs.size(); ++i) {
1756 const auto &pair{dms_pairs[i]};
1758 if(pair.type != 1 || getPersistence(pair) > persistenceThreshold) {
1762 const auto o{i - firstSadSadPair};
1763 const Cell birth{1, pair.birth};
1764 const Cell death{2, pair.death};
1766 for(
const auto child : s2Children[o]) {
1767 if(!simplifyS2[child]) {
1773 this->
printMsg(
"Skipping saddle connector " + birth.to_string() +
" -> "
1774 + death.to_string(),
1776 simplifyS2[o] =
false;
1781 this->discreteGradient_.getDescendingWall(death, mask, triangulation);
1783 std::vector<Cell> vpath{};
1784 this->discreteGradient_.getAscendingPathThroughWall(
1785 birth, death, isVisited, &vpath, triangulation,
true);
1787 if(vpath.back() == death) {
1788 this->discreteGradient_.reverseAscendingPathOnWall(vpath, triangulation);
1791 this->
printMsg(
"Could not return saddle connector " + birth.to_string()
1792 +
" -> " + death.to_string(),
1794 simplifyS2[o] =
false;
1798 this->
printMsg(
"Returned " + std::to_string(nReturned) +
" saddle connectors",
1799 1.0, tm.getElapsedTime(), this->threadNumber_);
#define TTK_PSORT(NTHREADS,...)
Parallel sort macro.
AbstractTriangulation is an interface class that defines an interface for efficient traversal methods...
virtual int preconditionCellNeighbors()
virtual int preconditionCellEdges()
virtual int setThreadNumber(const int threadNumber)
Minimalist debugging class.
virtual int setDebugLevel(const int &debugLevel)
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
int printErr(const std::string &msg, const debug::LineMode &lineMode=debug::LineMode::NEW, std::ostream &stream=std::cerr) const
TTK DiscreteMorseSandwich processing package.
TTK processing package for the computation of Morse-Smale complexes.
int setAscendingSegmentation(const std::vector< SimplexId > &maxima, SimplexId *const morseSmaleManifold, const triangulationType &triangulation) const
int setFinalSegmentation(const SimplexId numberOfMaxima, const SimplexId *const ascendingManifold, const SimplexId *const descendingManifold, SimplexId *const morseSmaleManifold, const triangulationType &triangulation) const
bool ComputeSaddleConnectors
void flattenSeparatricesVectors(std::vector< std::vector< Separatrix > > &separatrices) const
Flatten the vectors of vectors into their first component.
bool ComputeDescendingSeparatrices2
bool ComputeDescendingSegmentation
int getDualPolygon(const SimplexId edgeId, SimplexId *const polygon, const size_t polSize, const triangulationType &triangulation) const
void setComputeCriticalPoints(const bool state)
int getAscendingSeparatrices2(const std::vector< SimplexId > &saddles1, std::vector< Separatrix > &separatrices, std::vector< std::vector< SimplexId > > &separatricesSaddles, const triangulationType &triangulation) const
int getSaddleConnectors(const std::vector< SimplexId > &saddles2, std::vector< Separatrix > &separatrices, const triangulationType &triangulation) const
int getDescendingSeparatrices1(const std::vector< SimplexId > &saddles, std::vector< Separatrix > &separatrices, const triangulationType &triangulation) const
dcg::DiscreteGradient discreteGradient_
int returnSaddleConnectors(const double persistenceThreshold, const dataType *const scalars, const SimplexId *const offsets, const triangulationType &triangulation)
int setAscendingSeparatrices2(Output2Separatrices &outSeps2, const std::vector< Separatrix > &separatrices, const std::vector< std::vector< SimplexId > > &separatricesSaddles, const SimplexId *const offsets, const triangulationType &triangulation) const
bool ComputeAscendingSeparatrices1
bool ComputeDescendingSeparatrices1
int setDescendingSegmentation(const std::vector< SimplexId > &minima, SimplexId *const morseSmaleManifold, const triangulationType &triangulation) const
int setSeparatrices1(Output1Separatrices &outSeps1, const std::vector< Separatrix > &separatrices, const SimplexId *const offsets, const triangulationType &triangulation) const
bool ComputeAscendingSeparatrices2
bool ComputeAscendingSegmentation
void setComputeSegmentation(const bool doAscending, const bool doDescending, const bool doMorseSmale)
void setSaddleConnectorsPersistenceThreshold(const double threshold)
int sortDualPolygonVertices(SimplexId *const polygon, const size_t polSize, const triangulationType &triangulation) const
void setComputeSeparatrices2(const bool doAscending, const bool doDescending)
void preconditionTriangulation(AbstractTriangulation *const data)
int setDescendingSeparatrices2(Output2Separatrices &outSeps2, const std::vector< Separatrix > &separatrices, const std::vector< std::vector< SimplexId > > &separatricesSaddles, const SimplexId *const offsets, const triangulationType &triangulation) const
int getDescendingSeparatrices2(const std::vector< SimplexId > &saddles2, std::vector< Separatrix > &separatrices, std::vector< std::vector< SimplexId > > &separatricesSaddles, const triangulationType &triangulation) const
bool ComputeFinalSegmentation
void setComputeSeparatrices1(const bool doAscending, const bool doDescending, const bool doSaddleConnectors)
void setReturnSaddleConnectors(const bool state)
bool ComputeCriticalPoints
int execute(OutputCriticalPoints &outCP, Output1Separatrices &outSeps1, Output2Separatrices &outSeps2, OutputManifold &outManifold, const dataType *const scalars, const size_t scalarsMTime, const SimplexId *const offsets, const triangulationType &triangulation)
bool ReturnSaddleConnectors
double SaddleConnectorsPersistenceThreshold
int getAscendingSeparatrices1(const std::vector< SimplexId > &saddles, std::vector< Separatrix > &separatrices, const triangulationType &triangulation) const
TTK discreteGradient processing package.
void setInputOffsets(const SimplexId *const data)
int setManifoldSize(const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const SimplexId *const ascendingManifold, const SimplexId *const descendingManifold, std::vector< SimplexId > &manifoldSize) const
void setInputScalarField(const void *const data, const size_t mTime)
int getCriticalPoints(std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, const triangulationType &triangulation) const
void preconditionTriangulation(AbstractTriangulation *const data)
int buildGradient(const triangulationType &triangulation, bool bypassCache=false)
int setCriticalPoints(const std::array< std::vector< SimplexId >, 4 > &criticalCellsByDim, std::vector< std::array< float, 3 > > &points, std::vector< char > &cellDimensions, std::vector< SimplexId > &cellIds, std::vector< char > &isOnBoundary, std::vector< SimplexId > &PLVertexIdentifiers, const triangulationType &triangulation) const
int SimplexId
Identifier type for simplices of any dimension.
Persistence pair struct as exported by DiscreteGradient.
1-Separatrices point and cell data arrays
std::vector< ttk::SimplexId > destinationIds_
std::vector< SimplexId > sepFuncMinId_
std::vector< SimplexId > sepFuncMaxId_
SimplexId numberOfPoints_
std::vector< ttk::SimplexId > separatrixIds_
struct ttk::MorseSmaleComplex::Output1Separatrices::@1 cl
std::vector< char > isOnBoundary_
std::vector< ttk::SimplexId > connectivity_
std::vector< ttk::SimplexId > sourceIds_
std::vector< char > smoothingMask_
std::vector< char > separatrixTypes_
struct ttk::MorseSmaleComplex::Output1Separatrices::@0 pt
std::vector< char > cellDimensions_
std::vector< ttk::SimplexId > cellIds_
std::vector< float > points_
2-Separatrices point and cell data arrays
std::vector< ttk::SimplexId > offsets_
SimplexId numberOfPoints_
std::vector< float > points_
struct ttk::MorseSmaleComplex::Output2Separatrices::@2 pt
std::vector< SimplexId > sepFuncMinId_
std::vector< ttk::SimplexId > sourceIds_
std::vector< char > isOnBoundary_
struct ttk::MorseSmaleComplex::Output2Separatrices::@3 cl
std::vector< ttk::SimplexId > connectivity_
std::vector< char > separatrixTypes_
std::vector< ttk::SimplexId > separatrixIds_
std::vector< SimplexId > sepFuncMaxId_
Critical points data arrays.
std::vector< char > cellDimensions_
std::vector< char > isOnBoundary_
std::vector< SimplexId > cellIds_
std::vector< SimplexId > manifoldSize_
std::vector< SimplexId > PLVertexIdentifiers_
std::vector< std::array< float, 3 > > points_
Pointers to pre-allocated segmentation point data arrays.
std::vector< dcg::Cell > geometry_
Auto-cleaning re-usable graph propagations data structure.
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)