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);
424template <
typename dataType,
typename triangulationType>
429 const dataType *
const scalars,
430 const size_t scalarsMTime,
432 const triangulationType &triangulation) {
433#ifndef TTK_ENABLE_KAMIKAZE
434 if(scalars ==
nullptr) {
435 this->
printErr(
"Input scalar field pointer is null.");
439 if(offsets ==
nullptr) {
440 this->
printErr(
"Input offset field pointer is null.");
449 const auto dim = triangulation.getDimensionality();
461 const auto nVerts{triangulation.getNumberOfVertices()};
464 const auto pair{std::minmax_element(offsets, offsets + nVerts)};
466 const auto globmin = std::distance(offsets, pair.first);
467 const auto globmax = std::distance(offsets, pair.second);
468 persistenceThreshold *= (scalars[globmax] - scalars[globmin]);
469 this->
printMsg(
"Absolute saddle connectors persistence threshold is "
475 persistenceThreshold, scalars, offsets, triangulation);
478 std::array<std::vector<SimplexId>, 4> criticalPoints{};
482 this->
printMsg(
" Critical points extracted", 1.0, tm.getElapsedTime(),
487 std::vector<std::vector<Separatrix>> separatrices1{};
494 separatrices1.emplace_back();
497 criticalPoints[1], separatrices1.back(), triangulation);
499 this->
printMsg(
" Descending 1-separatrices computed", 1.0,
506 separatrices1.emplace_back();
509 criticalPoints[dim - 1], separatrices1.back(), triangulation);
511 this->
printMsg(
" Ascending 1-separatrices computed", 1.0,
519 separatrices1.emplace_back();
536 this->
printMsg(
" 1-separatrices set", 1.0, tmp.getElapsedTime(),
540 this->
printMsg(
"1-separatrices computed", 1.0, tm1sep.getElapsedTime(),
541 this->threadNumber_);
549 std::vector<Separatrix> separatrices;
550 std::vector<std::vector<SimplexId>> separatricesSaddles;
552 criticalPoints[2], separatrices, separatricesSaddles, triangulation);
554 outSeps2, separatrices, separatricesSaddles, offsets, triangulation);
556 this->
printMsg(
" Descending 2-separatrices computed", 1.0,
563 std::vector<Separatrix> separatrices;
564 std::vector<std::vector<SimplexId>> separatricesSaddles;
566 criticalPoints[1], separatrices, separatricesSaddles, triangulation);
568 outSeps2, separatrices, separatricesSaddles, offsets, triangulation);
570 this->
printMsg(
" Ascending 2-separatrices computed", 1.0,
577 this->
printMsg(
"2-separatrices computed", 1.0, tm2sep.getElapsedTime(),
578 this->threadNumber_);
586 criticalPoints[dim], outManifold.
ascending_, triangulation);
590 criticalPoints[0], outManifold.
descending_, triangulation);
600 "Segmentation computed", 1.0, tmp.
getElapsedTime(), this->threadNumber_);
617 +
" points) processed",
627template <
typename triangulationType>
629 const std::vector<SimplexId> &saddles,
630 std::vector<Separatrix> &separatrices,
631 const triangulationType &triangulation)
const {
633 const SimplexId numberOfSaddles = saddles.size();
636 const SimplexId numberOfSeparatrices = 2 * numberOfSaddles;
637 separatrices.resize(numberOfSeparatrices);
640#ifdef TTK_ENABLE_OPENMP
641#pragma omp parallel for num_threads(threadNumber_)
643 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
644 const Cell saddle{1, saddles[i]};
648 const Cell &saddle1 = saddle;
650 for(
int j = 0; j < 2; ++j) {
652 triangulation.getEdgeVertex(saddle1.
id_, j, vertexId);
654 std::vector<Cell> vpath;
655 vpath.push_back(saddle1);
656 discreteGradient_.getDescendingPath(
657 Cell(0, vertexId), vpath, triangulation);
659 const Cell &lastCell = vpath.back();
660 if(lastCell.
dim_ == 0 and discreteGradient_.isCellCritical(lastCell)) {
661 separatrices[2 * i + j].source_ = saddle;
662 separatrices[2 * i + j].destination_ = lastCell;
663 separatrices[2 * i + j].geometry_ = std::move(vpath);
672template <
typename triangulationType>
674 const std::vector<SimplexId> &saddles,
675 std::vector<Separatrix> &separatrices,
676 const triangulationType &triangulation)
const {
681 const auto dim{triangulation.getDimensionality()};
684 auto getFaceStarNumber = &triangulationType::getTriangleStarNumber;
685 auto getFaceStar = &triangulationType::getTriangleStar;
688 getFaceStarNumber = &triangulationType::getEdgeStarNumber;
689 getFaceStar = &triangulationType::getEdgeStar;
692 const SimplexId numberOfSaddles = saddles.size();
694 std::vector<std::vector<Separatrix>> sepsPerSaddle(numberOfSaddles);
697#ifdef TTK_ENABLE_OPENMP
698#pragma omp parallel for num_threads(threadNumber_)
700 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
701 const Cell saddle{dim - 1, saddles[i]};
704 const auto starNumber{(triangulation.*getFaceStarNumber)(saddle.id_)};
705 for(
SimplexId j = 0; j < starNumber; ++j) {
708 (triangulation.*getFaceStar)(saddle.id_, j, sId);
710 std::vector<Cell> vpath{saddle};
711 discreteGradient_.getAscendingPath(
Cell(dim, sId), vpath, triangulation);
713 const Cell &lastCell = vpath.back();
714 if(lastCell.
dim_ == dim and discreteGradient_.isCellCritical(lastCell)) {
715 sepsPerSaddle[i].emplace_back();
716 sepsPerSaddle[i].back().source_ = saddle;
717 sepsPerSaddle[i].back().destination_ = lastCell;
718 sepsPerSaddle[i].back().geometry_ = std::move(vpath);
723 this->flattenSeparatricesVectors(sepsPerSaddle);
725 separatrices = std::move(sepsPerSaddle[0]);
730template <
typename triangulationType>
732 const std::vector<SimplexId> &saddles2,
733 std::vector<Separatrix> &separatrices,
734 const triangulationType &triangulation)
const {
739 const auto nTriangles = triangulation.getNumberOfTriangles();
741 std::vector<bool> isVisited(nTriangles,
false);
742 std::vector<SimplexId> visitedTriangles{};
744 using Vpath = std::vector<Cell>;
746 const auto dim{triangulation.getDimensionality()};
748 std::vector<std::vector<Separatrix>> sepsByThread(saddles2.size());
749 std::vector<SimplexId> saddles1{};
751#ifdef TTK_ENABLE_OPENMP
752#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
753 firstprivate(isVisited, visitedTriangles, saddles1)
755 for(
size_t i = 0; i < saddles2.size(); ++i) {
756 const Cell s2{dim - 1, saddles2[i]};
759 discreteGradient_.getDescendingWall(
760 s2, mask, triangulation,
nullptr, &saddles1);
762 for(
const auto saddle1Id : saddles1) {
763 const Cell s1{1, saddle1Id};
766 const bool isMultiConnected
767 = discreteGradient_.getAscendingPathThroughWall(
768 s1, s2, isVisited, &vpath, triangulation);
774 const auto &last = vpath.back();
776 if(!isMultiConnected && last.dim_ == s2.dim_ && last.id_ == s2.id_) {
777 sepsByThread[i].emplace_back();
778 sepsByThread[i].back().source_ = s1;
779 sepsByThread[i].back().destination_ = s2;
780 sepsByThread[i].back().geometry_ = std::move(vpath);
785 this->flattenSeparatricesVectors(sepsByThread);
787 separatrices = std::move(sepsByThread[0]);
792template <
typename triangulationType>
795 const std::vector<Separatrix> &separatrices,
797 const triangulationType &triangulation)
const {
815 std::vector<size_t> geomPointsBegId{npoints};
817 std::vector<size_t> geomCellsBegId{ncells};
820 for(
size_t i = 0; i < separatrices.size(); ++i) {
821 const auto &sep = separatrices[i];
822 const auto sepSize = sep.geometry_.size();
824 ncells += sepSize - 1;
825 geomPointsBegId.emplace_back(npoints);
826 geomCellsBegId.emplace_back(ncells);
829 const int dimensionality = triangulation.getCellVertexNumber(0) - 1;
843 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
844 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
847#ifdef TTK_ENABLE_OPENMP
848#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic)
850 for(
size_t i = 0; i < separatrices.size(); ++i) {
851 const auto &sep = separatrices[i];
852 const auto &sepGeom = sep.geometry_;
853 const auto sepId = separatrixId + i;
860 const auto saddleConnector
861 = dimensionality == 3 && src.
dim_ == 1 && dst.
dim_ == 2;
863 = saddleConnector ? 1 : std::min(dst.
dim_, dimensionality - 1);
867 return offsets[a] < offsets[b];
869 const std::array<SimplexId, 2> gVerts{
870 discreteGradient_.getCellGreaterVertex(src, triangulation),
871 discreteGradient_.getCellGreaterVertex(dst, triangulation)};
872 const auto sepFuncMax
873 = *std::max_element(gVerts.begin(), gVerts.end(), vertsOrder);
874 const std::array<SimplexId, 2> lVerts{
875 discreteGradient_.getCellLowerVertex(src, triangulation),
876 discreteGradient_.getCellLowerVertex(dst, triangulation)};
877 const auto sepFuncMin
878 = *std::min_element(lVerts.begin(), lVerts.end(), vertsOrder);
879 separatrixFunctionMaxima[sepId] = sepFuncMax;
880 separatrixFunctionMinima[sepId] = sepFuncMin;
883 const auto onBoundary
884 =
static_cast<char>(discreteGradient_.isBoundary(src, triangulation))
885 +
static_cast<char>(discreteGradient_.isBoundary(dst, triangulation));
887 for(
size_t j = 0; j < sepGeom.size(); ++j) {
888 const auto &cell = sepGeom[j];
889 std::array<float, 3> pt{};
890 triangulation.getCellIncenter(cell.id_, cell.dim_, pt.data());
893 const auto k = geomPointsBegId[i] + j;
895 points[3 * k + 0] = pt[0];
896 points[3 * k + 1] = pt[1];
897 points[3 * k + 2] = pt[2];
900 = (j == 0 || j == sepGeom.size() - 1) ? 0 : 1;
909 const auto l = geomCellsBegId[i] + j - 1;
911 cellsConn[2 * l + 0] = k - 1;
912 cellsConn[2 * l + 1] = k;
933template <
typename triangulationType>
935 const std::vector<SimplexId> &saddles1,
936 std::vector<Separatrix> &separatrices,
937 std::vector<std::vector<SimplexId>> &separatricesSaddles,
938 const triangulationType &triangulation)
const {
939 const Cell emptyCell;
941 const SimplexId numberOfSaddles = saddles1.size();
945 const SimplexId numberOfSeparatrices = numberOfSaddles;
946 separatrices.resize(numberOfSeparatrices);
947 separatricesSaddles.resize(numberOfSeparatrices);
949 const auto nEdges = triangulation.getNumberOfEdges();
950 std::vector<bool> isVisited(nEdges,
false);
951 std::vector<SimplexId> visitedEdges{};
954#ifdef TTK_ENABLE_OPENMP
955#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
956 firstprivate(isVisited, visitedEdges)
958 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
959 const Cell saddle1{1, saddles1[i]};
961 std::vector<Cell> wall;
963 discreteGradient_.getAscendingWall(
964 saddle1, mask, triangulation, &wall, &separatricesSaddles[i]);
966 separatrices[i].source_ = saddle1;
967 separatrices[i].destination_ = emptyCell;
968 separatrices[i].geometry_ = std::move(wall);
974template <
typename triangulationType>
976 const std::vector<SimplexId> &saddles2,
977 std::vector<Separatrix> &separatrices,
978 std::vector<std::vector<SimplexId>> &separatricesSaddles,
979 const triangulationType &triangulation)
const {
980 const Cell emptyCell;
982 const SimplexId numberOfSaddles = saddles2.size();
986 const SimplexId numberOfSeparatrices = numberOfSaddles;
987 separatrices.resize(numberOfSeparatrices);
988 separatricesSaddles.resize(numberOfSeparatrices);
990 const auto nTriangles = triangulation.getNumberOfTriangles();
991 std::vector<bool> isVisited(nTriangles,
false);
992 std::vector<SimplexId> visitedTriangles{};
994 const auto dim{triangulation.getDimensionality()};
997#ifdef TTK_ENABLE_OPENMP
998#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) \
999 firstprivate(isVisited, visitedTriangles)
1001 for(
SimplexId i = 0; i < numberOfSaddles; ++i) {
1002 const Cell saddle2{dim - 1, saddles2[i]};
1004 std::vector<Cell> wall;
1006 discreteGradient_.getDescendingWall(
1007 saddle2, mask, triangulation, &wall, &separatricesSaddles[i]);
1009 separatrices[i].source_ = saddle2;
1010 separatrices[i].destination_ = emptyCell;
1011 separatrices[i].geometry_ = std::move(wall);
1017template <
typename triangulationType>
1021 const size_t polSize,
1022 const triangulationType &triangulation)
const {
1024 for(
size_t i = 0; i < polSize; ++i) {
1026 triangulation.getEdgeStar(edgeId, i, starId);
1027 polygon[i] = starId;
1033template <
typename triangulationType>
1036 const size_t polSize,
1037 const triangulationType &triangulation)
const {
1039 for(
size_t i = 1; i < polSize; ++i) {
1042 bool isFound =
false;
1044 for(; j < polSize; ++j) {
1047 k < triangulation.getCellNeighborNumber(polygon[i - 1]); ++k) {
1049 triangulation.getCellNeighbor(polygon[i - 1], k, neighborId);
1050 if(neighborId == polygon[j]) {
1061 std::swap(polygon[j], polygon[i]);
1068template <
typename triangulationType>
1071 const std::vector<Separatrix> &separatrices,
1072 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
1074 const triangulationType &triangulation)
const {
1092 const auto noldcells{ncells};
1096 std::vector<size_t> geomCellsBegId{ncells};
1099 for(
size_t i = 0; i < separatrices.size(); ++i) {
1100 const auto &sep = separatrices[i];
1101 ncells += sep.geometry_.size();
1102 geomCellsBegId.emplace_back(ncells);
1106 std::vector<SimplexId> sepSourceIds(separatrices.size());
1107 std::vector<SimplexId> sepIds(separatrices.size());
1108 std::vector<char> sepOnBoundary(separatrices.size());
1109 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
1110 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
1112 std::vector<SimplexId> polygonNTetras(ncells - noldcells);
1113 std::vector<SimplexId> polygonEdgeIds(ncells - noldcells);
1114 std::vector<SimplexId> polygonSepInfosIds(ncells - noldcells);
1116#ifdef TTK_ENABLE_OPENMP
1117#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic)
1119 for(
size_t i = 0; i < separatrices.size(); ++i) {
1120 const auto &sep = separatrices[i];
1121 const auto &sepGeom = sep.geometry_;
1122 const auto &sepSaddles = separatricesSaddles[i];
1123 const auto sepId = separatrixId + i;
1127 const auto sepFuncMin
1128 = discreteGradient_.getCellLowerVertex(src, triangulation);
1130 if(!sepSaddles.empty()) {
1132 const auto maxId = *std::max_element(
1133 sepSaddles.begin(), sepSaddles.end(),
1135 return offsets[discreteGradient_.getCellGreaterVertex(
1136 Cell{2, a}, triangulation)]
1137 < offsets[discreteGradient_.getCellGreaterVertex(
1138 Cell{2, b}, triangulation)];
1141 = discreteGradient_.getCellGreaterVertex(
Cell{2, maxId}, triangulation);
1145 const auto maxId = *std::max_element(
1146 sepGeom.begin(), sepGeom.end(),
1147 [&triangulation, offsets,
this](
const Cell &a,
const Cell &b) {
1148 return offsets[discreteGradient_.getCellGreaterVertex(
1150 < offsets[discreteGradient_.getCellGreaterVertex(
1153 sepFuncMax = discreteGradient_.getCellGreaterVertex(maxId, triangulation);
1157 const char onBoundary
1158 = (sepSaddles.empty()
1160 : std::count_if(sepSaddles.begin(), sepSaddles.end(),
1162 return triangulation.isTriangleOnBoundary(a);
1164 + triangulation.isEdgeOnBoundary(src.id_);
1167 sepSourceIds[i] = src.id_;
1168 separatrixFunctionMaxima[sepId] = sepFuncMax;
1169 separatrixFunctionMinima[sepId] = sepFuncMin;
1170 sepOnBoundary[i] = onBoundary;
1172 for(
size_t j = 0; j < sepGeom.size(); ++j) {
1173 const auto &cell = sepGeom[j];
1175 const auto k = geomCellsBegId[i] + j - noldcells;
1177 polygonNTetras[k] = triangulation.getEdgeStarNumber(cell.id_);
1179 if(polygonNTetras[k] > 2) {
1180 polygonEdgeIds[k] = cell.id_;
1181 polygonSepInfosIds[k] = i;
1187 std::vector<SimplexId> validTetraIds{};
1188 validTetraIds.reserve(polygonNTetras.size());
1190 for(
size_t i = 0; i < polygonNTetras.size(); ++i) {
1191 if(polygonNTetras[i] > 2) {
1192 validTetraIds.emplace_back(i);
1197 size_t nnewpoints{};
1198 std::vector<SimplexId> pointsPerCell(validTetraIds.size() + 1);
1199 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1200 nnewpoints += polygonNTetras[validTetraIds[i]];
1201 pointsPerCell[i + 1] = nnewpoints;
1205 outSeps2.cl.connectivity_.resize(firstCellId + nnewpoints);
1206 auto cellsConn = &outSeps2.cl.connectivity_[firstCellId];
1208 std::vector<SimplexId> cellVertsIds(nnewpoints);
1210#ifdef TTK_ENABLE_OPENMP
1211#pragma omp parallel for num_threads(threadNumber_)
1213 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1214 const auto k = validTetraIds[i];
1217 getDualPolygon(polygonEdgeIds[k], &cellVertsIds[pointsPerCell[i]],
1218 polygonNTetras[k], triangulation);
1220 sortDualPolygonVertices(
1221 &cellVertsIds[pointsPerCell[i]], polygonNTetras[k], triangulation);
1223 for(
SimplexId j = 0; j < polygonNTetras[k]; ++j) {
1224 cellsConn[pointsPerCell[i] + j] = cellVertsIds[pointsPerCell[i] + j];
1228 TTK_PSORT(this->threadNumber_, cellVertsIds.begin(), cellVertsIds.end());
1229 const auto last = std::unique(cellVertsIds.begin(), cellVertsIds.end());
1230 cellVertsIds.erase(last, cellVertsIds.end());
1233 std::vector<SimplexId> vertId2PointsId(triangulation.getNumberOfCells());
1235 const auto noldpoints{npoints};
1236 npoints += cellVertsIds.size();
1237 ncells = noldcells + validTetraIds.size();
1240 outSeps2.pt.points_.resize(3 * npoints);
1241 auto points = &outSeps2.pt.points_[3 * noldpoints];
1242 outSeps2.cl.offsets_.resize(ncells + 1);
1243 outSeps2.cl.offsets_[0] = 0;
1244 auto cellsOff = &outSeps2.cl.offsets_[noldcells];
1245 outSeps2.cl.sourceIds_.resize(ncells);
1246 outSeps2.cl.separatrixIds_.resize(ncells);
1247 outSeps2.cl.separatrixTypes_.resize(ncells);
1248 outSeps2.cl.isOnBoundary_.resize(ncells);
1250#ifdef TTK_ENABLE_OPENMP
1251#pragma omp parallel for num_threads(threadNumber_)
1253 for(
size_t i = 0; i < cellVertsIds.size(); ++i) {
1255 triangulation.getTetraIncenter(cellVertsIds[i], &points[3 * i]);
1257 vertId2PointsId[cellVertsIds[i]] = i + noldpoints;
1260#ifdef TTK_ENABLE_OPENMP
1261#pragma omp parallel for num_threads(threadNumber_)
1263 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1264 const auto m = validTetraIds[i];
1265 const auto k = pointsPerCell[i];
1266 for(
SimplexId j = 0; j < polygonNTetras[m]; ++j) {
1267 cellsConn[k + j] = vertId2PointsId[cellsConn[k + j]];
1269 const auto l = i + noldcells;
1270 const auto n = polygonSepInfosIds[m];
1271 outSeps2.cl.sourceIds_[l] = sepSourceIds[n];
1272 outSeps2.cl.separatrixIds_[l] = sepIds[n];
1273 outSeps2.cl.separatrixTypes_[l] = 1;
1274 outSeps2.cl.isOnBoundary_[l] = sepOnBoundary[n];
1277 for(
size_t i = 0; i < validTetraIds.size(); ++i) {
1279 cellsOff[i + 1] = cellsOff[i] + polygonNTetras[validTetraIds[i]];
1282 outSeps2.pt.numberOfPoints_ = npoints;
1283 outSeps2.cl.numberOfCells_ = ncells;
1288template <
typename triangulationType>
1291 const std::vector<Separatrix> &separatrices,
1292 const std::vector<std::vector<SimplexId>> &separatricesSaddles,
1294 const triangulationType &triangulation)
const {
1312 const auto noldcells{ncells};
1317 std::vector<size_t> geomCellsBegId{ncells};
1320 for(
size_t i = 0; i < separatrices.size(); ++i) {
1321 const auto &sep = separatrices[i];
1322 ncells += sep.geometry_.size();
1323 geomCellsBegId.emplace_back(ncells);
1330 firstCellId + 3 * (ncells - noldcells));
1331 auto cellsOff = &outSeps2.
cl.
offsets_[noldcells];
1336 separatrixFunctionMaxima.resize(separatrixId + separatrices.size());
1337 separatrixFunctionMinima.resize(separatrixId + separatrices.size());
1341 std::vector<SimplexId> cellVertsIds(3 * (ncells - noldcells));
1343#ifdef TTK_ENABLE_OPENMP
1344#pragma omp parallel for num_threads(threadNumber_)
1346 for(
size_t i = 0; i < separatrices.size(); ++i) {
1347 const auto &sep = separatrices[i];
1348 const auto &sepGeom = sep.geometry_;
1349 const auto &sepSaddles = separatricesSaddles[i];
1350 const auto sepId = separatrixId + i;
1352 const char sepType = 2;
1355 const auto sepFuncMax
1356 = discreteGradient_.getCellGreaterVertex(src, triangulation);
1358 if(!sepSaddles.empty()) {
1360 const auto minId = *std::min_element(
1361 sepSaddles.begin(), sepSaddles.end(),
1363 return offsets[discreteGradient_.getCellLowerVertex(
1364 Cell{1, a}, triangulation)]
1365 < offsets[discreteGradient_.getCellLowerVertex(
1366 Cell{1, b}, triangulation)];
1369 = discreteGradient_.getCellLowerVertex(
Cell{1, minId}, triangulation);
1373 const auto minId = *std::min_element(
1374 sepGeom.begin(), sepGeom.end(),
1375 [&triangulation, offsets,
this](
const Cell &a,
const Cell &b) {
1376 return offsets[discreteGradient_.getCellLowerVertex(a, triangulation)]
1377 < offsets[discreteGradient_.getCellLowerVertex(
1380 sepFuncMin = discreteGradient_.getCellLowerVertex(minId, triangulation);
1382 separatrixFunctionMaxima[sepId] = sepFuncMax;
1383 separatrixFunctionMinima[sepId] = sepFuncMin;
1386 const char onBoundary
1387 = (sepSaddles.empty()
1389 : std::count_if(sepSaddles.begin(), sepSaddles.end(),
1391 return triangulation.isEdgeOnBoundary(a);
1393 + triangulation.isTriangleOnBoundary(src.id_);
1395 for(
size_t j = 0; j < sepGeom.size(); ++j) {
1396 const auto &cell = sepGeom[j];
1400 triangulation.getTriangleVertex(cell.id_, 0, v0);
1401 triangulation.getTriangleVertex(cell.id_, 1, v1);
1402 triangulation.getTriangleVertex(cell.id_, 2, v2);
1405 const auto l = geomCellsBegId[i] + j;
1407 const auto m = l - noldcells;
1409 cellsConn[3 * m + 0] = v0;
1410 cellsConn[3 * m + 1] = v1;
1411 cellsConn[3 * m + 2] = v2;
1412 cellVertsIds[3 * m + 0] = v0;
1413 cellVertsIds[3 * m + 1] = v1;
1414 cellVertsIds[3 * m + 2] = v2;
1425 TTK_PSORT(this->threadNumber_, cellVertsIds.begin(), cellVertsIds.end());
1426 const auto last = std::unique(cellVertsIds.begin(), cellVertsIds.end());
1427 cellVertsIds.erase(last, cellVertsIds.end());
1430 std::vector<size_t> vertId2PointsId(triangulation.getNumberOfVertices());
1432 const auto noldpoints{npoints};
1433 npoints += cellVertsIds.size();
1434 outSeps2.pt.points_.resize(3 * npoints);
1435 auto points = &outSeps2.pt.points_[3 * noldpoints];
1437#ifdef TTK_ENABLE_OPENMP
1438#pragma omp parallel for num_threads(threadNumber_)
1440 for(
size_t i = 0; i < cellVertsIds.size(); ++i) {
1442 triangulation.getVertexPoint(
1443 cellVertsIds[i], points[3 * i + 0], points[3 * i + 1], points[3 * i + 2]);
1445 vertId2PointsId[cellVertsIds[i]] = i + noldpoints;
1448 const auto lastOffset = noldcells == 0 ? 0 : cellsOff[-1];
1450#ifdef TTK_ENABLE_OPENMP
1451#pragma omp parallel for num_threads(threadNumber_)
1453 for(
size_t i = 0; i < ncells - noldcells; ++i) {
1454 cellsOff[i] = 3 * i + lastOffset;
1455 cellsConn[3 * i + 0] = vertId2PointsId[cellsConn[3 * i + 0]];
1456 cellsConn[3 * i + 1] = vertId2PointsId[cellsConn[3 * i + 1]];
1457 cellsConn[3 * i + 2] = vertId2PointsId[cellsConn[3 * i + 2]];
1460 cellsOff[ncells - noldcells] = cellsOff[ncells - noldcells - 1] + 3;
1462 outSeps2.pt.numberOfPoints_ = npoints;
1463 outSeps2.cl.numberOfCells_ = ncells;
1472template <
typename triangulationType>
1474 const std::vector<SimplexId> &maxima,
1476 const triangulationType &triangulation)
const {
1478 if(morseSmaleManifold ==
nullptr) {
1479 this->printErr(
"Could not compute ascending segmentation");
1485 const auto nVerts{triangulation.getNumberOfVertices()};
1486 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, -1);
1487 if(maxima.empty()) {
1492 const auto nCells{triangulation.getNumberOfCells()};
1493 std::vector<SimplexId> morseSmaleManifoldOnCells(nCells, -1);
1496 for(
const auto &
id : maxima) {
1498 morseSmaleManifoldOnCells[id] = nMax++;
1501 const auto dim{triangulation.getDimensionality()};
1504 auto getFaceStarNumber = &triangulationType::getTriangleStarNumber;
1505 auto getFaceStar = &triangulationType::getTriangleStar;
1508 getFaceStarNumber = &triangulationType::getEdgeStarNumber;
1509 getFaceStar = &triangulationType::getEdgeStar;
1510 }
else if(dim == 1) {
1512 getFaceStarNumber = &triangulationType::getVertexStarNumber;
1513 getFaceStar = &triangulationType::getVertexStar;
1517 std::vector<SimplexId> visited{};
1519 std::vector<uint8_t> isMarked(nCells, 0);
1521#ifdef TTK_ENABLE_OPENMP
1522#pragma omp parallel for num_threads(threadNumber_) firstprivate(visited)
1525 if(isMarked[i] == 1) {
1530 while(morseSmaleManifoldOnCells[curr] == -1) {
1531 if(isMarked[curr] == 1) {
1535 const auto paired{this->discreteGradient_.getPairedCell(
1536 Cell{dim, curr}, triangulation,
true)};
1538 const auto nStars{(triangulation.*getFaceStarNumber)(paired)};
1540 (triangulation.*getFaceStar)(paired, j, next);
1546 visited.emplace_back(curr);
1553 for(
const auto el : visited) {
1554 morseSmaleManifoldOnCells[el] = morseSmaleManifoldOnCells[curr];
1559#ifdef TTK_ENABLE_OPENMP
1560#pragma omp parallel for num_threads(threadNumber_)
1563 if(triangulation.getVertexStarNumber(i) < 1) {
1568 triangulation.getVertexStar(i, 0, starId);
1570 morseSmaleManifold[i] = morseSmaleManifoldOnCells[starId];
1573 this->
printMsg(
" Ascending segmentation computed", 1.0, tm.getElapsedTime(),
1580template <
typename triangulationType>
1582 const std::vector<SimplexId> &minima,
1584 const triangulationType &triangulation)
const {
1586 if(morseSmaleManifold ==
nullptr) {
1587 this->printErr(
"Could not compute descending segmentation");
1593 const auto nVerts{triangulation.getNumberOfVertices()};
1595 if(minima.size() == 1) {
1597 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, 0);
1601 std::fill(morseSmaleManifold, morseSmaleManifold + nVerts, -1);
1604 for(
const auto &cp : minima) {
1606 morseSmaleManifold[cp] = nMin++;
1609 std::vector<SimplexId> visited{};
1611#ifdef TTK_ENABLE_OPENMP
1612#pragma omp parallel for num_threads(threadNumber_) firstprivate(visited)
1615 if(morseSmaleManifold[i] != -1) {
1620 while(morseSmaleManifold[curr] == -1) {
1622 const auto pairedEdge{
1623 this->discreteGradient_.getPairedCell(
Cell{0, curr}, triangulation)};
1625 triangulation.getEdgeVertex(pairedEdge, 0, next);
1627 triangulation.getEdgeVertex(pairedEdge, 1, next);
1629 visited.emplace_back(curr);
1632 for(
const auto el : visited) {
1633 morseSmaleManifold[el] = morseSmaleManifold[curr];
1637 this->
printMsg(
" Descending segmentation computed", 1.0, tm.getElapsedTime(),
1644template <
typename triangulationType>
1647 const SimplexId *
const ascendingManifold,
1648 const SimplexId *
const descendingManifold,
1650 const triangulationType &triangulation)
const {
1652 if(ascendingManifold ==
nullptr || descendingManifold ==
nullptr
1653 || morseSmaleManifold ==
nullptr) {
1654 this->printErr(
"Could not compute final segmentation");
1660 const size_t nVerts = triangulation.getNumberOfVertices();
1664#ifdef TTK_ENABLE_OPENMP
1665#pragma omp parallel for num_threads(threadNumber_)
1667 for(
size_t i = 0; i < nVerts; ++i) {
1668 const auto d = ascendingManifold[i];
1669 const auto a = descendingManifold[i];
1670 if(a == -1 || d == -1) {
1671 morseSmaleManifold[i] = -1;
1673 morseSmaleManifold[i] = a * numberOfMaxima + d;
1678 std::vector<SimplexId> sparseRegionIds(
1679 morseSmaleManifold, morseSmaleManifold + nVerts);
1683 this->threadNumber_, sparseRegionIds.begin(), sparseRegionIds.end());
1684 const auto last = std::unique(sparseRegionIds.begin(), sparseRegionIds.end());
1685 sparseRegionIds.erase(last, sparseRegionIds.end());
1688 std::map<SimplexId, size_t> sparseToDenseRegionId{};
1690 for(
size_t i = 0; i < sparseRegionIds.size(); ++i) {
1691 sparseToDenseRegionId[sparseRegionIds[i]] = i;
1696#ifdef TTK_ENABLE_OPENMP
1697#pragma omp parallel for num_threads(threadNumber_)
1699 for(
size_t i = 0; i < nVerts; ++i) {
1700 morseSmaleManifold[i] = sparseToDenseRegionId[morseSmaleManifold[i]];
1703 this->
printMsg(
" Final segmentation computed", 1.0, tm.getElapsedTime(),
1710template <
typename dataType,
typename triangulationType>
1712 const double persistenceThreshold,
1713 const dataType *
const scalars,
1715 const triangulationType &triangulation) {
1719 const auto dim{triangulation.getDimensionality()};
1721 this->printWrn(
"Can't return saddle connectors without a 3D dataset");
1728 dms.setDebugLevel(this->debugLevel_);
1729 dms.setGradient(std::move(this->discreteGradient_));
1733 std::vector<PersPairType> dms_pairs{};
1734 dms.computePersistencePairs(dms_pairs, offsets, triangulation,
false,
true);
1735 this->discreteGradient_ = dms.getGradient();
1737 this->discreteGradient_.setLocalGradient();
1739 const auto getPersistence
1740 = [
this, &triangulation, scalars](
const PersPairType &p) {
1741 return this->discreteGradient_.getPersistence(
1742 Cell{2, p.death},
Cell{1, p.birth}, scalars, triangulation);
1746 const auto firstSadSadPair{std::distance(
1748 std::find_if(dms_pairs.begin(), dms_pairs.end(),
1749 [](
const auto &pair) { return pair.type == 1; }))};
1753 std::vector<bool> isVisited(triangulation.getNumberOfTriangles(),
false);
1754 std::vector<SimplexId> visitedTriangles{};
1759 std::vector<std::tuple<size_t, dataType>> pairs;
1760 for(
size_t i = firstSadSadPair; i < dms_pairs.size(); ++i) {
1761 const auto &pair{dms_pairs[i]};
1762 pairs.emplace_back(std::make_tuple(i, getPersistence(pair)));
1764 const auto comparePersistence
1765 = [](
const std::tuple<size_t, dataType> &pair1,
1766 const std::tuple<size_t, dataType> &pair2) {
1767 return std::get<1>(pair1) < std::get<1>(pair2);
1770 this->threadNumber_, pairs.begin(), pairs.end(), comparePersistence);
1772 std::vector<std::tuple<dataType, SimplexId, SimplexId>> skippedPairsPers;
1775 for(
const auto &pairTup : pairs) {
1776 const auto &pairIndex = std::get<0>(pairTup);
1777 const auto &pair{dms_pairs[pairIndex]};
1778 const auto &pairPersistence = std::get<1>(pairTup);
1780 if(pair.type != 1 || pairPersistence > persistenceThreshold) {
1784 const Cell birth{1, pair.birth};
1785 const Cell death{2, pair.death};
1789 this->discreteGradient_.getDescendingWall(death, mask, triangulation);
1791 std::vector<Cell> vpath{};
1792 bool disableForkReversal = not ForceLoopFreeGradient;
1793 this->discreteGradient_.getAscendingPathThroughWall(
1794 birth, death, isVisited, &vpath, triangulation, disableForkReversal);
1796 if(vpath.back() == death) {
1797 this->discreteGradient_.reverseAscendingPathOnWall(vpath, triangulation);
1801 if(ForceLoopFreeGradient) {
1803 = this->discreteGradient_.detectGradientCycle(birth, triangulation);
1805 std::vector<bool> isVisitedUpdated(
1806 triangulation.getNumberOfTriangles(),
false);
1807 std::vector<SimplexId> visitedTrianglesUpdated{};
1808 VisitedMask maskUpdated{isVisitedUpdated, visitedTrianglesUpdated};
1809 discreteGradient_.getDescendingWall(death, maskUpdated, triangulation);
1810 discreteGradient_.getAscendingPathThroughWall(
1811 birth, death, isVisitedUpdated,
nullptr, triangulation,
false,
true,
1815 this->discreteGradient_.reverseAscendingPathOnWall(
1816 vpath, triangulation,
true);
1817 this->
printMsg(
"Could not return saddle connector " + birth.to_string()
1818 +
" -> " + death.to_string()
1819 +
" without creating a cycle.",
1822 skippedPairsPers.emplace_back(
1823 std::make_tuple(pairPersistence, pair.birth, pair.death));
1828 this->
printMsg(
"Could not return saddle connector " + birth.to_string()
1829 +
" -> " + death.to_string(),
1832 skippedPairsPers.emplace_back(
1833 std::make_tuple(pairPersistence, pair.birth, pair.death));
1839 this->threadNumber_, skippedPairsPers.begin(), skippedPairsPers.end());
1840 for(
unsigned int i = 0; i < skippedPairsPers.size(); ++i)
1848 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 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)
bool ForceLoopFreeGradient
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
std::string to_string(__int128)
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::BOLD+" | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|"+debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream)