diff --git a/.github/workflows/check.yml b/.github/workflows/check.yml index 2b06e9fad7..88d8e30e45 100644 --- a/.github/workflows/check.yml +++ b/.github/workflows/check.yml @@ -15,7 +15,7 @@ jobs: # Check source formatting # # ----------------------- # check-formatting: - runs-on: ubuntu-22.04 + runs-on: ubuntu-26.04 steps: - uses: actions/checkout@v4 diff --git a/core/base/common/welcomeMsg.inl b/core/base/common/welcomeMsg.inl index 109916832e..4ed3e2b1e4 100644 --- a/core/base/common/welcomeMsg.inl +++ b/core/base/common/welcomeMsg.inl @@ -12,7 +12,8 @@ printMsg( debug::LineMode::NEW, stream); printMsg(debug::output::BOLD - + "|_ _|_ _| |/ / / /__\\ \\ |___ \\ / _ \\___ \\ / /_" + + "|_ _|_ _| |/ / / /__\\ \\ |___ \\ / _ " + "\\___ \\ / /_" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, @@ -32,7 +33,8 @@ printMsg( debug::LineMode::NEW, stream); printMsg(debug::output::BOLD - + " |_| |_| |_|\\_\\ | |\\___| | |_____|\\___/_____|\\___/" + + " |_| |_| |_|\\_\\ | |\\___| | " + "|_____|\\___/_____|\\___/" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, diff --git a/core/base/continuousScatterPlot/ContinuousScatterPlot.h b/core/base/continuousScatterPlot/ContinuousScatterPlot.h index 626b1f05e6..536373acee 100644 --- a/core/base/continuousScatterPlot/ContinuousScatterPlot.h +++ b/core/base/continuousScatterPlot/ContinuousScatterPlot.h @@ -499,7 +499,7 @@ int ttk::ContinuousScatterPlot::execute( if(v < 0.0 or (u + v) > 1.0) continue; - // triangle/ray intersection below + // triangle/ray intersection below #ifdef TTK_ENABLE_OPENMP #pragma omp atomic update #endif diff --git a/core/base/contourTree/ContourTree.cpp b/core/base/contourTree/ContourTree.cpp index 4cb9c2ceeb..5a8dc2c118 100644 --- a/core/base/contourTree/ContourTree.cpp +++ b/core/base/contourTree/ContourTree.cpp @@ -2704,17 +2704,23 @@ int ContourTree::computeSkeleton(unsigned int arcResolution) { #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { SubLevelSetTree::computeSkeleton(arcResolution); } + { + SubLevelSetTree::computeSkeleton(arcResolution); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { mergeTree_.computeSkeleton(arcResolution); } + { + mergeTree_.computeSkeleton(arcResolution); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { splitTree_.computeSkeleton(arcResolution); } + { + splitTree_.computeSkeleton(arcResolution); + } } return 0; @@ -2728,17 +2734,23 @@ int ContourTree::smoothSkeleton(unsigned int skeletonSmoothing) { #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { SubLevelSetTree::smoothSkeleton(skeletonSmoothing); } + { + SubLevelSetTree::smoothSkeleton(skeletonSmoothing); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { mergeTree_.smoothSkeleton(skeletonSmoothing); } + { + mergeTree_.smoothSkeleton(skeletonSmoothing); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { splitTree_.smoothSkeleton(skeletonSmoothing); } + { + splitTree_.smoothSkeleton(skeletonSmoothing); + } } return 0; @@ -2752,17 +2764,23 @@ int ContourTree::clearSkeleton() { #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { SubLevelSetTree::clearSkeleton(); } + { + SubLevelSetTree::clearSkeleton(); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { mergeTree_.clearSkeleton(); } + { + mergeTree_.clearSkeleton(); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { splitTree_.clearSkeleton(); } + { + splitTree_.clearSkeleton(); + } } return 0; @@ -2893,12 +2911,16 @@ int ContourTree::simplify(const double &simplificationThreshold, #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { mergeTree_.simplify(simplificationThreshold, metric); } + { + mergeTree_.simplify(simplificationThreshold, metric); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { splitTree_.simplify(simplificationThreshold, metric); } + { + splitTree_.simplify(simplificationThreshold, metric); + } } return 0; diff --git a/core/base/delaunayRipsPersistenceDiagram/geoPHUtils.h b/core/base/delaunayRipsPersistenceDiagram/geoPHUtils.h index e8ca5202e9..87fe37be92 100644 --- a/core/base/delaunayRipsPersistenceDiagram/geoPHUtils.h +++ b/core/base/delaunayRipsPersistenceDiagram/geoPHUtils.h @@ -3,7 +3,7 @@ #include #include -#if((BOOST_VERSION / 100) % 1000) >= 81 +#if ((BOOST_VERSION / 100) % 1000) >= 81 #include #include #else @@ -11,7 +11,7 @@ #include #endif -#if((BOOST_VERSION / 100) % 1000) >= 84 +#if ((BOOST_VERSION / 100) % 1000) >= 84 #include #define TTK_CONCURRENT_HASHTABLE_AVAILABLE #endif @@ -47,7 +47,7 @@ namespace ttk::gph { template using PointCloud = std::vector>; -#if((BOOST_VERSION / 100) % 1000) >= 81 +#if ((BOOST_VERSION / 100) % 1000) >= 81 template using HashMap = boost::unordered_flat_map; template @@ -59,7 +59,7 @@ namespace ttk::gph { using HashSet = boost::unordered_set; #endif -#if((BOOST_VERSION / 100) % 1000) >= 84 +#if ((BOOST_VERSION / 100) % 1000) >= 84 template using ConcurrentHashMap = boost::concurrent_flat_map; #endif diff --git a/core/base/depthImageBasedGeometryApproximation/DepthImageBasedGeometryApproximation.h b/core/base/depthImageBasedGeometryApproximation/DepthImageBasedGeometryApproximation.h index 2c881db4e7..082f934f70 100644 --- a/core/base/depthImageBasedGeometryApproximation/DepthImageBasedGeometryApproximation.h +++ b/core/base/depthImageBasedGeometryApproximation/DepthImageBasedGeometryApproximation.h @@ -253,16 +253,16 @@ int ttk::DepthImageBasedGeometryApproximation::execute( triangleDistortions[triangleDistortionOffset++] = isNaN(i0Depth) || isNaN(i2Depth) || isNaN(i1Depth) ? myNan - : std::max( - absDiff(i0Depth, i1Depth), - std::max(absDiff(i1Depth, i2Depth), absDiff(i0Depth, i2Depth))); + : std::max(absDiff(i0Depth, i1Depth), + std::max(absDiff(i1Depth, i2Depth), + absDiff(i0Depth, i2Depth))); triangleDistortions[triangleDistortionOffset++] = isNaN(i1Depth) || isNaN(i2Depth) || isNaN(i3Depth) ? myNan - : std::max( - absDiff(i1Depth, i3Depth), - std::max(absDiff(i3Depth, i2Depth), absDiff(i2Depth, i1Depth))); + : std::max(absDiff(i1Depth, i3Depth), + std::max(absDiff(i3Depth, i2Depth), + absDiff(i2Depth, i1Depth))); } } } diff --git a/core/base/discreteGradient/DiscreteGradient.h b/core/base/discreteGradient/DiscreteGradient.h index bf653879d1..0fe3984748 100644 --- a/core/base/discreteGradient/DiscreteGradient.h +++ b/core/base/discreteGradient/DiscreteGradient.h @@ -428,6 +428,14 @@ user in the gradient. const triangulationType &triangulation, const bool enableCycleDetector = false) const; + /** + * Return all VPath coming from the given cell. + */ + template + int getAllAscendingPaths(const Cell &cell, + std::vector> &vpaths, + const triangulationType &triangulation) const; + /** * Return the VPath terminating at the given cell. */ @@ -436,6 +444,14 @@ user in the gradient. std::vector &vpath, const triangulationType &triangulation) const; + /** + * Return all VPath terminating at the given cell. + */ + template + int getAllDescendingPaths(const Cell &cell, + std::vector> &vpaths, + const triangulationType &triangulation) const; + /** * Return the VPath terminating at the given 2-saddle restricted to the 2-separatrice of the 1-saddle. diff --git a/core/base/discreteGradient/DiscreteGradient_Template.h b/core/base/discreteGradient/DiscreteGradient_Template.h index 25dbd7b7c1..57553fb3b0 100644 --- a/core/base/discreteGradient/DiscreteGradient_Template.h +++ b/core/base/discreteGradient/DiscreteGradient_Template.h @@ -1548,6 +1548,124 @@ int DiscreteGradient::getDescendingPath( } } while(connectedEdgeId != -1); + } else { + printWrn("Descending path not implemented for this simplex dimension!"); + } + + return 0; +} + +template +int DiscreteGradient::getAllDescendingPaths( + const Cell &cell, + std::vector> &vpaths, + const triangulationType &triangulation) const { + + vpaths.clear(); + + using vPath = std::vector; + + // Each stack entry carries the current partial path and the current simplex. + // We use DFS to enumerate all paths (branching is possible at each step). + struct StackEntry { + vPath partialPath_; + Cell currentCell_; + }; + + std::stack stack; + + // Bootstrap: push the starting simplex onto the stack + { + StackEntry stackEntry; + stackEntry.currentCell_ = cell; + stackEntry.partialPath_.push_back(cell); + stack.push(std::move(stackEntry)); + } + + while(!stack.empty()) { + + StackEntry stackEntry = std::move(stack.top()); + stack.pop(); + + const Cell ¤tCell = stackEntry.currentCell_; + + // 1. Follow the gradient arrow out of `currentCell` + const SimplexId pairedId = getPairedCell(currentCell, triangulation); + + if(pairedId == -1) { + // currentCell is a critical simplex: this path has terminated. + vpaths.push_back(stackEntry.partialPath_); + continue; + } + + // The gradient arrow takes us to a (dim+1)-simplex + const int pairedDim = currentCell.dim_ + 1; + Cell paired; + paired.dim_ = pairedDim; + paired.id_ = pairedId; + + // Append the paired simplex to the path + stackEntry.partialPath_.push_back(paired); + + // 2. Enumerate all facets of `paired` of dimension `dim` + // The discrete gradient arrow *entering* `paired` came from `currentCell`. + // We continue the path by following gradient arrows out of the OTHER + // facets of `paired` (i.e., cofacets of `paired` in dimension dim + // that are themselves paired to a simplex of dimension dim+1, + // or terminate if critical). + // + // Standard discrete Morse theory: we look at all dim-faces of `paired`, + // exclude `currentCell` itself, and for each remaining face that is paired + // (i.e., its gradient arrow points to some (dim+1)-simplex, not back + // to `paired`), we spawn a new path branch. + + SimplexId numFacets = 0; + if(pairedDim == 1) { + numFacets = 2; + } else if(pairedDim == 2) { + numFacets = 3; + } else if(pairedDim == 3) { + numFacets = 4; + } + + bool anyBranch = false; + + for(SimplexId f = 0; f < numFacets; ++f) { + SimplexId facetId = -1; + + if(pairedDim == 1) { + triangulation.getEdgeVertex(pairedId, f, facetId); + } else if(pairedDim == 2) { + triangulation.getTriangleEdge(pairedId, f, facetId); + } else if(pairedDim == 3) { + triangulation.getCellTriangle(pairedId, f, facetId); + } + + if(facetId == -1) + continue; + + // Skip the facet we just came from + if(facetId == currentCell.id_) + continue; + + // This facet is a dim-simplex; start a new branch of the V-path from it + Cell nextSimplex; + nextSimplex.dim_ = currentCell.dim_; + nextSimplex.id_ = facetId; + + StackEntry newStackEntry; + // copy currentCellentCell_ path + newStackEntry.partialPath_ = stackEntry.partialPath_; + newStackEntry.partialPath_.push_back(nextSimplex); + newStackEntry.currentCell_ = nextSimplex; + stack.push(std::move(newStackEntry)); + anyBranch = true; + } + + // If no other facet was found (degenerate case), terminate the path here + if(!anyBranch) { + vpaths.push_back(stackEntry.partialPath_); + } } return 0; @@ -1719,6 +1837,8 @@ int DiscreteGradient::getAscendingPath(const Cell &cell, // stop at convergence caused by boundary effect } while(currentId != oldId); + } else { + printWrn("Ascending path not implemented for this simplex dimension!"); } } else if(dimensionality_ == 3) { if(cell.dim_ == 3) { @@ -1780,6 +1900,155 @@ int DiscreteGradient::getAscendingPath(const Cell &cell, // stop at convergence caused by boundary effect } while(currentId != oldId); + } else { + printWrn("Ascending path not implemented for this simplex dimension!"); + } + } else { + printWrn("Ascending path not implemented for this input dimension!"); + } + + return 0; +} + +template +int DiscreteGradient::getAllAscendingPaths( + const Cell &cell, + std::vector> &vpaths, + const triangulationType &triangulation) const { + + vpaths.clear(); + + using vPath = std::vector; + + struct StackEntry { + vPath partialPath_; + Cell currentCell_; + }; + + std::stack stack; + + { + StackEntry stackEntry; + stackEntry.currentCell_ = cell; + stackEntry.partialPath_.push_back(cell); + stack.push(std::move(stackEntry)); + } + + while(!stack.empty()) { + + StackEntry stackEntry = std::move(stack.top()); + stack.pop(); + + const Cell ¤tCell = stackEntry.currentCell_; + const SimplexId pairedCofacetId = getPairedCell(currentCell, triangulation); + + if((currentCell.id_ != cell.id_) && (isCellCritical(currentCell))) { + // currentCell is a critical simplex: this path has terminated. + // the simplex has already been added to the stack path + vpaths.push_back(stackEntry.partialPath_); + continue; + } + + if(currentCell.dim_ != cell.dim_) + continue; + + if(currentCell.dim_ >= dimensionality_) { + // currentCell is a maximal simplex: it admits no cofacet at all (in + // particular, the star of a triangle is only defined in 3D). the + // ascending path terminates here. + vpaths.push_back(stackEntry.partialPath_); + continue; + } + + // check all cofacets + SimplexId cofacetNumber = -1; + + switch(currentCell.dim_) { + case 1: + cofacetNumber = triangulation.getEdgeTriangleNumber(currentCell.id_); + break; + case 2: + cofacetNumber = triangulation.getTriangleStarNumber(currentCell.id_); + break; + default: + cofacetNumber = triangulation.getVertexEdgeNumber(currentCell.id_); + break; + } + + bool hasProgressed = false; + + for(SimplexId i = 0; i < cofacetNumber; i++) { + SimplexId cofacetId = -1; + switch(currentCell.dim_) { + case 1: + triangulation.getEdgeTriangle(currentCell.id_, i, cofacetId); + break; + case 2: + triangulation.getTriangleStar(currentCell.id_, i, cofacetId); + break; + default: + triangulation.getVertexEdge(currentCell.id_, i, cofacetId); + break; + } + if(cofacetId != pairedCofacetId) { + + // we don't want to go down the v-path, we want to go backwards + Cell cofacet; + cofacet.dim_ = currentCell.dim_ + 1; + cofacet.id_ = cofacetId; + + // the path shared by all the branches below (one per face of the + // cofacet): each branch extends its own copy of it + std::vector cofacetPath = stackEntry.partialPath_; + cofacetPath.push_back(cofacet); + + // now find the simplex we came from + int simplexNumber = -1; + + simplexNumber = cofacet.dim_ + 1; + + for(int j = 0; j < simplexNumber; j++) { + SimplexId simplexId = -1; + switch(cofacet.dim_) { + case 1: + triangulation.getEdgeVertex(cofacet.id_, j, simplexId); + break; + case 2: + triangulation.getTriangleEdge(cofacet.id_, j, simplexId); + break; + default: + triangulation.getCellTriangle(cofacet.id_, j, simplexId); + break; + } + + Cell simplex; + simplex.id_ = simplexId; + simplex.dim_ = cofacet.dim_ - 1; + const SimplexId simplexPair = getPairedCell(simplex, triangulation); + + if(isCellCritical(simplex)) { + // always terminate here — don't continue the path through a + // critical cell + std::vector criticalPath = cofacetPath; + criticalPath.push_back(simplex); + vpaths.push_back(std::move(criticalPath)); + hasProgressed = true; // prevent the fallback push too + // do NOT push to stack + } else if(simplexPair == cofacet.id_) { + StackEntry newStackEntry; + newStackEntry.partialPath_ = cofacetPath; + newStackEntry.partialPath_.push_back(simplex); + newStackEntry.currentCell_ = simplex; + stack.push(std::move(newStackEntry)); + hasProgressed = true; + } + } + } + } + if(!hasProgressed) { + // example: boundary edge paired with its interior cofacet, we stop the + // backward vpath here. + vpaths.push_back(stackEntry.partialPath_); } } diff --git a/core/base/discreteMorseSandwich/DiscreteMorseSandwich.cpp b/core/base/discreteMorseSandwich/DiscreteMorseSandwich.cpp index ed0499fb2b..f136add1f4 100644 --- a/core/base/discreteMorseSandwich/DiscreteMorseSandwich.cpp +++ b/core/base/discreteMorseSandwich/DiscreteMorseSandwich.cpp @@ -131,9 +131,10 @@ void ttk::DiscreteMorseSandwich::displayStats( std::count_if(pairs.begin(), pairs.end(), [](const PersistencePair &a) { return a.type == 0; }))}, {" #Saddle-saddle pairs", - std::to_string(dim == 3 ? std::count_if( - pairs.begin(), pairs.end(), - [](const PersistencePair &a) { return a.type == 1; }) + std::to_string(dim == 3 ? std::count_if(pairs.begin(), pairs.end(), + [](const PersistencePair &a) { + return a.type == 1; + }) : 0)}, {" #Saddle-max pairs", std::to_string(std::count_if( diff --git a/core/base/discreteMorseSandwich/DiscreteMorseSandwich.h b/core/base/discreteMorseSandwich/DiscreteMorseSandwich.h index 93a5584c1b..a0370575e5 100644 --- a/core/base/discreteMorseSandwich/DiscreteMorseSandwich.h +++ b/core/base/discreteMorseSandwich/DiscreteMorseSandwich.h @@ -699,31 +699,32 @@ void ttk::DiscreteMorseSandwich::getMaxSaddlePairs( const auto dim = this->dg_.getDimensionality(); auto saddle2ToMaxima - = dim == 3 - ? getSaddle2ToMaxima( - criticalSaddles, - [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) { - return triangulation.getTriangleStar(a, i, r); - }, - [&triangulation](const SimplexId a) { - return triangulation.getTriangleStarNumber(a); - }, - [&triangulation](const SimplexId a) { - return triangulation.isTriangleOnBoundary(a); - }, - triangulation) - : getSaddle2ToMaxima( - criticalSaddles, - [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) { - return triangulation.getEdgeStar(a, i, r); - }, - [&triangulation](const SimplexId a) { - return triangulation.getEdgeStarNumber(a); - }, - [&triangulation](const SimplexId a) { - return triangulation.isEdgeOnBoundary(a); - }, - triangulation); + = dim == 3 ? getSaddle2ToMaxima( + criticalSaddles, + [&triangulation]( + const SimplexId a, const SimplexId i, SimplexId &r) { + return triangulation.getTriangleStar(a, i, r); + }, + [&triangulation](const SimplexId a) { + return triangulation.getTriangleStarNumber(a); + }, + [&triangulation](const SimplexId a) { + return triangulation.isTriangleOnBoundary(a); + }, + triangulation) + : getSaddle2ToMaxima( + criticalSaddles, + [&triangulation]( + const SimplexId a, const SimplexId i, SimplexId &r) { + return triangulation.getEdgeStar(a, i, r); + }, + [&triangulation](const SimplexId a) { + return triangulation.getEdgeStarNumber(a); + }, + [&triangulation](const SimplexId a) { + return triangulation.isEdgeOnBoundary(a); + }, + triangulation); Timer tmseq{}; diff --git a/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.cpp b/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.cpp index afb7fbdb83..da2c9fc218 100644 --- a/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.cpp +++ b/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.cpp @@ -31,9 +31,10 @@ void ttk::DiscreteMorseSandwichMPI::displayStats( std::count_if(pairs.begin(), pairs.end(), [](const PersistencePair &a) { return a.type == 0; }))}, {" #Saddle-saddle pairs", - std::to_string(dim == 3 ? std::count_if( - pairs.begin(), pairs.end(), - [](const PersistencePair &a) { return a.type == 1; }) + std::to_string(dim == 3 ? std::count_if(pairs.begin(), pairs.end(), + [](const PersistencePair &a) { + return a.type == 1; + }) : 0)}, {" #Saddle-max pairs", std::to_string(std::count_if( diff --git a/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.h b/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.h index 038e912a0f..14d8d6cd9b 100644 --- a/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.h +++ b/core/base/discreteMorseSandwichMPI/DiscreteMorseSandwichMPI.h @@ -2197,8 +2197,8 @@ int ttk::DiscreteMorseSandwichMPI::getSaddle1ToMinima( }; ttk::SimplexId elementNumber = 0; // follow vpaths from 1-saddles to minima -#pragma omp parallel shared(extremaLocks, saddleAtomic) reduction(+: elementNumber) \ - num_threads(localThreadNumber) +#pragma omp parallel shared(extremaLocks, saddleAtomic) \ + reduction(+ : elementNumber) num_threads(localThreadNumber) { int threadNumber = omp_get_thread_num(); #pragma omp for schedule(static) @@ -2428,7 +2428,8 @@ void ttk::DiscreteMorseSandwichMPI::getSaddle2ToMaxima( ttk::SimplexId totalFinishedElement{0}; ttk::SimplexId totalElement{0}; #ifdef TTK_ENABLE_OPENMP -#pragma omp parallel for num_threads(localThreadNumber) reduction(+:totalFinishedElement) +#pragma omp parallel for num_threads(localThreadNumber) \ + reduction(+ : totalFinishedElement) #endif for(size_t i = 0; i < criticalSaddles.size(); ++i) { totalFinishedElement += getFaceStarNumber(criticalSaddles[i]); @@ -2519,8 +2520,8 @@ void ttk::DiscreteMorseSandwichMPI::getSaddle2ToMaxima( }; // follow vpaths from 2-saddles to maxima char saddleLocalId; -#pragma omp parallel shared(extremaLocks, saddleAtomic) reduction(+: elementNumber) \ - num_threads(localThreadNumber) +#pragma omp parallel shared(extremaLocks, saddleAtomic) \ + reduction(+ : elementNumber) num_threads(localThreadNumber) { int threadNumber = omp_get_thread_num(); #pragma omp for schedule(static) @@ -2767,7 +2768,9 @@ void ttk::DiscreteMorseSandwichMPI::unpackGhostPresence( // Add the entry to the map if(lid == -1 || getSimplexRank(lid) != ttk::MPIrank_) { #pragma omp critical - { localGhostPresenceMap[vp.extremaId_] = ghost; } + { + localGhostPresenceMap[vp.extremaId_] = ghost; + } } else { lid = localTriangToLocalVectExtrema.find(lid)->second; extremaLocks[lid].lock(); @@ -2867,9 +2870,11 @@ void ttk::DiscreteMorseSandwichMPI::getMinSaddlePairs( if(criticalExtremasNumber > 0) { // extracts the global min -#pragma omp declare reduction(get_min : std::pair :omp_out = omp_out.second < omp_in.second ? omp_out : omp_in) -#pragma omp parallel for reduction(get_min \ - : localMin) num_threads(localThreadNumber) +#pragma omp declare reduction( \ + get_min : std::pair : omp_out \ + = omp_out.second < omp_in.second ? omp_out : omp_in) +#pragma omp parallel for reduction(get_min : localMin) \ + num_threads(localThreadNumber) for(ttk::SimplexId i = 0; i < criticalExtremasNumber; i++) { if(offsets[criticalExtremas[i]] < localMin.second) { localMin.first = criticalExtremas[i]; @@ -2925,9 +2930,10 @@ void ttk::DiscreteMorseSandwichMPI::getMinSaddlePairs( MPI_IN_PLACE, &totalNumberOfPairs, 1, MPI_SimplexId, MPI_SUM, MPIcomm); globalToLocalSaddle.reserve(criticalEdgesNumber); -#pragma omp declare reduction (merge : std::vector: omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) -#pragma omp parallel for reduction(merge \ - : extremasGid) schedule(static) \ +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp parallel for reduction(merge : extremasGid) schedule(static) \ shared(saddles) num_threads(localThreadNumber) for(ttk::SimplexId i = 0; i < criticalEdgesNumber; ++i) { auto &mins = saddle1ToMinima[i]; @@ -2954,8 +2960,8 @@ void ttk::DiscreteMorseSandwichMPI::getMinSaddlePairs( extremasGid.size(), std::vector()); std::vector extremaLocks(extremasGid.size(), 0); std::vector> extremas(extremasGid.size(), extremaNode<1>()); -#pragma omp parallel master shared(extremaLocks, extremas, globalMinLid, \ - ghostPresence, saddles, globalMinOffset) \ +#pragma omp parallel master shared(extremaLocks, extremas, globalMinLid, \ + ghostPresence, saddles, globalMinOffset) \ num_threads(localThreadNumber) { #pragma omp task @@ -3204,9 +3210,10 @@ void ttk::DiscreteMorseSandwichMPI::computeMaxSaddlePairs( MPI_Allreduce( MPI_IN_PLACE, &totalNumberOfPairs, 1, MPI_SimplexId, MPI_SUM, MPIcomm); globalToLocalSaddle.reserve(criticalSaddlesNumber); -#pragma omp declare reduction (merge : std::vector: omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) -#pragma omp parallel for reduction(merge \ - : extremasGid) schedule(static) \ +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp parallel for reduction(merge : extremasGid) schedule(static) \ shared(saddles, saddle2ToMaxima) num_threads(localThreadNumber) for(ttk::SimplexId i = 0; i < criticalSaddlesNumber; ++i) { auto &maxs = saddle2ToMaxima[i]; @@ -3254,7 +3261,7 @@ void ttk::DiscreteMorseSandwichMPI::computeMaxSaddlePairs( std::vector> extremas( extremasGid.size(), extremaNode()); #pragma omp parallel master shared(extremaLocks, extremas, ghostPresence, \ - saddles) num_threads(localThreadNumber) + saddles) num_threads(localThreadNumber) { #pragma omp task { @@ -3409,8 +3416,7 @@ void ttk::DiscreteMorseSandwichMPI::getMaxSaddlePairs( ttk::SimplexId globalMaxOffset{0}; if(criticalExtremasNumber > 0) { // extracts the global max -#pragma omp parallel for reduction(max \ - : globalMaxOffset) \ +#pragma omp parallel for reduction(max : globalMaxOffset) \ num_threads(localThreadNumber) for(ttk::SimplexId i = 0; i < vertexNumber; i++) { if(globalMaxOffset < offsets[i]) { @@ -3422,9 +3428,10 @@ void ttk::DiscreteMorseSandwichMPI::getMaxSaddlePairs( MPI_Allreduce( MPI_IN_PLACE, &globalMaxOffset, 1, MPI_SimplexId, MPI_MAX, MPIcomm); -#pragma omp declare reduction (merge : std::vector : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) -#pragma omp parallel for reduction(merge \ - : localMaxId) \ +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp parallel for reduction(merge : localMaxId) \ num_threads(localThreadNumber) for(ttk::SimplexId i = 0; i < criticalExtremasNumber; i++) { if(triangulation.getCellRank(criticalExtremas[i]) == ttk::MPIrank_) { @@ -3840,9 +3847,10 @@ void ttk::DiscreteMorseSandwichMPI::extractPairs( ttk::SimplexId saddleNumber = saddleToPairedExtrema.size(); #ifdef TTK_ENABLE_OPENMP -#pragma omp declare reduction (merge : std::vector : omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) -#pragma omp parallel for reduction(merge \ - : pairs) schedule(static) \ +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp parallel for reduction(merge : pairs) schedule(static) \ num_threads(localThreadNumber) #endif for(ttk::SimplexId i = 0; i < saddleNumber; i++) { @@ -3870,7 +3878,8 @@ ttk::SimplexId ttk::DiscreteMorseSandwichMPI::computePairNumbers( ttk::SimplexId saddleNumber = saddleToPairedExtrema.size(); ttk::SimplexId computedSaddleNumber{0}; #ifdef TTK_ENABLE_OPENMP -#pragma omp parallel for reduction(+ : computedSaddleNumber) schedule(static) num_threads(localThreadNumber) +#pragma omp parallel for reduction(+ : computedSaddleNumber) schedule(static) \ + num_threads(localThreadNumber) #endif for(ttk::SimplexId i = 0; i < saddleNumber; i++) { if(saddleToPairedExtrema[i] > -1 && saddles[i].rank_ == ttk::MPIrank_) { @@ -5710,7 +5719,9 @@ void ttk::DiscreteMorseSandwichMPI::getSaddleSaddlePairs( ttk::Timer t_mpi; ttk::startMPITimer(t_mpi, ttk::MPIrank_, ttk::MPIsize_); #endif -#pragma omp declare reduction (merge : std::vector: omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) #pragma omp parallel for reduction(merge : saddles1Gid) schedule(static) for(size_t i = 0; i < critical1Saddles.size(); i++) { const auto s1 = critical1Saddles[i]; @@ -5838,9 +5849,9 @@ void ttk::DiscreteMorseSandwichMPI::getSaddleSaddlePairs( = std::min(saddle2Number + 1, static_cast(10)); ttk::SimplexId taskNum = static_cast(saddle2Number / taskSize) + 1; -#pragma omp parallel num_threads(threadNumber_) shared( \ - onBoundaryThread, s1Locks, s2Locks, s2GlobalBoundaries, s2LocalBoundaries, \ - localEdgeToSaddle1_, saddles2, edgeTrianglePartner) +#pragma omp parallel num_threads(threadNumber_) shared( \ + onBoundaryThread, s1Locks, s2Locks, s2GlobalBoundaries, s2LocalBoundaries, \ + localEdgeToSaddle1_, saddles2, edgeTrianglePartner) { #pragma omp single nowait { @@ -6014,7 +6025,7 @@ void ttk::DiscreteMorseSandwichMPI::getSaddleSaddlePairs( ->second; #pragma omp task firstprivate(lid) \ shared(s2GlobalBoundaries, s2LocalBoundaries, edgeTrianglePartner, s1Locks, \ - s2Locks, saddles2) + s2Locks, saddles2) { ttk::SimplexId lidBlock; ttk::SimplexId lidElement; @@ -6052,7 +6063,9 @@ void ttk::DiscreteMorseSandwichMPI::getSaddleSaddlePairs( Timer tmseq{}; // extract saddle-saddle pairs from computed boundaries -#pragma omp declare reduction (merge : std::vector: omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end())) +#pragma omp declare reduction( \ + merge : std::vector : omp_out.insert( \ + omp_out.end(), omp_in.begin(), omp_in.end())) #pragma omp parallel for reduction(merge : pairs) schedule(static) for(size_t i = 0; i < edgeTrianglePartner.size(); ++i) { if(edgeTrianglePartner[i] != -1) { diff --git a/core/base/discreteVectorField/DiscreteVectorField.h b/core/base/discreteVectorField/DiscreteVectorField.h index 25d78b87f2..510df06be9 100644 --- a/core/base/discreteVectorField/DiscreteVectorField.h +++ b/core/base/discreteVectorField/DiscreteVectorField.h @@ -60,8 +60,8 @@ namespace ttk { const std::array &lowVerts, const std::array &lowVertWeights, const std::array &faces) - : Cell{dim, id}, lowVerts_{lowVerts}, - lowVertWeights_{lowVertWeights}, faces_{faces} { + : Cell{dim, id}, lowVerts_{lowVerts}, lowVertWeights_{lowVertWeights}, + faces_{faces} { } // ID values for Outward vertices in current Outward star diff --git a/core/base/distanceMatrixDistortion/DistanceMatrixDistortion.cpp b/core/base/distanceMatrixDistortion/DistanceMatrixDistortion.cpp index e5e1fb63ac..e91c515987 100644 --- a/core/base/distanceMatrixDistortion/DistanceMatrixDistortion.cpp +++ b/core/base/distanceMatrixDistortion/DistanceMatrixDistortion.cpp @@ -41,9 +41,8 @@ int ttk::DistanceMatrixDistortion::execute( } #ifdef TTK_ENABLE_OPENMP -#pragma omp parallel for num_threads(this->threadNumber_) reduction(max \ - : maxi) \ - schedule(dynamic) +#pragma omp parallel for num_threads(this->threadNumber_) \ + reduction(max : maxi) schedule(dynamic) #endif // TTK_ENABLE_OPENMP for(size_t i = 0; i < n; i++) { for(size_t j = i + 1; j < n; j++) { diff --git a/core/base/ftmTree/FTMAtomicVector.h b/core/base/ftmTree/FTMAtomicVector.h index 80477a73a3..59e7df268d 100644 --- a/core/base/ftmTree/FTMAtomicVector.h +++ b/core/base/ftmTree/FTMAtomicVector.h @@ -69,7 +69,9 @@ namespace ttk { // WARNING: In parallel we do not want to make reserve as it can lead // to data race, we should not enter here #pragma omp critical(AtomicUFReserve) - { std::vector::resize(newSize, defaultValue); } + { + std::vector::resize(newSize, defaultValue); + } } else #endif diff --git a/core/base/ftmTree/FTMTree_CT_Template.h b/core/base/ftmTree/FTMTree_CT_Template.h index 2ff214b909..11c0d372bb 100644 --- a/core/base/ftmTree/FTMTree_CT_Template.h +++ b/core/base/ftmTree/FTMTree_CT_Template.h @@ -24,7 +24,9 @@ namespace ttk { #ifdef TTK_ENABLE_OPENMP4 #pragma omp single nowait #endif - { leafSearch(mesh); } + { + leafSearch(mesh); + } } printTime(precomputeTime, "leafSearch", 3); } @@ -100,61 +102,61 @@ namespace ttk { this->printMsg({"- final number of nodes :", nbNodes}); } } -// clang-format on -// clang format fail to use the right indentation level -// here, but it break the code if not disabled... + // clang-format on + // clang format fail to use the right indentation level + // here, but it break the code if not disabled... -// ------------------------------------------------------------------------ + // ------------------------------------------------------------------------ -template -int FTMTree_CT::leafSearch(const triangulationType *mesh) { - const auto nbScalars = scalars_->size; - const auto chunkSize = getChunkSize(); - const auto chunkNb = getChunkCount(); + template + int FTMTree_CT::leafSearch(const triangulationType *mesh) { + const auto nbScalars = scalars_->size; + const auto chunkSize = getChunkSize(); + const auto chunkNb = getChunkCount(); - // Extrema extract and launch tasks - for(SimplexId chunkId = 0; chunkId < chunkNb; ++chunkId) { + // Extrema extract and launch tasks + for(SimplexId chunkId = 0; chunkId < chunkNb; ++chunkId) { #ifdef TTK_ENABLE_OPENMP4 #pragma omp task firstprivate(chunkId) #endif - { - const SimplexId lowerBound = chunkId * chunkSize; - const SimplexId upperBound - = std::min(nbScalars, (chunkId + 1) * chunkSize); - for(SimplexId v = lowerBound; v < upperBound; ++v) { - const auto &neighNumb = mesh->getVertexNeighborNumber(v); - valence upval = 0; - valence downval = 0; - - for(valence n = 0; n < neighNumb; ++n) { - SimplexId neigh{-1}; - mesh->getVertexNeighbor(v, n, neigh); - if(scalars_->isLower(neigh, v)) { - ++downval; - } else { - ++upval; + { + const SimplexId lowerBound = chunkId * chunkSize; + const SimplexId upperBound + = std::min(nbScalars, (chunkId + 1) * chunkSize); + for(SimplexId v = lowerBound; v < upperBound; ++v) { + const auto &neighNumb = mesh->getVertexNeighborNumber(v); + valence upval = 0; + valence downval = 0; + + for(valence n = 0; n < neighNumb; ++n) { + SimplexId neigh{-1}; + mesh->getVertexNeighbor(v, n, neigh); + if(scalars_->isLower(neigh, v)) { + ++downval; + } else { + ++upval; + } + } + + jt_.setValence(v, downval); + st_.setValence(v, upval); + + if(!downval) { + jt_.makeNode(v); + } + + if(!upval) { + st_.makeNode(v); + } } } - - jt_.setValence(v, downval); - st_.setValence(v, upval); - - if(!downval) { - jt_.makeNode(v); - } - - if(!upval) { - st_.makeNode(v); - } } - } - } #ifdef TTK_ENABLE_OPENMP4 #pragma omp taskwait #endif - return 0; -} + return 0; + } -} // namespace ftm + } // namespace ftm } // namespace ttk diff --git a/core/base/ftmTree/FTMTree_MT.cpp b/core/base/ftmTree/FTMTree_MT.cpp index 4a3b355799..db7d8a3dd8 100644 --- a/core/base/ftmTree/FTMTree_MT.cpp +++ b/core/base/ftmTree/FTMTree_MT.cpp @@ -918,7 +918,9 @@ vector FTMTree_MT::sortedNodes(const bool para) { #ifdef TTK_ENABLE_OPENMP #pragma omp single #endif - { std::sort(sortedNodes.begin(), sortedNodes.end(), indirect_sort); } + { + std::sort(sortedNodes.begin(), sortedNodes.end(), indirect_sort); + } } return sortedNodes; diff --git a/core/base/ftrGraph/FTRGraphPrivate_Template.h b/core/base/ftrGraph/FTRGraphPrivate_Template.h index a117c8b0bd..2eeba2f0f8 100644 --- a/core/base/ftrGraph/FTRGraphPrivate_Template.h +++ b/core/base/ftrGraph/FTRGraphPrivate_Template.h @@ -38,7 +38,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif @@ -72,7 +74,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif @@ -214,7 +218,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif @@ -244,7 +250,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif @@ -336,7 +344,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif @@ -382,7 +392,9 @@ void ttk::ftr::FTRGraph::growthFromSeed( --nbProp_; } #pragma omp critical(stats) - { curTime = sweepStart_.getElapsedTime(); } + { + curTime = sweepStart_.getElapsedTime(); + } propTimes_[curProp - 1] = curTime; } #endif diff --git a/core/base/geometry/Geometry.cpp b/core/base/geometry/Geometry.cpp index eaa1e6dd1c..2dd1e719ce 100644 --- a/core/base/geometry/Geometry.cpp +++ b/core/base/geometry/Geometry.cpp @@ -534,6 +534,34 @@ T Geometry::magnitude(const T *o, const T *d) { return sqrt(mag); } +template +int Geometry::normalizeBarycentricWeights(std::vector &baryCentrics) { + + if(baryCentrics.empty()) + return -1; + + T sum = 0; + + for(size_t i = 0; i < baryCentrics.size(); i++) { + // clamp the negative weights induced by numerical inaccuracies + if(baryCentrics[i] < 0) + baryCentrics[i] = 0; + sum += baryCentrics[i]; + } + + if(!(sum > 0)) { + // degenerated weights: fall back on the barycenter + for(size_t i = 0; i < baryCentrics.size(); i++) + baryCentrics[i] = 1.0 / baryCentrics.size(); + return -2; + } + + for(size_t i = 0; i < baryCentrics.size(); i++) + baryCentrics[i] /= sum; + + return 0; +} + template void Geometry::projectOnTrianglePlane(const T *p, const T *a, @@ -885,6 +913,8 @@ void Geometry::transposeMatrix(const std::vector> &a, template TYPE Geometry::magnitudeFlatten( \ std::vector> const &); \ template TYPE Geometry::magnitude(TYPE const *, TYPE const *); \ + template int Geometry::normalizeBarycentricWeights( \ + std::vector &); \ template void Geometry::projectOnTrianglePlane( \ TYPE const *, TYPE const *, TYPE const *, TYPE *); \ template void Geometry::projectOnEdge( \ diff --git a/core/base/geometry/Geometry.h b/core/base/geometry/Geometry.h index 346630536d..125f63b9df 100644 --- a/core/base/geometry/Geometry.h +++ b/core/base/geometry/Geometry.h @@ -377,6 +377,16 @@ namespace ttk { template T magnitude(const T *o, const T *d); + /// Normalize a list of barycentric weights: the negative weights + /// (induced by numerical inaccuracies) are clamped to zero and the + /// remaining ones are re-scaled, such that they sum up to one. + /// \param baryCentrics Input/output barycentric weights. + /// \return Returns 0 upon success, negative values otherwise (in + /// particular, -2 if the input weights are degenerated, in which case the + /// output weights are those of the barycenter). + template + int normalizeBarycentricWeights(std::vector &baryCentrics); + /// Compute the integer power of a floating-point value /// (std::pow is optimised for floating-point exponents) template diff --git a/core/base/implicitTriangulation/ImplicitTriangulation.h b/core/base/implicitTriangulation/ImplicitTriangulation.h index 2824b66cad..eda8e9be3d 100644 --- a/core/base/implicitTriangulation/ImplicitTriangulation.h +++ b/core/base/implicitTriangulation/ImplicitTriangulation.h @@ -110,29 +110,33 @@ namespace ttk { virtual int getTetrahedronEdge(const SimplexId &tetId, const int &id, - SimplexId &edgeId) const = 0; + SimplexId &edgeId) const + = 0; int getTetrahedronEdges(std::vector> &edges) const; virtual int getTetrahedronTriangle(const SimplexId &tetId, const int &id, - SimplexId &triangleId) const = 0; + SimplexId &triangleId) const + = 0; int getTetrahedronTriangles( std::vector> &triangles) const; virtual int getTetrahedronNeighbor(const SimplexId &tetId, const int &localNeighborId, - SimplexId &neighborId) const = 0; + SimplexId &neighborId) const + = 0; - virtual SimplexId - getTetrahedronNeighborNumber(const SimplexId &tetId) const = 0; + virtual SimplexId getTetrahedronNeighborNumber(const SimplexId &tetId) const + = 0; int getTetrahedronNeighbors(std::vector> &neighbors); virtual int getTetrahedronVertex(const SimplexId &tetId, const int &localVertexId, - SimplexId &vertexId) const = 0; + SimplexId &vertexId) const + = 0; SimplexId getTriangleEdgeNumberInternal( const SimplexId & /*triangleId*/) const override { @@ -155,10 +159,12 @@ namespace ttk { virtual int getTriangleNeighbor(const SimplexId &triangleId, const int &localNeighborId, - SimplexId &neighborId) const = 0; + SimplexId &neighborId) const + = 0; virtual SimplexId - getTriangleNeighborNumber(const SimplexId &triangleId) const = 0; + getTriangleNeighborNumber(const SimplexId &triangleId) const + = 0; int getTriangleNeighbors(std::vector> &neighbors); diff --git a/core/base/integralLines/IntegralLines.h b/core/base/integralLines/IntegralLines.h index 077793af87..c8ef149c76 100644 --- a/core/base/integralLines/IntegralLines.h +++ b/core/base/integralLines/IntegralLines.h @@ -713,8 +713,8 @@ int ttk::IntegralLines::execute(triangulationType *triangulation) { int const taskNumber = (int)seedNumber_ / chunkSize_; #ifdef TTK_ENABLE_OPENMP4 #ifdef TTK_ENABLE_MPI -#pragma omp parallel shared( \ - ttk::intgl::finishedElement_, toSend_, ttk::intgl::addedElement_) \ +#pragma omp parallel shared( \ + ttk::intgl::finishedElement_, toSend_, ttk::intgl::addedElement_) \ num_threads(threadNumber_) { #else diff --git a/core/base/lowestCommonAncestor/LowestCommonAncestor.cpp b/core/base/lowestCommonAncestor/LowestCommonAncestor.cpp index 2fbf81c706..f6798aae04 100644 --- a/core/base/lowestCommonAncestor/LowestCommonAncestor.cpp +++ b/core/base/lowestCommonAncestor/LowestCommonAncestor.cpp @@ -57,7 +57,7 @@ int ttk::LowestCommonAncestor::RMQuery(const int &i, const int &j) const { // Position of the min in the blocs between the bloc of i and j min_pos[1] = ((blocJ - blocI) > 1) ? blocMinimumPosition_[blocMinimumValueRMQ_.query( - blocI + 1, blocJ - 1)] + blocI + 1, blocJ - 1)] : INT_MAX; // Position of the min in the bloc containing the jth case min_pos[2] diff --git a/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.cpp b/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.cpp index b9ce20b458..b6ebaef4ab 100644 --- a/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.cpp +++ b/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.cpp @@ -1016,11 +1016,15 @@ int MandatoryCriticalPoints::enumerateMandatorySaddles( #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { upperLca.preprocess(); } + { + upperLca.preprocess(); + } #ifdef TTK_ENABLE_OPENMP #pragma omp section #endif - { lowerLca.preprocess(); } + { + lowerLca.preprocess(); + } } // Link lists for each thread diff --git a/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.h b/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.h index 5ef92d0448..8d3b5c5cc6 100644 --- a/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.h +++ b/core/base/mandatoryCriticalPoints/MandatoryCriticalPoints.h @@ -1021,13 +1021,17 @@ int ttk::MandatoryCriticalPoints::buildSubTrees( #ifdef TTK_ENABLE_OPENMP #pragma omp critical #endif - { lowerMinimumList_.push_back(i); } + { + lowerMinimumList_.push_back(i); + } } if(isUpperMax) { #ifdef TTK_ENABLE_OPENMP #pragma omp critical #endif - { upperMaximumList_.push_back(i); } + { + upperMaximumList_.push_back(i); + } } } diff --git a/core/base/mergeTreeClustering/BranchMappingDistance.h b/core/base/mergeTreeClustering/BranchMappingDistance.h index de3e47c835..ab1e25d4bf 100644 --- a/core/base/mergeTreeClustering/BranchMappingDistance.h +++ b/core/base/mergeTreeClustering/BranchMappingDistance.h @@ -336,14 +336,15 @@ namespace ttk { if(tree1->getNumberOfChildren(curr1) == 0) { memT[curr1 + l * dim2 + nn2 * dim3 + 0 * dim4] = this->baseMetric_ == 0 ? editCost_Wasserstein1( - curr1, parent1, -1, -1, tree1, tree2) - : this->baseMetric_ == 1 ? editCost_Wasserstein2( - curr1, parent1, -1, -1, tree1, tree2) + curr1, parent1, -1, -1, tree1, tree2) + : this->baseMetric_ == 1 + ? editCost_Wasserstein2( + curr1, parent1, -1, -1, tree1, tree2) : this->baseMetric_ == 2 ? editCost_Persistence( - curr1, parent1, -1, -1, tree1, tree2) + curr1, parent1, -1, -1, tree1, tree2) : editCost_Shifting( - curr1, parent1, -1, -1, tree1, tree2); + curr1, parent1, -1, -1, tree1, tree2); } //----------------------------------------------------------------------- // If first subtree has more than one branch, try all decompositions @@ -378,14 +379,15 @@ namespace ttk { if(tree2->getNumberOfChildren(curr2) == 0) { memT[nn1 + 0 * dim2 + curr2 * dim3 + l * dim4] = this->baseMetric_ == 0 ? editCost_Wasserstein1( - -1, -1, curr2, parent2, tree1, tree2) - : this->baseMetric_ == 1 ? editCost_Wasserstein2( - -1, -1, curr2, parent2, tree1, tree2) + -1, -1, curr2, parent2, tree1, tree2) + : this->baseMetric_ == 1 + ? editCost_Wasserstein2( + -1, -1, curr2, parent2, tree1, tree2) : this->baseMetric_ == 2 ? editCost_Persistence( - -1, -1, curr2, parent2, tree1, tree2) + -1, -1, curr2, parent2, tree1, tree2) : editCost_Shifting( - -1, -1, curr2, parent2, tree1, tree2); + -1, -1, curr2, parent2, tree1, tree2); } //----------------------------------------------------------------------- // If first subtree has more than one branch, try all decompositions @@ -433,15 +435,17 @@ namespace ttk { if(tree1->getNumberOfChildren(curr1) == 0 and tree2->getNumberOfChildren(curr2) == 0) { memT[curr1 + l1 * dim2 + curr2 * dim3 + l2 * dim4] - = this->baseMetric_ == 0 ? editCost_Wasserstein1( - curr1, parent1, curr2, parent2, tree1, tree2) - : this->baseMetric_ == 1 ? editCost_Wasserstein2( - curr1, parent1, curr2, parent2, tree1, tree2) + = this->baseMetric_ == 0 + ? editCost_Wasserstein1( + curr1, parent1, curr2, parent2, tree1, tree2) + : this->baseMetric_ == 1 + ? editCost_Wasserstein2( + curr1, parent1, curr2, parent2, tree1, tree2) : this->baseMetric_ == 2 ? editCost_Persistence( - curr1, parent1, curr2, parent2, tree1, tree2) + curr1, parent1, curr2, parent2, tree1, tree2) : editCost_Shifting( - curr1, parent1, curr2, parent2, tree1, tree2); + curr1, parent1, curr2, parent2, tree1, tree2); } //--------------------------------------------------------------------------- // If first tree only has one branch, try all decompositions of @@ -655,12 +659,14 @@ namespace ttk { matchedNodes[m.first.first] = m.second.first; matchedNodes[m.first.second] = m.second.second; matchedCost[m.first.first] - = this->baseMetric_ == 0 ? editCost_Wasserstein1( - m.first.first, m.first.second, m.second.first, m.second.second, - tree1, tree2) - : this->baseMetric_ == 1 ? editCost_Wasserstein2( - m.first.first, m.first.second, m.second.first, - m.second.second, tree1, tree2) + = this->baseMetric_ == 0 + ? editCost_Wasserstein1(m.first.first, m.first.second, + m.second.first, + m.second.second, tree1, tree2) + : this->baseMetric_ == 1 + ? editCost_Wasserstein2(m.first.first, m.first.second, + m.second.first, + m.second.second, tree1, tree2) : this->baseMetric_ == 2 ? editCost_Persistence(m.first.first, m.first.second, m.second.first, diff --git a/core/base/mergeTreeClustering/MergeTreeClustering.h b/core/base/mergeTreeClustering/MergeTreeClustering.h index c6ec23903f..ccacb10eeb 100644 --- a/core/base/mergeTreeClustering/MergeTreeClustering.h +++ b/core/base/mergeTreeClustering/MergeTreeClustering.h @@ -288,7 +288,7 @@ namespace ttk { #ifdef TTK_ENABLE_OPENMP4 #pragma omp parallel for schedule(dynamic) \ shared(centroids, centroids2, oldCentroids_, oldCentroids2_) \ - num_threads(this->threadNumber_) if(parallelize_) + num_threads(this->threadNumber_) if(parallelize_) #endif for(unsigned int i = 0; i < centroids.size(); ++i) { std::vector> matching, @@ -341,7 +341,7 @@ namespace ttk { for(unsigned int i = 0; i < trees.size(); ++i) identified[i] = (upperBound_[i] <= centroidScore[bestCentroid_[i]]); - // Step 3 + // Step 3 #ifdef TTK_ENABLE_OPENMP4 #pragma omp parallel for schedule(dynamic) shared(centroids, centroids2) \ num_threads(this->threadNumber_) if(parallelize_) diff --git a/core/base/mergeTreeClustering/MergeTreeDistance.h b/core/base/mergeTreeClustering/MergeTreeDistance.h index 331d5d0b2a..8cdbb84b81 100644 --- a/core/base/mergeTreeClustering/MergeTreeDistance.h +++ b/core/base/mergeTreeClustering/MergeTreeDistance.h @@ -935,7 +935,7 @@ namespace ttk { #ifdef TTK_ENABLE_OPENMP4 #pragma omp task firstprivate(taskQueue, nodeT) UNTIED() \ shared(treeTable, forestTable, treeBackTable, forestBackTable, \ - treeChildDone, treeNodeDone) if(isTree1) + treeChildDone, treeNodeDone) if(isTree1) { #endif const ftm::FTMTree_MT *treeT = (isTree1) ? tree1 : tree2; @@ -1089,7 +1089,7 @@ namespace ttk { #ifdef TTK_ENABLE_OPENMP4 #pragma omp task firstprivate(nodeT) UNTIED() \ shared(treeTable, forestTable, treeBackTable, forestBackTable, \ - treeChildDone, treeNodeDone) + treeChildDone, treeNodeDone) { #endif while((int)nodeT != -1) { diff --git a/core/base/numericalIntegralLines/CMakeLists.txt b/core/base/numericalIntegralLines/CMakeLists.txt new file mode 100644 index 0000000000..90a520a0c6 --- /dev/null +++ b/core/base/numericalIntegralLines/CMakeLists.txt @@ -0,0 +1,9 @@ +ttk_add_base_library(numericalIntegralLines + SOURCES + NumericalIntegralLines.cpp + HEADERS + NumericalIntegralLines.h + DEPENDS + geometry + triangulation + ) diff --git a/core/base/numericalIntegralLines/NumericalIntegralLines.cpp b/core/base/numericalIntegralLines/NumericalIntegralLines.cpp new file mode 100644 index 0000000000..a1119d8682 --- /dev/null +++ b/core/base/numericalIntegralLines/NumericalIntegralLines.cpp @@ -0,0 +1,11 @@ +#include + +using namespace std; +using namespace ttk; +using namespace nil; + +NumericalIntegralLines::NumericalIntegralLines() { + this->setDebugMsgPrefix("NumericalIntegralLines"); +} + +NumericalIntegralLines::~NumericalIntegralLines() = default; diff --git a/core/base/numericalIntegralLines/NumericalIntegralLines.h b/core/base/numericalIntegralLines/NumericalIntegralLines.h new file mode 100644 index 0000000000..8dd25d4f52 --- /dev/null +++ b/core/base/numericalIntegralLines/NumericalIntegralLines.h @@ -0,0 +1,1127 @@ +/// \ingroup base +/// \class ttk::NumericalIntegralLines +/// \author Julien Tierny +/// \date May 2026 +/// \date NumericalIntegralLines extractor wrapping the DiscreteGradient class. +/// +/// \brief TTK convenience class wrapping the DiscreteGradient class for +/// the easy extraction of vpaths. +/// +/// Given a simplexId and dimension, this class returns a descending (or +/// ascending) vpath started in the given input simplex. +/// +/// \sa NumericalIntegralLines.cpp %for an alternative integral line backend. +/// \sa DiscreteGradient.cpp %for the core mechanisms. +/// \sa ttkNumericalIntegralLines.cpp %for a usage example. +/// + +#pragma once + +// base code includes +#include +#include +// std includes +#include +#include +#include +#include + +namespace ttk { + namespace nil { + + struct PathPoint { + SimplexId simplexId_; + int simplexDimension_; + std::vector barycentricWeights_; + }; + + /// Status of an elementary advection step (see doGradientStep()). + enum StepStatus { + /// The advection carries on (possibly in another simplex). + REGULAR_STEP = 0, + /// The advection reached the boundary of the domain. + BOUNDARY_REACHED = 1, + /// The advection reached a maximum (a minimum if backward). + EXTREMUM_REACHED = 2 + }; + + class NumericalIntegralLines : virtual public Debug { + + public: + NumericalIntegralLines(); + ~NumericalIntegralLines() override; + + /** + * @brief Compute a single numerical integral line. + * + * @param seed (SimplexId, dimension) + * @param barycentricWeights Weights for the input seed. + * @param output Output integral line (vector of 3D points). + * @param isForawrd Forward or backward line (default: forward). + */ + template + int computeIntegralLine(const triangulationType *triangulation, + const std::pair &seed, + const std::vector &barycentricWeights, + std::vector &output, + const bool &isForward = false) const; + + /** + * @brief Compute the gradient of the piecewise linear scalar field, + * restricted to the affine hull of the input simplex. + * + * @param simplexDimension Dimension of the input simplex. + * @param simplexId Identifier of the input simplex. + * @param gradient Output 3D gradient vector. + * @param barycentricGradient Optional output expression of the gradient + * in the edge basis (v1 - v0, ... vd - v0) of the simplex. This is also + * the variation of the barycentric weights (but for the first one) + * induced by a displacement along the gradient. + */ + template + int computeNumericalGradient(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector &gradient, + std::vector *barycentricGradient + = nullptr) const; + + /** + * @brief Elementary step of advection. + * + * Since the gradient of a piecewise linear scalar field is constant + * within a simplex, the integral line is a straight segment there. + * Hence, this step is integrated exactly: the current point is advected + * within its simplex until it reaches its boundary. Then, the simplex in + * which the advection carries on is identified (along with the + * barycentric coordinates of the advected point within it). + * + * @param current Input point (simplex, dimension, barycentric weights). + * @param isForward Forward or backward advection. + * @param next Output point (simplex, dimension, barycentric weights). + * @return StepStatus upon success (negative values otherwise). When the + * advection cannot carry on (boundary of the domain or extremum), the + * output simplex is the input one (with updated barycentric weights). + */ + template + int doGradientStep(const triangulationType *triangulation, + const PathPoint ¤t, + const bool &isForward, + PathPoint &next) const; + + /** + * @brief Compute numerical integral lines. + * + * @param output Vector storing the output vpaths (1 entry per seed, + * with possibly multiple v-path per seed). + * @param isForward Forward or backward vpath (default: forward). + */ + template + int execute(const triangulationType *triangulation, + const std::vector> &seeds, + std::vector> &output, + const bool &isForward = false) const; + + /** + * @brief Compute the variation of the barycentric weights of a point + * advected within the input simplex, per unit of arc length. + * + * @param slope Optional output slope of the scalar field along the + * (unit) advection direction. + * @return 0 if the advection can be carried on within the simplex + * (negative values otherwise, in particular if the restriction of the + * scalar field to the simplex is uniform). + */ + template + int getBarycentricVelocity(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const bool &isForward, + std::vector &velocity, + float *slope = nullptr) const; + + /** + * @brief Identify the simplex which takes the flow over at the input + * point. + * + * The flow is taken over by the coface of the input simplex which + * maximizes the slope of the scalar field (the steepest one) among the + * cofaces which admit the advection. Note that the gradient restricted + * to a face of a simplex is the projection of the gradient of the + * simplex: the slope of a coface is therefore always larger than (or + * equal to) that of its faces. + * + * @param point Input point (simplex, dimension, barycentric weights). + * @param isForward Forward or backward advection. + * @param next Output point (same location, expressed in the coface). + * @return 0 if a coface takes the flow over (negative values otherwise, + * in which case \p next is left untouched). + */ + template + int getFlowSimplex(const triangulationType *triangulation, + const PathPoint &point, + const bool &isForward, + PathPoint &next) const; + + /** + * @brief Retrieve the identifier of the face of the input simplex which + * is spanned by the input vertices. + */ + template + int getFaceIdentifier(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const std::vector &faceVertices, + const int &faceDimension, + SimplexId &faceIdentifier) const; + + /** + * @brief Compute the 3D coordinates of a path point. + */ + template + int getPointCoordinates(const triangulationType *triangulation, + const PathPoint &point, + std::array &coordinates) const; + + /** + * @brief Retrieve the cofaces of the input simplex (i.e. the simplices + * of the star of the input simplex, of higher dimension), as a list of + * (identifier, dimension) pairs. + */ + template + int getCofaces(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector> &cofaces) const; + + /** + * @brief Retrieve the face of the input simplex which supports the point + * of input barycentric weights (i.e. the face spanned by the vertices of + * non-zero weight). + */ + template + int getSubSimplex(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const std::vector &barycentricWeights, + PathPoint &subSimplex) const; + + template + int getVertexIdentifiers(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector &vertexIdentifiers) const; + + /** + * @brief Check if the input simplex is on the boundary of the domain. + */ + template + bool isOnDomainBoundary(const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId) const; + + /** + * @brief Check if an advection of input velocity can be carried on + * within a simplex, from a point of input barycentric weights (i.e. the + * advection does not immediately leave the simplex). + */ + static inline bool + isMotionAdmissible(const std::vector &barycentricWeights, + const std::vector &velocity) { + + float maximumVelocity = 0; + for(int i = 0; i < (int)velocity.size(); i++) + if(std::abs(velocity[i]) > maximumVelocity) + maximumVelocity = std::abs(velocity[i]); + + if(!(maximumVelocity > 0)) + return false; + + for(int i = 0; i < (int)barycentricWeights.size(); i++) + if((barycentricWeights[i] <= barycentricEpsilon_) + && (velocity[i] < -relativeEpsilon_ * maximumVelocity)) + // the advection immediately exits through the i-th face + return false; + + return true; + } + + /** + * @brief Express the barycentric weights of a point, given for a + * simplex, in the basis of one of its cofaces. + */ + static inline int + mapBarycentricWeights(const std::vector &sourceVertices, + const std::vector &sourceWeights, + const std::vector &targetVertices, + std::vector &targetWeights) { + + targetWeights.clear(); + targetWeights.resize(targetVertices.size(), 0); + + for(int i = 0; i < (int)sourceVertices.size(); i++) { + bool isFound = false; + for(int j = 0; j < (int)targetVertices.size(); j++) { + if(targetVertices[j] == sourceVertices[i]) { + targetWeights[j] = sourceWeights[i]; + isFound = true; + break; + } + } + if(!isFound) + // the source simplex is not a face of the target one + return -1; + } + + return 0; + } + + /** + * @brief Triangulation preconditioning. + */ + inline void + preconditionTriangulation(AbstractTriangulation *triangulation) { + + if(triangulation == nullptr) + return; + + // precondition simplex2face + triangulation->preconditionEdges(); + triangulation->preconditionCellEdges(); + + // precondition face2cofacets + triangulation->preconditionVertexEdges(); + triangulation->preconditionVertexStars(); + triangulation->preconditionEdgeStars(); + + // precondition boundary + triangulation->preconditionBoundaryVertices(); + if(triangulation->getDimensionality() > 1) + // in 1D, edges are cells (and the boundary is made of vertices) + triangulation->preconditionBoundaryEdges(); + + if(triangulation->getDimensionality() == 3) { + triangulation->preconditionTriangles(); + triangulation->preconditionTriangleEdges(); + triangulation->preconditionCellTriangles(); + triangulation->preconditionVertexTriangles(); + triangulation->preconditionEdgeTriangles(); + triangulation->preconditionTriangleStars(); + triangulation->preconditionBoundaryTriangles(); + } + } + + inline void setInputScalarField(const void *const scalars) { + scalars_ = scalars; + } + + protected: + /// Below this value, a barycentric weight is considered as null. + static constexpr float barycentricEpsilon_{1e-6}; + /// Relative tolerance used for the null tests on the velocity. + static constexpr float relativeEpsilon_{1e-6}; + + int maximumIterationNumber_{1000000000}; + /// Number of consecutive steps without any motion after which the + /// advection is considered as arbitrarily close to an extremum. + int maximumStalledStepNumber_{8}; + const void *scalars_{}; + }; + } // namespace nil +} // namespace ttk + +template +int ttk::nil::NumericalIntegralLines::computeIntegralLine( + const triangulationType *triangulation, + const std::pair &seed, + const std::vector &startBarycentricWeights, + std::vector &output, + const bool &isForward) const { + + output.clear(); + +#ifndef TTK_ENABLE_KAMIKAZE + if(triangulation == nullptr) + return -1; + if(scalars_ == nullptr) + return -2; + if((seed.second < 0) || (seed.second > triangulation->getDimensionality())) + return -3; +#endif + + PathPoint current; + current.simplexId_ = seed.first; + current.simplexDimension_ = seed.second; + current.barycentricWeights_ = startBarycentricWeights; + + if((int)current.barycentricWeights_.size() != seed.second + 1) + // no valid input coordinates: start from the barycenter of the seed + current.barycentricWeights_.assign( + seed.second + 1, 1.0 / (seed.second + 1)); + ttk::Geometry::normalizeBarycentricWeights(current.barycentricWeights_); + + // the seed simplex is where the integral line starts, not necessarily the + // simplex within which the advection takes place. generically, a point in + // the interior of a face is immediately advected within one of its cofaces + // (for instance, the mid-point of an edge is taken over by one of the + // triangles of its star). hand the flow over right away: this does not move + // the point, it only re-expresses it in the coface. if no coface admits the + // advection, the seed simplex constrains the flow and is kept as is. + PathPoint flowPoint; + if(getFlowSimplex( + triangulation, current, isForward, flowPoint) + == 0) + current = flowPoint; + + output.push_back(current); + + std::array previousCoordinates{}, currentCoordinates{}; + getPointCoordinates(triangulation, current, previousCoordinates); + + // an advection step may legitimately not move the current point (it can + // simply update the simplex supporting it, for instance when leaving a + // vertex for one of the tetrahedra of its star). however, a point which no + // longer moves is arbitrarily close to an extremum. + int stalledStepNumber = 0; + float pathLength = 0; + + int step = 0, status = REGULAR_STEP; + + for(step = 0; step < maximumIterationNumber_; step++) { + + PathPoint next; + + status = doGradientStep( + triangulation, current, isForward, next); + + if(status < 0) + return status; + + getPointCoordinates(triangulation, next, currentCoordinates); + + const float stepLength = Geometry::distance( + previousCoordinates.data(), currentCoordinates.data()); + + if(stepLength > pathLength * std::numeric_limits::epsilon()) { + output.push_back(next); + pathLength += stepLength; + stalledStepNumber = 0; + } else { + // the point did not move: only update the simplex supporting it + output.back() = next; + stalledStepNumber++; + } + + current = next; + previousCoordinates = currentCoordinates; + + if(status != REGULAR_STEP) + // the advection either left the domain or reached an extremum + break; + + if(stalledStepNumber > maximumStalledStepNumber_) { + // the advection is arbitrarily close to an extremum + status = EXTREMUM_REACHED; + break; + } + } + + if(step == maximumIterationNumber_) { +#ifdef TTK_ENABLE_OPENMP +#pragma omp critical +#endif + printWrn("Maximum iteration number reached for seed-#" + + std::to_string(seed.first) + + " (dim: " + std::to_string(seed.second) + ")."); + } + + return 0; +} + +// TODO +// move that function to the geometry class + +template +int ttk::nil::NumericalIntegralLines::computeNumericalGradient( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector &gradient, + std::vector *barycentricGradient) const { + + gradient = {0, 0, 0}; + + if(barycentricGradient) + barycentricGradient->assign(simplexDimension, 0); + + if(!simplexDimension) + return -1; + + std::vector vertexIdentifiers; + + getVertexIdentifiers( + triangulation, simplexDimension, simplexId, vertexIdentifiers); + + const int vertexNumber = vertexIdentifiers.size(); + + std::vector> vertexPoints(vertexNumber); + std::vector vertexScalars(vertexNumber); + + for(int i = 0; i < (int)vertexNumber; i++) { + triangulation->getVertexPoint(vertexIdentifiers[i], vertexPoints[i][0], + vertexPoints[i][1], vertexPoints[i][2]); + vertexScalars[i] = ((const dataType *)scalars_)[vertexIdentifiers[i]]; + } + + // build edge vectors and corresponding differences, wrt v0 + std::vector> edgeVectors(simplexDimension); + std::vector edgeDifferences(simplexDimension); + + for(int i = 0; i < simplexDimension; i++) { + for(int c = 0; c < 3; c++) + edgeVectors[i][c] = vertexPoints[i + 1][c] - vertexPoints[0][c]; + edgeDifferences[i] + = ((float)vertexScalars[i + 1]) - ((float)vertexScalars[0]); + } + + // Gram matrix gramMatrix[i][j] = edgeVectors[i] . edgeVectors[j] + std::vector> gramMatrix( + simplexDimension, std::vector(simplexDimension, 0)); + + float maximumDiagonalEntry = 0; + + for(int i = 0; i < simplexDimension; i++) { + for(int j = 0; j < simplexDimension; j++) + gramMatrix[i][j] = ttk::Geometry::dotProduct( + edgeVectors[i].data(), edgeVectors[j].data()); + + if(gramMatrix[i][i] > maximumDiagonalEntry) + maximumDiagonalEntry = gramMatrix[i][i]; + } + + if(!(maximumDiagonalEntry > 0)) + // degenerated simplex + return -2; + + // Gaussian elimintation + std::vector> augmentedMatrix( + simplexDimension, std::vector(simplexDimension + 1)); + for(int i = 0; i < simplexDimension; ++i) { + for(int j = 0; j < simplexDimension; ++j) + augmentedMatrix[i][j] = gramMatrix[i][j]; + augmentedMatrix[i][simplexDimension] = edgeDifferences[i]; + } + + for(int col = 0; col < simplexDimension; col++) { + // Partial pivot + int pivot = col; + for(int row = col + 1; row < simplexDimension; row++) + if(std::abs(augmentedMatrix[row][col]) + > std::abs(augmentedMatrix[pivot][col])) + pivot = row; + std::swap(augmentedMatrix[col], augmentedMatrix[pivot]); + + const float diagVal = augmentedMatrix[col][col]; + if(std::abs(diagVal) < powf(10, -FLT_DIG) * maximumDiagonalEntry) + // degenerated simplex + return -2; + + for(int row = 0; row < simplexDimension; row++) { + if(row == col) + continue; + const float factor = augmentedMatrix[row][col] / diagVal; + for(int j = col; j <= simplexDimension; ++j) + augmentedMatrix[row][j] -= factor * augmentedMatrix[col][j]; + } + } + + std::vector alpha(simplexDimension); + for(int i = 0; i < simplexDimension; i++) + alpha[i] = augmentedMatrix[i][simplexDimension] / augmentedMatrix[i][i]; + + // reconstruct the 3D gradient + for(int i = 0; i < simplexDimension; i++) + for(int c = 0; c < 3; c++) + gradient[c] += alpha[i] * edgeVectors[i][c]; + + if(barycentricGradient) + *barycentricGradient = alpha; + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::doGradientStep( + const triangulationType *triangulation, + const PathPoint ¤t, + const bool &isForward, + PathPoint &next) const { + + next = current; + + // 1) advection within the current simplex. + // the gradient of a piecewise linear scalar field is constant within a + // simplex. hence, the integral line is a straight segment there, which can + // be integrated exactly: the point is advected until it reaches the boundary + // of the simplex. + std::vector weights = current.barycentricWeights_; + std::vector velocity; + bool hasMoved = false; + + if(getBarycentricVelocity( + triangulation, current.simplexDimension_, current.simplexId_, isForward, + velocity) + == 0) { + + float maximumVelocity = 0; + for(int i = 0; i <= current.simplexDimension_; i++) + if(std::abs(velocity[i]) > maximumVelocity) + maximumVelocity = std::abs(velocity[i]); + + // largest arc length which maintains the point within the simplex + float travelDistance = std::numeric_limits::infinity(); + for(int i = 0; i <= current.simplexDimension_; i++) { + if(velocity[i] < -relativeEpsilon_ * maximumVelocity) { + const float distance = current.barycentricWeights_[i] / (-velocity[i]); + if(distance < travelDistance) + travelDistance = distance; + } + } + + if((travelDistance > 0) + && (travelDistance < std::numeric_limits::infinity())) { + + for(int i = 0; i <= current.simplexDimension_; i++) + weights[i] + = current.barycentricWeights_[i] + travelDistance * velocity[i]; + ttk::Geometry::normalizeBarycentricWeights(weights); + + hasMoved = true; + } + } + + // the advected point is now supported by a face of the current simplex + PathPoint exitPoint; + if(getSubSimplex(triangulation, current.simplexDimension_, current.simplexId_, + weights, exitPoint) + < 0) + return -1; + + // 2) identify the simplex in which the advection carries on. + if(getFlowSimplex( + triangulation, exitPoint, isForward, next) + == 0) + return REGULAR_STEP; + + // no coface takes the flow over. + const bool isOnBoundary = isOnDomainBoundary( + triangulation, exitPoint.simplexDimension_, exitPoint.simplexId_); + + if((!isOnBoundary) && (exitPoint.simplexDimension_ > 0) + && ((exitPoint.simplexDimension_ != current.simplexDimension_) + || (exitPoint.simplexId_ != current.simplexId_))) { + + // in the interior of the domain, the flow is then constrained to the exit + // face itself (typically, two cells whose gradients both point towards + // their common face). + if(getBarycentricVelocity( + triangulation, exitPoint.simplexDimension_, exitPoint.simplexId_, + isForward, velocity) + == 0) { + + if(isMotionAdmissible(exitPoint.barycentricWeights_, velocity)) { + next = exitPoint; + return REGULAR_STEP; + } + } + } + + // the advection stops here: report the advected point within the current + // simplex (same simplex, different barycentric weights). + next.simplexId_ = current.simplexId_; + next.simplexDimension_ = current.simplexDimension_; + next.barycentricWeights_ = weights; + + if(hasMoved && isOnBoundary) + return BOUNDARY_REACHED; + + return EXTREMUM_REACHED; +} + +template +int ttk::nil::NumericalIntegralLines::execute( + const triangulationType *triangulation, + const std::vector> &seeds, + std::vector> &output, + const bool &isForward) const { + + Timer t; + + output.resize(seeds.size()); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) +#endif + for(int i = 0; i < (int)seeds.size(); i++) { + std::vector barycentricWeights( + seeds[i].second + 1, 1.0 / (seeds[i].second + 1)); + computeIntegralLine( + triangulation, seeds[i], barycentricWeights, output[i], isForward); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp critical +#endif + printMsg(" - Seed-#" + std::to_string(seeds[i].first) + + " (dim: " + std::to_string(seeds[i].second) + + ", f: " + std::to_string(isForward) + + "): " + std::to_string(output[i].size()) + " point(s).", + debug::Priority::DETAIL); + } + + printMsg("Computed from " + std::to_string(output.size()) + " seed(s)", 1, + t.getElapsedTime(), threadNumber_); + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getBarycentricVelocity( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const bool &isForward, + std::vector &velocity, + float *slope) const { + + velocity.clear(); + + if(slope) + *slope = 0; + + if(simplexDimension < 1) + // no motion within a vertex + return -1; + + std::vector gradient, barycentricGradient; + + if(computeNumericalGradient( + triangulation, simplexDimension, simplexId, gradient, + &barycentricGradient) + < 0) + return -2; + + // along the unit advection direction, the variation of the scalar field is + // given by the magnitude of the gradient + const float magnitude = ttk::Geometry::magnitude(gradient.data()); + + if(!(magnitude > 0)) + // uniform scalar field: no motion + return -3; + + if(slope) + *slope = magnitude; + + // unit speed advection (the integration variable is the arc length) + const float scale = (isForward ? 1.0 : -1.0) / magnitude; + + velocity.resize(simplexDimension + 1, 0); + for(int i = 0; i < simplexDimension; i++) { + velocity[i + 1] = scale * barycentricGradient[i]; + // the barycentric weights sum up to 1 + velocity[0] -= scale * barycentricGradient[i]; + } + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getFlowSimplex( + const triangulationType *triangulation, + const PathPoint &point, + const bool &isForward, + PathPoint &next) const { + + // the flow is taken over by the coface of the input simplex which maximizes + // the slope of the scalar field (the steepest one) among the cofaces which + // admit the advection. note that the gradient restricted to a face of a + // simplex is the projection of the gradient of the simplex: the slope of a + // coface is therefore always larger than (or equal to) that of its faces. + std::vector> cofaces; + getCofaces(triangulation, point.simplexDimension_, point.simplexId_, cofaces); + + std::vector pointVertices, cofaceVertices; + getVertexIdentifiers( + triangulation, point.simplexDimension_, point.simplexId_, pointVertices); + + std::vector velocity; + + PathPoint bestPoint; + bestPoint.simplexDimension_ = -1; + float bestSlope = 0; + + for(int i = 0; i < (int)cofaces.size(); i++) { + + PathPoint candidate; + candidate.simplexId_ = cofaces[i].first; + candidate.simplexDimension_ = cofaces[i].second; + + getVertexIdentifiers(triangulation, candidate.simplexDimension_, + candidate.simplexId_, cofaceVertices); + + // express the advected point in the barycentric basis of the coface + if(mapBarycentricWeights(pointVertices, point.barycentricWeights_, + cofaceVertices, candidate.barycentricWeights_) + < 0) + continue; + + float slope = 0; + if(getBarycentricVelocity( + triangulation, candidate.simplexDimension_, candidate.simplexId_, + isForward, velocity, &slope) + < 0) + continue; + + if(!isMotionAdmissible(candidate.barycentricWeights_, velocity)) + // the advection would immediately leave this coface + continue; + + // steepest slope: along a unit direction, the variation of the scalar + // field is given by the magnitude of the gradient. + // ties (the gradient of the coface is aligned with one of its faces) are + // settled in favor of the coface of highest dimension (i.e. the least + // constrained advection). + if((slope > bestSlope) + || ((slope > bestSlope * (1 - relativeEpsilon_)) + && (candidate.simplexDimension_ > bestPoint.simplexDimension_))) { + bestSlope = slope; + bestPoint = candidate; + } + } + + if(bestPoint.simplexDimension_ < 0) + // no coface takes the flow over + return -1; + + next = bestPoint; + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getFaceIdentifier( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const std::vector &faceVertices, + const int &faceDimension, + SimplexId &faceIdentifier) const { + + faceIdentifier = -1; + + if(faceDimension == simplexDimension) { + faceIdentifier = simplexId; + return 0; + } + + if(!faceDimension) { + faceIdentifier = faceVertices[0]; + return 0; + } + + const int cellDimension = triangulation->getDimensionality(); + + int faceNumber = 0; + if(faceDimension == 1) { + if(simplexDimension == cellDimension) + faceNumber = triangulation->getCellEdgeNumber(simplexId); + else + // edges of a triangle + faceNumber = 3; + } else if(faceDimension == 2) + // triangles of a tetrahedron + faceNumber = triangulation->getCellTriangleNumber(simplexId); + else + return -1; + + // identify the face of the simplex which spans the input vertices + std::vector candidateVertices; + + for(int i = 0; i < faceNumber; i++) { + + SimplexId candidateId = -1; + + if(faceDimension == 1) { + if(simplexDimension == cellDimension) + triangulation->getCellEdge(simplexId, i, candidateId); + else + triangulation->getTriangleEdge(simplexId, i, candidateId); + } else + triangulation->getCellTriangle(simplexId, i, candidateId); + + getVertexIdentifiers( + triangulation, faceDimension, candidateId, candidateVertices); + + bool isMatching = true; + for(int j = 0; j < (int)faceVertices.size(); j++) { + bool isFound = false; + for(int k = 0; k < (int)candidateVertices.size(); k++) { + if(candidateVertices[k] == faceVertices[j]) { + isFound = true; + break; + } + } + if(!isFound) { + isMatching = false; + break; + } + } + + if(isMatching) { + faceIdentifier = candidateId; + return 0; + } + } + + return -2; +} + +template +int ttk::nil::NumericalIntegralLines::getPointCoordinates( + const triangulationType *triangulation, + const PathPoint &point, + std::array &coordinates) const { + + coordinates = {0, 0, 0}; + + std::vector vertexIdentifiers; + + if(getVertexIdentifiers(triangulation, point.simplexDimension_, + point.simplexId_, vertexIdentifiers) + < 0) + return -1; + + for(int i = 0; i < (int)vertexIdentifiers.size(); i++) { + std::array vertexPoint; + triangulation->getVertexPoint( + vertexIdentifiers[i], vertexPoint[0], vertexPoint[1], vertexPoint[2]); + for(int c = 0; c < 3; c++) + coordinates[c] += point.barycentricWeights_[i] * vertexPoint[c]; + } + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getCofaces( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector> &cofaces) const { + + cofaces.clear(); + + const int cellDimension = triangulation->getDimensionality(); + + if((simplexDimension < 0) || (simplexDimension >= cellDimension)) + // a top dimensional cell has no coface + return -1; + + SimplexId cofaceId = -1; + + if(!simplexDimension) { + if(cellDimension > 1) { + // the edges of the star of the vertex + // (in 1D, edges are cells: they are collected below) + const SimplexId edgeNumber + = triangulation->getVertexEdgeNumber(simplexId); + for(SimplexId i = 0; i < edgeNumber; i++) { + triangulation->getVertexEdge(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, 1); + } + } + + if(cellDimension == 3) { + // the triangles of the star of the vertex + const SimplexId triangleNumber + = triangulation->getVertexTriangleNumber(simplexId); + for(SimplexId i = 0; i < triangleNumber; i++) { + triangulation->getVertexTriangle(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, 2); + } + } + + // the cells of the star of the vertex + const SimplexId starNumber = triangulation->getVertexStarNumber(simplexId); + for(SimplexId i = 0; i < starNumber; i++) { + triangulation->getVertexStar(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, cellDimension); + } + + return 0; + } + + if(simplexDimension == 1) { + if(cellDimension == 3) { + // the triangles of the star of the edge + const SimplexId triangleNumber + = triangulation->getEdgeTriangleNumber(simplexId); + for(SimplexId i = 0; i < triangleNumber; i++) { + triangulation->getEdgeTriangle(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, 2); + } + } + + // the cells of the star of the edge + const SimplexId starNumber = triangulation->getEdgeStarNumber(simplexId); + for(SimplexId i = 0; i < starNumber; i++) { + triangulation->getEdgeStar(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, cellDimension); + } + + return 0; + } + + // the cells of the star of the triangle + const SimplexId starNumber = triangulation->getTriangleStarNumber(simplexId); + for(SimplexId i = 0; i < starNumber; i++) { + triangulation->getTriangleStar(simplexId, i, cofaceId); + cofaces.emplace_back(cofaceId, cellDimension); + } + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getSubSimplex( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + const std::vector &barycentricWeights, + PathPoint &subSimplex) const { + + std::vector vertexIdentifiers; + + if(getVertexIdentifiers( + triangulation, simplexDimension, simplexId, vertexIdentifiers) + < 0) + return -1; + + // the point is supported by the face spanned by the vertices of non-zero + // barycentric weight + std::vector faceVertices; + std::vector faceWeights; + + for(int i = 0; i < (int)vertexIdentifiers.size(); i++) { + if(barycentricWeights[i] > barycentricEpsilon_) { + faceVertices.push_back(vertexIdentifiers[i]); + faceWeights.push_back(barycentricWeights[i]); + } + } + + if(faceVertices.empty()) { + // degenerated weights: fall back on the closest vertex + int closestVertex = 0; + for(int i = 1; i < (int)vertexIdentifiers.size(); i++) + if(barycentricWeights[i] > barycentricWeights[closestVertex]) + closestVertex = i; + faceVertices = {vertexIdentifiers[closestVertex]}; + faceWeights = {1}; + } + + const int faceDimension = faceVertices.size() - 1; + + if(faceDimension == simplexDimension) { + subSimplex.simplexId_ = simplexId; + subSimplex.simplexDimension_ = simplexDimension; + subSimplex.barycentricWeights_ = barycentricWeights; + ttk::Geometry::normalizeBarycentricWeights(subSimplex.barycentricWeights_); + return 0; + } + + SimplexId faceIdentifier = -1; + + if(getFaceIdentifier(triangulation, simplexDimension, simplexId, faceVertices, + faceDimension, faceIdentifier) + < 0) + return -2; + + subSimplex.simplexId_ = faceIdentifier; + subSimplex.simplexDimension_ = faceDimension; + + // the vertices of the face are not necessarily ordered as in the simplex + std::vector subVertexIdentifiers; + getVertexIdentifiers( + triangulation, faceDimension, faceIdentifier, subVertexIdentifiers); + + if(mapBarycentricWeights(faceVertices, faceWeights, subVertexIdentifiers, + subSimplex.barycentricWeights_) + < 0) + return -3; + + ttk::Geometry::normalizeBarycentricWeights(subSimplex.barycentricWeights_); + + return 0; +} + +template +int ttk::nil::NumericalIntegralLines::getVertexIdentifiers( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId, + std::vector &vertexIdentifiers) const { + + if((simplexDimension < 0) || (simplexDimension > 3)) + return -1; + + vertexIdentifiers.resize(simplexDimension + 1); + + if(simplexDimension == triangulation->getDimensionality()) { + // top dimensional simplex: use the (faster) cell accessors + // (in 2D, triangles are cells) + for(int i = 0; i < simplexDimension + 1; i++) + triangulation->getCellVertex(simplexId, i, vertexIdentifiers[i]); + return 0; + } + + switch(simplexDimension) { + case 0: + vertexIdentifiers[0] = simplexId; + break; + case 1: + triangulation->getEdgeVertex(simplexId, 0, vertexIdentifiers[0]); + triangulation->getEdgeVertex(simplexId, 1, vertexIdentifiers[1]); + break; + case 2: + triangulation->getTriangleVertex(simplexId, 0, vertexIdentifiers[0]); + triangulation->getTriangleVertex(simplexId, 1, vertexIdentifiers[1]); + triangulation->getTriangleVertex(simplexId, 2, vertexIdentifiers[2]); + break; + default: + return -1; + break; + } + + return 0; +} + +template +bool ttk::nil::NumericalIntegralLines::isOnDomainBoundary( + const triangulationType *triangulation, + const int &simplexDimension, + const SimplexId &simplexId) const { + + if(simplexDimension == triangulation->getDimensionality()) + // a top dimensional cell is never on the boundary of the domain + return false; + + switch(simplexDimension) { + case 0: + return triangulation->isVertexOnBoundary(simplexId); + case 1: + return triangulation->isEdgeOnBoundary(simplexId); + case 2: + return triangulation->isTriangleOnBoundary(simplexId); + } + + return false; +} diff --git a/core/base/periodicImplicitTriangulation/PeriodicImplicitTriangulation.h b/core/base/periodicImplicitTriangulation/PeriodicImplicitTriangulation.h index d0110a4536..c1f8b71dda 100644 --- a/core/base/periodicImplicitTriangulation/PeriodicImplicitTriangulation.h +++ b/core/base/periodicImplicitTriangulation/PeriodicImplicitTriangulation.h @@ -113,20 +113,23 @@ namespace ttk { virtual int getTetrahedronEdge(const SimplexId &tetId, const int &id, - SimplexId &edgeId) const = 0; + SimplexId &edgeId) const + = 0; int getTetrahedronEdges(std::vector> &edges) const; virtual int getTetrahedronTriangle(const SimplexId &tetId, const int &id, - SimplexId &triangleId) const = 0; + SimplexId &triangleId) const + = 0; int getTetrahedronTriangles( std::vector> &triangles) const; virtual int getTetrahedronNeighbor(const SimplexId &tetId, const int &localNeighborId, - SimplexId &neighborId) const = 0; + SimplexId &neighborId) const + = 0; SimplexId getTetrahedronNeighborNumber(const SimplexId &tetId) const; @@ -134,7 +137,8 @@ namespace ttk { virtual int getTetrahedronVertex(const SimplexId &tetId, const int &localVertexId, - SimplexId &vertexId) const = 0; + SimplexId &vertexId) const + = 0; SimplexId getTriangleEdgeNumberInternal( const SimplexId &ttkNotUsed(triangleId)) const override { @@ -157,7 +161,8 @@ namespace ttk { virtual int getTriangleNeighbor(const SimplexId &triangleId, const int &localNeighborId, - SimplexId &neighborId) const = 0; + SimplexId &neighborId) const + = 0; SimplexId getTriangleNeighborNumber(const SimplexId &triangleId) const; diff --git a/core/base/persistenceDiagram/PersistenceDiagram.h b/core/base/persistenceDiagram/PersistenceDiagram.h index 8508647398..65e2e80d7e 100644 --- a/core/base/persistenceDiagram/PersistenceDiagram.h +++ b/core/base/persistenceDiagram/PersistenceDiagram.h @@ -934,9 +934,8 @@ int ttk::PersistenceDiagram::executeDiscreteMorseSandwichMPI( && triangulation->getSimplexRank(lid, simplexType) == ttk::MPIrank_) { // Add the relevant data - struct dataResponse res { - .lid_ = element.lid_, .isBirth_ = element.isBirth_ - }; + struct dataResponse res{ + .lid_ = element.lid_, .isBirth_ = element.isBirth_}; ttk::SimplexId vLid = dmsMPI_.getCellGreaterVertex( Cell{element.dim_ + (1 - element.isBirth_), lid}, *triangulation); res.vertexGid_ = triangulation->getVertexGlobalId(vLid); diff --git a/core/base/planarGraphLayout/PlanarGraphLayout.h b/core/base/planarGraphLayout/PlanarGraphLayout.h index 55053adb33..c6b1ddd921 100644 --- a/core/base/planarGraphLayout/PlanarGraphLayout.h +++ b/core/base/planarGraphLayout/PlanarGraphLayout.h @@ -263,7 +263,9 @@ int ttk::PlanarGraphLayout::computeDotString( // --------------------------------------------------------------------------- // Build Dot String - { dotString = headString + nodeString + edgeString + rankString + "}"; } + { + dotString = headString + nodeString + edgeString + rankString + "}"; + } // Print Status this->printMsg("Generating DOT string", 1, t.getElapsedTime()); diff --git a/core/base/regularGridTriangulation/RegularGridTriangulation.cpp b/core/base/regularGridTriangulation/RegularGridTriangulation.cpp index e931676a65..26a20f2bb9 100644 --- a/core/base/regularGridTriangulation/RegularGridTriangulation.cpp +++ b/core/base/regularGridTriangulation/RegularGridTriangulation.cpp @@ -424,11 +424,9 @@ int ttk::RegularGridTriangulation::preconditionDistributedVertices() { localBBox_y_max{this->localGridOffset_[1]}, localBBox_z_max{this->localGridOffset_[2]}; #ifdef TTK_ENABLE_OPENMP -#pragma omp parallel for reduction( \ - min \ - : localBBox_x_min, localBBox_y_min, localBBox_z_min) \ - reduction(max \ - : localBBox_x_max, localBBox_y_max, localBBox_z_max) +#pragma omp parallel for reduction( \ + min : localBBox_x_min, localBBox_y_min, localBBox_z_min) \ + reduction(max : localBBox_x_max, localBBox_y_max, localBBox_z_max) #endif for(SimplexId lvid = 0; lvid < nLocVertices; ++lvid) { // only keep non-ghost vertices diff --git a/core/base/regularGridTriangulation/RegularGridTriangulation.h b/core/base/regularGridTriangulation/RegularGridTriangulation.h index 751e4a333a..d0a97c03f5 100644 --- a/core/base/regularGridTriangulation/RegularGridTriangulation.h +++ b/core/base/regularGridTriangulation/RegularGridTriangulation.h @@ -88,16 +88,20 @@ namespace ttk { float spacing_[3]; // virtual void vertexToPosition2d(const SimplexId vertex, - SimplexId p[2]) const = 0; - virtual void vertexToPosition(const SimplexId vertex, - SimplexId p[3]) const = 0; + SimplexId p[2]) const + = 0; + virtual void vertexToPosition(const SimplexId vertex, SimplexId p[3]) const + = 0; virtual void triangleToPosition2d(const SimplexId triangle, - SimplexId p[2]) const = 0; + SimplexId p[2]) const + = 0; virtual void triangleToPosition(const SimplexId triangle, const int k, - SimplexId p[3]) const = 0; + SimplexId p[3]) const + = 0; virtual void tetrahedronToPosition(const SimplexId tetrahedron, - SimplexId p[3]) const = 0; + SimplexId p[3]) const + = 0; SimplexId findEdgeFromVertices(const SimplexId v0, const SimplexId v1) const; diff --git a/core/base/ripsPersistenceDiagram/FastRipsPersistenceDiagram2.cpp b/core/base/ripsPersistenceDiagram/FastRipsPersistenceDiagram2.cpp index 96a4627a55..34c77397a2 100644 --- a/core/base/ripsPersistenceDiagram/FastRipsPersistenceDiagram2.cpp +++ b/core/base/ripsPersistenceDiagram/FastRipsPersistenceDiagram2.cpp @@ -152,8 +152,8 @@ void FastRipsPersistenceDiagram2::computeRips0And1Persistence( e.d)) { // RNG edge critical.push_back(e); rng_.push_back(e); - if constexpr(std::is_same_v || std::is_same_v) + if constexpr(std::is_same_v + || std::is_same_v) ph[1].emplace_back(e.e); } else { // not RNG edge : merge neighboring polygons const int poly1 = UF.find(e.f1); diff --git a/core/base/topologicalCompression/TopologicalCompression.cpp b/core/base/topologicalCompression/TopologicalCompression.cpp index 171d5bb2b1..f71d2455fc 100644 --- a/core/base/topologicalCompression/TopologicalCompression.cpp +++ b/core/base/topologicalCompression/TopologicalCompression.cpp @@ -639,9 +639,10 @@ int ttk::TopologicalCompression::WriteToFile(FILE *fp, numberOfVertices *= (1 + dataExtent[2 * i + 1] - dataExtent[2 * i]); NbVertices = numberOfVertices; - int const totalSize = usePersistence ? ComputeTotalSizeForPersistenceDiagram( - getMapping(), getCriticalConstraints(), zfpOnly, - getNbSegments(), getNbVertices(), zfpTolerance) + int const totalSize = usePersistence + ? ComputeTotalSizeForPersistenceDiagram( + getMapping(), getCriticalConstraints(), zfpOnly, + getNbSegments(), getNbVertices(), zfpTolerance) : useOther ? ComputeTotalSizeForOther() : 0; diff --git a/core/base/triangulation/Triangulation.cpp b/core/base/triangulation/Triangulation.cpp index 617b6c9b94..ebdfc9e612 100644 --- a/core/base/triangulation/Triangulation.cpp +++ b/core/base/triangulation/Triangulation.cpp @@ -43,9 +43,9 @@ Triangulation::Triangulation(const Triangulation &rhs) Triangulation::Triangulation(Triangulation &&rhs) noexcept : AbstractTriangulation( - std::move(*static_cast(&rhs))), - abstractTriangulation_{nullptr}, explicitTriangulation_{std::move( - rhs.explicitTriangulation_)}, + std::move(*static_cast(&rhs))), + abstractTriangulation_{nullptr}, + explicitTriangulation_{std::move(rhs.explicitTriangulation_)}, implicitTriangulation_{std::move(rhs.implicitTriangulation_)}, periodicImplicitTriangulation_{ std::move(rhs.periodicImplicitTriangulation_)}, diff --git a/core/base/vPaths/CMakeLists.txt b/core/base/vPaths/CMakeLists.txt new file mode 100644 index 0000000000..f63bf59778 --- /dev/null +++ b/core/base/vPaths/CMakeLists.txt @@ -0,0 +1,10 @@ +ttk_add_base_library(vPaths + SOURCES + VPaths.cpp + HEADERS + VPaths.h + DEPENDS + discreteGradient + geometry + triangulation + ) diff --git a/core/base/vPaths/VPaths.cpp b/core/base/vPaths/VPaths.cpp new file mode 100644 index 0000000000..fdde6cc6e1 --- /dev/null +++ b/core/base/vPaths/VPaths.cpp @@ -0,0 +1,11 @@ +#include + +using namespace std; +using namespace ttk; +using namespace vp; + +VPaths::VPaths() { + this->setDebugMsgPrefix("VPaths"); +} + +VPaths::~VPaths() = default; diff --git a/core/base/vPaths/VPaths.h b/core/base/vPaths/VPaths.h new file mode 100644 index 0000000000..d8c57b1950 --- /dev/null +++ b/core/base/vPaths/VPaths.h @@ -0,0 +1,122 @@ +/// \ingroup base +/// \class ttk::VPaths +/// \author Julien Tierny +/// \date May 2026 +/// \date VPaths extractor wrapping the DiscreteGradient class. +/// +/// \brief TTK convenience class wrapping the DiscreteGradient class for +/// the easy extraction of vpaths. +/// +/// Given a simplexId and dimension, this class returns a descending (or +/// ascending) vpath started in the given input simplex. +/// +/// \sa VPaths.cpp %for an alternative integral line backend. +/// \sa DiscreteGradient.cpp %for the core mechanisms. +/// \sa ttkVPaths.cpp %for a usage example. +/// + +#pragma once + +// base code includes +#include +#include +// std includes + +namespace ttk { + namespace vp { + + class VPaths : virtual public Debug { + + public: + VPaths(); + ~VPaths() override; + + // template + // int execute(triangulationType *triangulation); + + /** + * @brief Extract vpaths. + * + * @param output Vector storing the output vpaths (1 entry per seed, + * with possibly multiple v-path per seed). + * @param isForward Forward or backward vpath (default: forward). + */ + template + int execute(const triangulationType *triangulation, + const std::vector &seeds, + std::vector>> &output, + const bool &isForward = false); + + /** + * @brief Triangulation preconditioning. + */ + inline void + preconditionTriangulation(AbstractTriangulation *triangulation) { + + // see dms precondition + dcg_.preconditionTriangulation(triangulation); + } + + inline void setInputOffsets(const SimplexId *const offsets) { + this->dcg_.setInputOffsets(offsets); + } + + inline void setInputScalarField(const void *const scalars, + const size_t &mTime) { + this->dcg_.setInputScalarField(scalars, mTime); + } + + protected: + dcg::DiscreteGradient dcg_{}; + }; + } // namespace vp +} // namespace ttk + +template +int ttk::vp::VPaths::execute( + const triangulationType *triangulation, + const std::vector &seeds, + std::vector>> &output, + const bool &isForward) { + + // fetching discrete gradient (or pre-computing it) + dcg_.setDebugLevel(debugLevel_); + dcg_.setThreadNumber(threadNumber_); + dcg_.buildGradient(*triangulation, false, nullptr); + + Timer t; + + output.resize(seeds.size()); + + /* + * NOTE: + * when considering seeds of non-zero dimension, mutliple v-paths may exist + * for a given seed. + */ + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) schedule(dynamic) +#endif + for(int i = 0; i < (int)seeds.size(); i++) { + if(!isForward) { + dcg_.getAllDescendingPaths(seeds[i], output[i], *triangulation); + } else { + dcg_.getAllAscendingPaths(seeds[i], output[i], *triangulation); + } + +#ifdef TTK_ENABLE_OPENMP +#pragma omp critical +#endif + printMsg(" - Seed-#" + std::to_string(seeds[i].id_) + + " (dim: " + std::to_string(seeds[i].dim_) + + ", f: " + std::to_string(isForward) + + "): " + std::to_string(output[i].size()) + " path(s).", + debug::Priority::DETAIL); + } + + printMsg( + "Computed v-path(s) from " + std::to_string(output.size()) + " seed(s)", 1, + t.getElapsedTime(), threadNumber_); + + return 0; +} diff --git a/core/base/vectorSimplification/VectorSimplification.cpp b/core/base/vectorSimplification/VectorSimplification.cpp index 2e181eb45f..5bd12905f1 100644 --- a/core/base/vectorSimplification/VectorSimplification.cpp +++ b/core/base/vectorSimplification/VectorSimplification.cpp @@ -21,9 +21,10 @@ void ttk::VectorSimplification::displayStats( std::count_if(pairs.begin(), pairs.end(), [](const CandidatePair &a) { return a.type == 0; }))}, {" #Saddle-saddle pairs", - std::to_string(dim == 3 ? std::count_if( - pairs.begin(), pairs.end(), - [](const CandidatePair &a) { return a.type == 1; }) + std::to_string(dim == 3 ? std::count_if(pairs.begin(), pairs.end(), + [](const CandidatePair &a) { + return a.type == 1; + }) : 0)}, {" #Saddle-max pairs", std::to_string(std::count_if( diff --git a/core/base/vectorSimplification/VectorSimplification.h b/core/base/vectorSimplification/VectorSimplification.h index 8987941abc..af7b5a8e5f 100644 --- a/core/base/vectorSimplification/VectorSimplification.h +++ b/core/base/vectorSimplification/VectorSimplification.h @@ -654,31 +654,32 @@ void ttk::VectorSimplification::getAscSaddlePairs( const auto dim = this->dcvf_.getDimensionality(); auto saddle2ToMaxima - = dim == 3 - ? getSaddle2ToAscPair( - criticalSaddles, - [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) { - return triangulation.getTriangleStar(a, i, r); - }, - [&triangulation](const SimplexId a) { - return triangulation.getTriangleStarNumber(a); - }, - [&triangulation](const SimplexId a) { - return triangulation.isTriangleOnBoundary(a); - }, - triangulation, static_cast(0.0)) - : getSaddle2ToAscPair( - criticalSaddles, - [&triangulation](const SimplexId a, const SimplexId i, SimplexId &r) { - return triangulation.getEdgeStar(a, i, r); - }, - [&triangulation](const SimplexId a) { - return triangulation.getEdgeStarNumber(a); - }, - [&triangulation](const SimplexId a) { - return triangulation.isEdgeOnBoundary(a); - }, - triangulation, static_cast(0.0)); + = dim == 3 ? getSaddle2ToAscPair( + criticalSaddles, + [&triangulation]( + const SimplexId a, const SimplexId i, SimplexId &r) { + return triangulation.getTriangleStar(a, i, r); + }, + [&triangulation](const SimplexId a) { + return triangulation.getTriangleStarNumber(a); + }, + [&triangulation](const SimplexId a) { + return triangulation.isTriangleOnBoundary(a); + }, + triangulation, static_cast(0.0)) + : getSaddle2ToAscPair( + criticalSaddles, + [&triangulation]( + const SimplexId a, const SimplexId i, SimplexId &r) { + return triangulation.getEdgeStar(a, i, r); + }, + [&triangulation](const SimplexId a) { + return triangulation.getEdgeStarNumber(a); + }, + [&triangulation](const SimplexId a) { + return triangulation.isEdgeOnBoundary(a); + }, + triangulation, static_cast(0.0)); for(size_t i = 0; i < saddle2ToMaxima.size(); ++i) { auto &maxs = saddle2ToMaxima[i]; diff --git a/core/vtk/ttkContourAroundPoint/ttkContourAroundPoint.cpp b/core/vtk/ttkContourAroundPoint/ttkContourAroundPoint.cpp index 5db49f2145..6d36843366 100644 --- a/core/vtk/ttkContourAroundPoint/ttkContourAroundPoint.cpp +++ b/core/vtk/ttkContourAroundPoint/ttkContourAroundPoint.cpp @@ -122,7 +122,7 @@ bool ttkContourAroundPoint::preprocessPts(vtkUnstructuredGrid *nodes, if(!scalarBuf || !codeBuf) return false; - // ---- Cell data ---- // + // ---- Cell data ---- // #ifndef NDEBUG // each arc should of course be defined by exactly two vertices auto cells = arcs->GetCells(); diff --git a/core/vtk/ttkContourTreeAlignment/ttkContourTreeAlignment.cpp b/core/vtk/ttkContourTreeAlignment/ttkContourTreeAlignment.cpp index d76d664d63..d7453acee1 100644 --- a/core/vtk/ttkContourTreeAlignment/ttkContourTreeAlignment.cpp +++ b/core/vtk/ttkContourTreeAlignment/ttkContourTreeAlignment.cpp @@ -50,7 +50,9 @@ int ttkContourTreeAlignment::RequestData(vtkInformation *ttkNotUsed(request), //================================================================================================================== // Print status - { this->printMsg("RequestData"); } + { + this->printMsg("RequestData"); + } //================================================================================================================== // Prepare input diff --git a/core/vtk/ttkImportEmbeddingFromTable/ttkImportEmbeddingFromTable.cpp b/core/vtk/ttkImportEmbeddingFromTable/ttkImportEmbeddingFromTable.cpp index b54e4827d8..964df589a4 100644 --- a/core/vtk/ttkImportEmbeddingFromTable/ttkImportEmbeddingFromTable.cpp +++ b/core/vtk/ttkImportEmbeddingFromTable/ttkImportEmbeddingFromTable.cpp @@ -73,15 +73,15 @@ int ttkImportEmbeddingFromTable::RequestData( vtkDataArray *xarr = XColumn.empty() ? nullptr : vtkDataArray::SafeDownCast( - inputTable->GetColumnByName(XColumn.data())); + inputTable->GetColumnByName(XColumn.data())); vtkDataArray *yarr = YColumn.empty() ? nullptr : vtkDataArray::SafeDownCast( - inputTable->GetColumnByName(YColumn.data())); + inputTable->GetColumnByName(YColumn.data())); vtkDataArray *zarr = ZColumn.empty() ? nullptr : vtkDataArray::SafeDownCast( - inputTable->GetColumnByName(ZColumn.data())); + inputTable->GetColumnByName(ZColumn.data())); if(xarr == nullptr or yarr == nullptr or zarr == nullptr) { printErr("invalid input columns."); diff --git a/core/vtk/ttkIntegralLines/ttk.module b/core/vtk/ttkIntegralLines/ttk.module index daed81a33a..b50b6e3731 100644 --- a/core/vtk/ttkIntegralLines/ttk.module +++ b/core/vtk/ttkIntegralLines/ttk.module @@ -6,4 +6,6 @@ HEADERS ttkIntegralLines.h DEPENDS integralLines + numericalIntegralLines + vPaths ttkAlgorithm diff --git a/core/vtk/ttkIntegralLines/ttkIntegralLines.cpp b/core/vtk/ttkIntegralLines/ttkIntegralLines.cpp index a7432747aa..79177e41f5 100644 --- a/core/vtk/ttkIntegralLines/ttkIntegralLines.cpp +++ b/core/vtk/ttkIntegralLines/ttkIntegralLines.cpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include #include #include #include @@ -225,6 +227,315 @@ int ttkIntegralLines::RequestData(vtkInformation *ttkNotUsed(request), vtkDataArray *inputOffsets = this->GetOrderArray( domain, 0, triangulation, false, 1, ForceInputOffsetScalarField); + bool isRunningWithMPI = false; + +#ifdef TTK_ENABLE_MPI + isRunningWithMPI = ttk::isRunningWithMPI(); +#endif + + std::vector idSpareStorage{}; + ttk::SimplexId *identifiers = this->GetIdentifierArrayPtr( + ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds, + idSpareStorage); + + if(!isRunningWithMPI) { + + if(BackEnd == BACKEND::NUMERICAL) { + printMsg("Selected `numerical` backend."); + + ttk::nil::NumericalIntegralLines num; + + num.setDebugLevel(debugLevel_); + num.setThreadNumber(threadNumber_); + + // setup the mesh + num.preconditionTriangulation(triangulation); + + // setup the data + num.setInputScalarField(inputScalars->GetVoidPointer(0)); + + // setup the seeds (simplexId, dimension) + std::vector> seedCells( + seeds->GetNumberOfCells()); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel num_threads(threadNumber_) +#endif + { + // GetCell(vtkIdType) returns a shared internal cell object: only the + // vtkGenericCell overload is thread-safe. + vtkNew cell{}; + +#ifdef TTK_ENABLE_OPENMP +#pragma omp for +#endif + for(int i = 0; i < (int)seedCells.size(); i++) { + seeds->GetCell(i, cell); + seedCells[i].first = identifiers[i]; + seedCells[i].second = cell->GetCellDimension(); + } + } + + std::vector> outputPaths; + + int status{}; + ttkVtkTemplateMacro(inputScalars->GetDataType(), triangulation->getType(), + (status = num.execute( + static_cast(triangulation->getData()), + seedCells, outputPaths, + // isForward? + Direction == 0))); + + if(status) + return status; + + int pointNumber{0}; + for(auto &path : outputPaths) { + pointNumber += path.size(); + } + + vtkNew outputPathGeometry; + + vtkNew pointCoords{}; + vtkNew vertexSeedId{}; + vtkNew cellSeedId{}; + vtkNew vertexSimplexId{}; + vtkNew vertexSimplexDimension{}; + vtkNew cellSimplexNumber{}; + vtkNew vertexDistanceFromSeed{}; + vtkNew outputMaskField{}; + + pointCoords->SetNumberOfComponents(3); + pointCoords->SetNumberOfTuples(pointNumber); + + vertexSeedId->SetNumberOfComponents(1); + vertexSeedId->SetNumberOfTuples(pointNumber); + vertexSeedId->SetName("SeedIdentifier"); + + vertexSimplexId->SetNumberOfComponents(1); + vertexSimplexId->SetNumberOfTuples(pointNumber); + vertexSimplexId->SetName("SimplexIdentifier"); + + vertexSimplexDimension->SetNumberOfComponents(1); + vertexSimplexDimension->SetNumberOfTuples(pointNumber); + vertexSimplexDimension->SetName("SimplexDimension"); + + vertexDistanceFromSeed->SetNumberOfComponents(1); + vertexDistanceFromSeed->SetNumberOfTuples(pointNumber); + vertexDistanceFromSeed->SetName("DistanceFromSeed"); + + outputMaskField->SetNumberOfComponents(1); + outputMaskField->SetNumberOfTuples(pointNumber); + outputMaskField->SetName(ttk::MaskScalarFieldName); + + cellSeedId->SetName("SeedIdentifier"); + cellSimplexNumber->SetName("SimplexNumber"); + + int pointId = 0; + for(int i = 0; i < (int)outputPaths.size(); i++) { + + const int simplexNumber = outputPaths[i].size(); + double distanceFromSeed = 0; + std::array point{}, previousPoint{}; + + for(int j = 0; j < simplexNumber; j++) { + + const ttk::nil::PathPoint &pathPoint = outputPaths[i][j]; + + // the path points are expressed by their barycentric coordinates + // within their simplex + num.getPointCoordinates(triangulation, pathPoint, point); + + if(j) + distanceFromSeed + += ttk::Geometry::distance(previousPoint.data(), point.data()); + + pointCoords->SetTuple3(pointId, point[0], point[1], point[2]); + vertexSeedId->SetTuple1(pointId, (int)seedCells[i].first); + vertexSimplexId->SetTuple1(pointId, (int)pathPoint.simplexId_); + vertexSimplexDimension->SetTuple1( + pointId, pathPoint.simplexDimension_); + vertexDistanceFromSeed->SetTuple1(pointId, distanceFromSeed); + // mask out the extremities of the integral line + outputMaskField->SetTuple1( + pointId, ((!j) || (j == simplexNumber - 1)) ? 0 : 1); + pointId++; + + if(j) { + vtkIdType edgeIds[2] = {pointId - 2, pointId - 1}; + outputPathGeometry->InsertNextCell(VTK_LINE, 2, edgeIds); + cellSeedId->InsertNextValue((int)seedCells[i].first); + cellSimplexNumber->InsertNextValue(simplexNumber); + } + + previousPoint = point; + } + } + + vtkNew pointSet{}; + pointSet->SetData(pointCoords); + outputPathGeometry->SetPoints(pointSet); + outputPathGeometry->GetPointData()->AddArray(vertexSeedId); + outputPathGeometry->GetPointData()->AddArray(outputMaskField); + outputPathGeometry->GetPointData()->AddArray(vertexSimplexId); + outputPathGeometry->GetPointData()->AddArray(vertexSimplexDimension); + outputPathGeometry->GetPointData()->AddArray(vertexDistanceFromSeed); + outputPathGeometry->GetCellData()->AddArray(cellSeedId); + outputPathGeometry->GetCellData()->AddArray(cellSimplexNumber); + + output->ShallowCopy(outputPathGeometry); + + return 1; + } else if(BackEnd == BACKEND::DISCRETE) { + printMsg("Selected `discrete` backend."); + + ttk::vp::VPaths vpaths; + + vpaths.setDebugLevel(debugLevel_); + vpaths.setThreadNumber(threadNumber_); + + // setup the mesh + vpaths.preconditionTriangulation(triangulation); + + // setup the data + vpaths.setInputScalarField( + inputScalars->GetVoidPointer(0), inputScalars->GetMTime()); + vpaths.setInputOffsets( + static_cast(ttkUtils::GetVoidPointer(inputOffsets))); + + std::vector seedCells(seeds->GetNumberOfCells()); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel num_threads(threadNumber_) +#endif + { + // GetCell(vtkIdType) returns a shared internal cell object: only the + // vtkGenericCell overload is thread-safe. + vtkNew cell{}; + +#ifdef TTK_ENABLE_OPENMP +#pragma omp for +#endif + for(int i = 0; i < (int)seedCells.size(); i++) { + seeds->GetCell(i, cell); + seedCells[i].dim_ = cell->GetCellDimension(); + seedCells[i].id_ = identifiers[i]; + } + } + + std::vector>> outputPaths; + + int status{}; + ttkTemplateMacro( + triangulation->getType(), + status = vpaths.execute(static_cast(triangulation->getData()), + seedCells, outputPaths, + // isForward? + Direction == 0)); + + if(status) + return status; + + int pointNumber{0}; + for(auto &seedPaths : outputPaths) { + for(auto &path : seedPaths) { + pointNumber += path.size(); + } + } + + vtkNew outputPathGeometry; + + vtkNew pointCoords{}; + vtkNew vertexSeedId{}; + vtkNew cellSeedId{}; + vtkNew cellForkId{}; + vtkNew vertexSimplexId{}; + vtkNew vertexSimplexDimension{}; + vtkNew cellSimplexNumber{}; + vtkNew outputMaskField{}; + + pointCoords->SetNumberOfComponents(3); + pointCoords->SetNumberOfTuples(pointNumber); + + vertexSeedId->SetNumberOfComponents(1); + vertexSeedId->SetNumberOfTuples(pointNumber); + vertexSeedId->SetName("SeedIdentifier"); + + vertexSimplexId->SetNumberOfComponents(1); + vertexSimplexId->SetNumberOfTuples(pointNumber); + vertexSimplexId->SetName("SimplexIdentifier"); + + vertexSimplexDimension->SetNumberOfComponents(1); + vertexSimplexDimension->SetNumberOfTuples(pointNumber); + vertexSimplexDimension->SetName("SimplexDimension"); + + outputMaskField->SetNumberOfComponents(1); + outputMaskField->SetNumberOfTuples(pointNumber); + outputMaskField->SetName(ttk::MaskScalarFieldName); + + cellSeedId->SetName("SeedIdentifier"); + cellForkId->SetName("ForkIdentifier"); + cellSimplexNumber->SetName("SimplexNumber"); + + int pointId = 0; + int localSeedId = 0; + int pathPointId = 0; + for(auto &seedPaths : outputPaths) { + + int forkId = 0; + for(auto &path : seedPaths) { + + pathPointId = 0; + + int simplexNumber = path.size(); + + for(auto &c : path) { + float point[3]; + triangulation->getCellIncenter(c.id_, c.dim_, point); + pointCoords->SetTuple3(pointId, point[0], point[1], point[2]); + vertexSeedId->SetTuple1(pointId, (int)seedCells[localSeedId].id_); + vertexSimplexId->SetTuple1(pointId, (int)c.id_); + vertexSimplexDimension->SetTuple1(pointId, (int)c.dim_); + // mask out the extremities of the integral line + outputMaskField->SetTuple1( + pointId, + ((!pathPointId) || (pathPointId == simplexNumber - 1)) ? 0 : 1); + pointId++; + pathPointId++; + + if(pathPointId > 1) { + vtkIdType edgeIds[2] = {pointId - 2, pointId - 1}; + outputPathGeometry->InsertNextCell(VTK_LINE, 2, edgeIds); + cellSeedId->InsertNextValue((int)seedCells[localSeedId].id_); + cellForkId->InsertNextValue((int)forkId); + cellSimplexNumber->InsertNextValue((simplexNumber)); + } + } + forkId++; + } + localSeedId++; + } + + vtkNew pointSet{}; + pointSet->SetData(pointCoords); + outputPathGeometry->SetPoints(pointSet); + outputPathGeometry->GetPointData()->AddArray(vertexSeedId); + outputPathGeometry->GetPointData()->AddArray(outputMaskField); + outputPathGeometry->GetPointData()->AddArray(vertexSimplexId); + outputPathGeometry->GetPointData()->AddArray(vertexSimplexDimension); + outputPathGeometry->GetCellData()->AddArray(cellSeedId); + outputPathGeometry->GetCellData()->AddArray(cellForkId); + outputPathGeometry->GetCellData()->AddArray(cellSimplexNumber); + + output->ShallowCopy(outputPathGeometry); + + return 1; + } + } else { + if(BackEnd != BACKEND::ONESKELETON) + printWrn("Distributed run, defaulting to the `OneSkeleton` backend."); + } + const ttk::SimplexId numberOfPointsInDomain = domain->GetNumberOfPoints(); this->setVertexNumber(numberOfPointsInDomain); int numberOfPointsInSeeds = seeds->GetNumberOfPoints(); @@ -281,11 +592,11 @@ int ttkIntegralLines::RequestData(vtkInformation *ttkNotUsed(request), } numberOfPointsInSeeds = inputIdentifiers.size(); } else { - std::vector idSpareStorage{}; + std::vector lIdSpareStorage{}; ttk::SimplexId *inputIdentifierGlobalId; inputIdentifierGlobalId = this->GetIdentifierArrayPtr( ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds, - idSpareStorage); + lIdSpareStorage); ttk::SimplexId localId = 0; for(int i = 0; i < numberOfPointsInSeeds; i++) { localId = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]); @@ -302,21 +613,17 @@ int ttkIntegralLines::RequestData(vtkInformation *ttkNotUsed(request), this->setGlobalElementCounter(numberOfPointsInSeeds); inputIdentifiers.resize(numberOfPointsInSeeds); totalSeeds = numberOfPointsInSeeds; - std::vector idSpareStorage{}; + std::vector lIdSpareStorage{}; ttk::SimplexId *inputIdentifierGlobalId; inputIdentifierGlobalId = this->GetIdentifierArrayPtr( ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds, - idSpareStorage); + lIdSpareStorage); for(int i = 0; i < numberOfPointsInSeeds; i++) { inputIdentifiers.at(i) = triangulation->getVertexLocalId(inputIdentifierGlobalId[i]); } } #else - std::vector idSpareStorage{}; - ttk::SimplexId *identifiers = this->GetIdentifierArrayPtr( - ForceInputVertexScalarField, 2, ttk::VertexScalarFieldName, seeds, - idSpareStorage); std::unordered_set isSeed; for(ttk::SimplexId k = 0; k < numberOfPointsInSeeds; ++k) { isSeed.insert(identifiers[k]); diff --git a/core/vtk/ttkIntegralLines/ttkIntegralLines.h b/core/vtk/ttkIntegralLines/ttkIntegralLines.h index d3ffee41f9..60b710cd59 100644 --- a/core/vtk/ttkIntegralLines/ttkIntegralLines.h +++ b/core/vtk/ttkIntegralLines/ttkIntegralLines.h @@ -70,7 +70,10 @@ // ttk code includes #include +#include +#include #include +#include class vtkUnstructuredGrid; @@ -82,6 +85,12 @@ class TTKINTEGRALLINES_EXPORT ttkIntegralLines : public ttkAlgorithm, vtkTypeMacro(ttkIntegralLines, ttkAlgorithm); + enum class BACKEND { + ONESKELETON = 0, + NUMERICAL = 1, + DISCRETE = 2, + }; + vtkGetMacro(Direction, int); vtkSetMacro(Direction, int); @@ -91,6 +100,9 @@ class TTKINTEGRALLINES_EXPORT ttkIntegralLines : public ttkAlgorithm, vtkSetMacro(ForceInputOffsetScalarField, bool); vtkGetMacro(ForceInputOffsetScalarField, bool); + ttkSetEnumMacro(BackEnd, BACKEND); + vtkGetEnumMacro(BackEnd, BACKEND); + vtkSetMacro(EnableForking, bool); vtkGetMacro(EnableForking, bool); @@ -123,6 +135,7 @@ class TTKINTEGRALLINES_EXPORT ttkIntegralLines : public ttkAlgorithm, vtkInformationVector *outputVector) override; private: + BACKEND BackEnd{BACKEND::ONESKELETON}; int Direction{0}; bool ForceInputVertexScalarField{false}; bool ForceInputOffsetScalarField{false}; diff --git a/core/vtk/ttkMergeTreePrincipalGeodesicsDecoding/ttkMergeTreePrincipalGeodesicsDecoding.cpp b/core/vtk/ttkMergeTreePrincipalGeodesicsDecoding/ttkMergeTreePrincipalGeodesicsDecoding.cpp index 474cc6f4aa..a7d80a59d3 100644 --- a/core/vtk/ttkMergeTreePrincipalGeodesicsDecoding/ttkMergeTreePrincipalGeodesicsDecoding.cpp +++ b/core/vtk/ttkMergeTreePrincipalGeodesicsDecoding/ttkMergeTreePrincipalGeodesicsDecoding.cpp @@ -499,7 +499,7 @@ int ttkMergeTreePrincipalGeodesicsDecoding::runOutput( and (inputMTrees.empty() or baryMatchings_.empty())) printWrn("Please provide input trees and correlation matrix to transfer " "input trees information."); - // TODO fix if an interpolation is empty + // TODO fix if an interpolation is empty #ifdef TTK_ENABLE_OPENMP #pragma omp parallel for schedule(dynamic) num_threads(this->threadNumber_) #endif diff --git a/core/vtk/ttkPlanarGraphLayout/ttkMergeTreeVisualization.h b/core/vtk/ttkPlanarGraphLayout/ttkMergeTreeVisualization.h index 0c4b9adee0..83b47eeed7 100644 --- a/core/vtk/ttkPlanarGraphLayout/ttkMergeTreeVisualization.h +++ b/core/vtk/ttkPlanarGraphLayout/ttkMergeTreeVisualization.h @@ -1281,7 +1281,7 @@ class ttkMergeTreeVisualization : public ttk::MergeTreeVisualization { // Insert point // -------------- auto getPoint - = [&](vtkUnstructuredGrid *vtu, int pointID, double(&point)[3]) { + = [&](vtkUnstructuredGrid *vtu, int pointID, double (&point)[3]) { if(not vtu) return; if(not isPersistenceDiagram or convertedToDiagram) { diff --git a/core/vtk/ttkTriangulationRequest/ttkTriangulationRequest.cpp b/core/vtk/ttkTriangulationRequest/ttkTriangulationRequest.cpp index 8f141bd02e..18b49b66a5 100644 --- a/core/vtk/ttkTriangulationRequest/ttkTriangulationRequest.cpp +++ b/core/vtk/ttkTriangulationRequest/ttkTriangulationRequest.cpp @@ -260,7 +260,8 @@ int ttkTriangulationRequest::RequestData(vtkInformation *ttkNotUsed(request), case SIMPLEX::VERTEX: { const auto vid = addVertex(si); cells->InsertNextCell(VTK_VERTEX, 1, &vid); - cellIds->InsertNextTuple1(vid); + // report the identifier of the vertex (not that of the point) + cellIds->InsertNextTuple1(si); cellDims->InsertNextTuple1(0); } break; @@ -493,7 +494,8 @@ int ttkTriangulationRequest::RequestData(vtkInformation *ttkNotUsed(request), if(triangulation->isVertexOnBoundary(v)) { const auto vid = addVertex(v); cells->InsertNextCell(VTK_VERTEX, 1, &vid); - cellIds->InsertNextTuple1(vid); + // report the identifier of the vertex (not that of the point) + cellIds->InsertNextTuple1(v); cellDims->InsertNextTuple1(0); } } diff --git a/paraview/xmls/IntegralLines.xml b/paraview/xmls/IntegralLines.xml index 461e96dbfc..19b93cf4e5 100644 --- a/paraview/xmls/IntegralLines.xml +++ b/paraview/xmls/IntegralLines.xml @@ -13,13 +13,13 @@ lines of the gradient of an input scalar field."> The filter takes on its input a scalar field attached as point data to an input geometry (either 2D or 3D, either regular grids or triangulations) -and computes the forward or backward integral lines along the edges of the -input mesh, given a list of input sources. +and computes the forward or backward integral lines, given a list of input sources. Several backends are available: -The sources are specified with a vtkPointSet on which is attached as point -data a scalar field that represent the vertex identifiers of the sources in -the input geometry (use the Identifiers plugin prior to the selection of the -sources if necessary to create such an identifier field). +- the 'OneSkeleton' backend follows the flow along the edges of the mesh. This backend supports MPI computation. +- the 'Numerical' backend computes integral lines following the piecewise constant gradient of each simplex (i.e., it is going through the simplices) +- the 'Discrete' backend computes a discrete integral line (a.k.a. v-path), following the discrete Morse theory model. + +The sources are specified with a vtkPointSet, containing the simplices from which to start the integral lines. Each simplex is associated to its global identifier, passed as a point data array for vertices or cell data array in general for simplices of arbitrary dimension. Use PointAndCellIds to compute those in advance (the name of the generated field must then be provided to this filter). Online examples: @@ -110,7 +110,7 @@ sources if necessary to create such an identifier field). - Select the vertex identifier scalar field in the sources. + Select the seed identifier scalar field in the sources. @@ -178,6 +178,26 @@ plateaus). + + + + + + + + Backend for the computation of the integral lines. + - the 'OneSkeleton' backend follows the flow along the edges of the mesh. This backend supports MPI computation. +- the 'Numerical' backend computes integral lines following the piecewise constant gradient of each simplex (i.e., it is going through the simplices) +- the 'Discrete' backend computes a discrete integral line (a.k.a. v-path), following the discrete Morse theory model. + + + + + + Enables forking when the integral line comes accross a saddle vertex. In that case, it will spawn as many new integral lines as there are @@ -202,6 +225,7 @@ plateaus). +