diff --git a/.github/workflows/package.yml b/.github/workflows/package.yml index cc3d4baece..3c149e572e 100644 --- a/.github/workflows/package.yml +++ b/.github/workflows/package.yml @@ -180,7 +180,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-12, macos-14] + os: [macos-13, macos-14, macos-15] env: DYLD_LIBRARY_PATH: /usr/local/lib @@ -188,35 +188,36 @@ jobs: - uses: actions/checkout@v4 name: Checkout TTK source code - - name: Remove hosted Python - run: | - sudo rm -rf /Library/Frameworks/Python.framework/ - sudo rm -rf /usr/local/Frameworks/Python.framework/ + - uses: actions/setup-python@v5 + with: + python-version: '3.12' - name: Install macOS dependencies run: | # ParaView dependencies - brew reinstall python brew install --cask xquartz brew install llvm mesa glew qt@5 ninja # TTK dependencies brew install boost eigen graphviz spectra sqlite zlib numpy qhull + - name: Install and setup sccache + uses: f3d-app/sccache-setup-action@v1 + with: + key: macos-0 + - name: Install optional dependencies uses: ./.github/actions/install-deps-unix - name: Fetch & install TTK-ParaView run: | - if [[ "${{ matrix.os }}" == "macos-12" ]]; then - wget https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz - tar xzf ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz - fi - if [[ "${{ matrix.os }}" == "macos-14" ]]; then - wget https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}-arm64.tar.gz - tar xzf ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}-arm64.tar.gz - fi + wget https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz + tar xzf ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz sudo cp -r ttk-paraview/* /usr/local + # pvpython does not embed the correct PYTHONPATH + echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV pvpython -m pip install --break-system-packages scikit-learn + # pvpython is expecting a vtkpython executable at this path + sudo ln -s /usr/local/bin/pvpython /Library/Frameworks/Python.framework/Versions/3.12/vtkpython - name: Set compilers as environment variables run: | @@ -230,8 +231,6 @@ jobs: cmake \ -DCMAKE_BUILD_TYPE=Release \ -DQt5_DIR=$(brew --prefix qt@5)/lib/cmake/Qt5 \ - -DPython3_FIND_STRATEGY=LOCATION \ - -DPython3_ROOT_DIR=$(brew --prefix python) \ -DTTK_BUILD_PARAVIEW_PLUGINS=TRUE \ -DTTK_BUILD_VTK_WRAPPERS=TRUE \ -DTTK_BUILD_STANDALONE_APPS=TRUE \ @@ -266,36 +265,36 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-12, macos-14] + os: [macos-13, macos-14, macos-15] env: DYLD_LIBRARY_PATH: /usr/local/lib steps: - uses: actions/checkout@v4 name: Checkout TTK source code - - name: Remove hosted Python - run: | - sudo rm -rf /Library/Frameworks/Python.framework - sudo rm -rf /usr/local/Frameworks/Python.framework + - uses: actions/setup-python@v5 + with: + python-version: '3.12' - name: Install macOS dependencies run: | # ParaView dependencies - brew reinstall python brew install --cask xquartz brew install llvm mesa glew qt@5 ninja # TTK dependencies brew install boost eigen graphviz spectra sqlite zlib numpy qhull + - name: Install and setup sccache + uses: f3d-app/sccache-setup-action@v1 + with: + key: macos-0 + + - name: Install optional dependencies + uses: ./.github/actions/install-deps-unix + - name: Fetch TTK-ParaView run: | - if [[ "${{ matrix.os }}" == "macos-12" ]]; then - wget -O ttk-paraview.tar.gz https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz - fi - if [[ "${{ matrix.os }}" == "macos-14" ]]; then - wget -O ttk-paraview.tar.gz https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}-arm64.tar.gz - fi - + wget -O ttk-paraview.tar.gz https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-${{ env.PV_TAG }}-${{ matrix.os }}.tar.gz - name: Fetch TTK .tar.gz artifact uses: actions/download-artifact@v4.1.7 @@ -317,6 +316,8 @@ jobs: echo "CMAKE_PREFIX_PATH=$(brew --prefix qt@5)/lib/cmake:$CMAKE_PREFIX_PATH" >> $GITHUB_ENV # pvpython does not embed the correct PYTHONPATH echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV + # pvpython is expecting a vtkpython executable at this path + sudo ln -s /usr/local/bin/pvpython /Library/Frameworks/Python.framework/Versions/3.12/vtkpython - name: Run TTK tests uses: ./.github/actions/test-ttk-unix @@ -337,7 +338,7 @@ jobs: # remove examples which fill up the memory rm python/topologicalOptimization_darkSky.py # some cases fail only with version 12 - if [[ "${{ matrix.os }}" == "macos-12" ]]; then + if [[ "${{ matrix.os }}" == "macos-13" ]]; then rm python/contourTreeAlignment.py rm python/geometryApproximation.py rm python/harmonicSkeleton.py @@ -659,13 +660,13 @@ jobs: file: ttk-ubuntu-24.04.deb/ttk.deb asset_name: ttk-$tag-ubuntu-24.04.deb - - name: Upload MacOS 12 binary archives as Release Asset + - name: Upload MacOS 13 binary archives as Release Asset uses: svenstaro/upload-release-action@v2 with: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} - file: ttk-macos-12.tar.gz/ttk.tar.gz - asset_name: ttk-$tag-macos-12.tar.gz + file: ttk-macos-13.tar.gz/ttk.tar.gz + asset_name: ttk-$tag-macos-13.tar.gz - name: Upload MacOS 14 binary archives as Release Asset uses: svenstaro/upload-release-action@v2 @@ -673,7 +674,15 @@ jobs: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} file: ttk-macos-14.tar.gz/ttk.tar.gz - asset_name: ttk-$tag-macos-14-arm64.tar.gz + asset_name: ttk-$tag-macos-14.tar.gz + + - name: Upload MacOS 15 binary archives as Release Asset + uses: svenstaro/upload-release-action@v2 + with: + repo_token: ${{ secrets.GITHUB_TOKEN }} + tag: ${{ github.ref }} + file: ttk-macos-15.tar.gz/ttk.tar.gz + asset_name: ttk-$tag-macos-15.tar.gz - name: Upload Windows .exe as Release Asset uses: svenstaro/upload-release-action@v2 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index d0a9693994..aa2b8b3c15 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -271,7 +271,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-12, macos-14] + os: [macos-13, macos-14, macos-15] if: ${{ github.repository_owner == 'topology-tool-kit' || !contains(github.ref, 'heads') }} env: DYLD_LIBRARY_PATH: /usr/local/lib @@ -280,14 +280,13 @@ jobs: - uses: actions/checkout@v4 name: Checkout TTK source code - - name: Remove hosted Python - run: | - sudo rm -rf /Library/Frameworks/Python.framework/ - sudo rm -rf /usr/local/Frameworks/Python.framework/ + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + - name: Install macOS dependencies run: | # ParaView dependencies - brew reinstall python brew install --cask xquartz brew install llvm ninja open-mpi # TTK dependencies @@ -303,20 +302,18 @@ jobs: - name: Fetch TTK-ParaView headless macOS binary archive run: | - if [[ "${{ matrix.os }}" == "macos-12" ]]; then wget -O ttk-paraview-headless.tar.gz \ https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-headless-${{ matrix.os }}.tar.gz - fi - if [[ "${{ matrix.os }}" == "macos-14" ]]; then - wget -O ttk-paraview-headless.tar.gz \ - https://github.com/${{ env.PV_REPO }}/releases/download/${{ env.PV_TAG }}/ttk-paraview-headless-${{ matrix.os }}-arm64.tar.gz - fi - name: Install ParaView run: | tar xzf ttk-paraview-headless.tar.gz sudo cp -r ttk-paraview/* /usr/local + # pvpython does not embed the correct PYTHONPATH + echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV pvpython -m pip install --break-system-packages scikit-learn + # pvpython is expecting a vtkpython executable at this path + sudo ln -s /usr/local/bin/pvpython /Library/Frameworks/Python.framework/Versions/3.12/vtkpython - name: Set compilers as environment variables run: | @@ -355,11 +352,6 @@ jobs: shell: bash run: sccache --show-stats - - name: Set PYTHONPATH for macOS pvpython - run: | - # pvpython does not embed the correct PYTHONPATH - echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV - - name: Run TTK tests uses: ./.github/actions/test-ttk-unix @@ -371,7 +363,6 @@ jobs: name: Checkout ttk-data - name: Run ttk-data states [NOT ENFORCED] - if: false id: validate continue-on-error: true run: | @@ -380,6 +371,8 @@ jobs: rm ttk-data/states/nestedTrackingFromOverlap.py # remove examples which fill up the memory rm ttk-data/states/topologicalOptimization_darkSky.pvsm + # remove examples that take too long to render under mac + rm ttk-data/states/persistentGenerators_householdAnalysis.pvsm cd ttk-data/tests mkdir output_screenshots pvpython -u validate.py || (tar zcf screenshots.tar.gz output_screenshots && false) @@ -390,7 +383,7 @@ jobs: if: steps.validate.outcome == 'failure' uses: actions/upload-artifact@v4 with: - name: screenshots-macOS.tar.gz + name: screenshots-${{ matrix.os }}.tar.gz path: ttk-data/tests/screenshots.tar.gz retention-days: 10 diff --git a/core/base/assignmentSolver/AssignmentAuction.h b/core/base/assignmentSolver/AssignmentAuction.h index 341074fe5a..c68102bbe3 100644 --- a/core/base/assignmentSolver/AssignmentAuction.h +++ b/core/base/assignmentSolver/AssignmentAuction.h @@ -275,6 +275,7 @@ namespace ttk { initFirstRound(); while(not stoppingCriterion(this->costMatrix)) { initBiddersAndGoods(); + runAuctionRound(this->costMatrix); dataType cost = getMatchingDistance(this->costMatrix); diff --git a/core/base/common/welcomeMsg.inl b/core/base/common/welcomeMsg.inl index a3aae7ea9e..cdc132bdb5 100644 --- a/core/base/common/welcomeMsg.inl +++ b/core/base/common/welcomeMsg.inl @@ -2,53 +2,44 @@ // Julien Tierny // January 2020. -// "TTK (c) 2024" +// "TTK (c) 2025" printMsg( debug::output::BOLD - //+ " _____ _____ _ __ __ __ ____ ___ ____ _____" - + " _____ _____ _ __ __ __ ____ ___ ____ _ _" + + " _____ _____ _ __ __ __ ____ ___ ____ ____" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream); printMsg(debug::output::BOLD - //+ "|_ _|_ _| |/ / / /__\\ \\ |___ \\ / _ " - //"\\___ \\|___ /" - + "|_ _|_ _| |/ / / /__\\ \\ |___ \\ / _ " - "\\___ \\| || |" + + "|_ _|_ _| |/ / / /__\\ \\ |___ \\ / _ " + "\\___ \\| ___|" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream); printMsg( debug::output::BOLD - //+ " | | | | | ' / | |/ __| | __) | | | |__) | |_ - //\\" - + " | | | | | ' / | |/ __| | __) | | | |__) | || |_" + + " | | | | | ' / | |/ __| | __) | | | |__) |___ \\" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream); printMsg( debug::output::BOLD - // + " | | | | | . \\ | | (__| | / __/| |_| / __/ - // ___) |" - + " | | | | | . \\ | | (__| | / __/| |_| / __/|__ _|" + + " | | | | | . \\ | | (__| | / __/| |_| / __/ ___) |" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream); printMsg(debug::output::BOLD - // + " |_| |_| |_|\\_\\ | |\\___| | " - // "|_____|\\___/_____|____/" - + " |_| |_| |_|\\_\\ | |\\___| | " - "|_____|\\___/_____| |_|" + + " |_| |_| |_|\\_\\ | |\\___| | " + "|_____|\\___/_____|____/" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, stream); -printMsg(debug::output::BOLD + " \\_\\ /_/" +printMsg(debug::output::BOLD + " \\_\\ /_/" + debug::output::ENDCOLOR, debug::Priority::PERFORMANCE, debug::LineMode::NEW, diff --git a/core/base/persistenceDiagram/PersistenceDiagram.h b/core/base/persistenceDiagram/PersistenceDiagram.h index b7af9e5258..fe49c02f71 100644 --- a/core/base/persistenceDiagram/PersistenceDiagram.h +++ b/core/base/persistenceDiagram/PersistenceDiagram.h @@ -521,6 +521,7 @@ int ttk::PersistenceDiagram::executeDiscreteMorseSandwich( Timer const tm{}; const auto dim = triangulation->getDimensionality(); + dms_.setThreadNumber(this->threadNumber_); dms_.buildGradient( inputScalars, scalarsMTime, inputOffsets, *triangulation, updateMask); diff --git a/core/base/trackingFromCriticalPoints/CMakeLists.txt b/core/base/trackingFromCriticalPoints/CMakeLists.txt new file mode 100644 index 0000000000..5eb4252b2c --- /dev/null +++ b/core/base/trackingFromCriticalPoints/CMakeLists.txt @@ -0,0 +1,10 @@ +ttk_add_base_library(trackingFromCriticalPoints + SOURCES + TrackingFromCriticalPoints.cpp + HEADERS + TrackingFromCriticalPoints.h + DEPENDS + persistenceDiagram + triangulation + assignmentSolver + ) diff --git a/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.cpp b/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.cpp new file mode 100644 index 0000000000..33afb2dcc0 --- /dev/null +++ b/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.cpp @@ -0,0 +1,426 @@ +#include +#include +#include +#include + +// Compute L_p distance betweem (p,f(p)) and (q,f(q)) where p and q are critical +// points + +double ttk::TrackingFromCriticalPoints::criticalPointDistance( + const std::array &coords_p1, + const double &sfValue_p1, + const std::array &coords_p2, + const double &sfValue_p2, + const int &p = 2) { + return std::pow(xWeight * std::pow(coords_p1[0] - coords_p2[0], p) + + yWeight * std::pow(coords_p1[1] - coords_p2[1], p) + + zWeight * std::pow(coords_p1[2] - coords_p2[2], p) + + fWeight * std::pow(sfValue_p1 - sfValue_p2, p), + 1.0 / p); +} + +// Sort the critical points by types + +void ttk::TrackingFromCriticalPoints::sortCriticalPoint( + const DiagramType &d, + const double minimumRelevantPersistence, + std::vector> &maxCoords, + std::vector> &sad_1Coords, + std::vector> &sad_2Coords, + std::vector> &minCoords, + std::vector &maxScalar, + std::vector &sad_1Scalar, + std::vector &sad_2Scalar, + std::vector &minScalar, + std::vector &mapMax, + std::vector &mapSad_1, + std::vector &mapSad_2, + std::vector &mapMin) { + for(unsigned int i = 0; i < d.size(); i++) { + std::array birthCoords = d[i].birth.coords; + std::array deathCoords = d[i].death.coords; + if(std::abs(d[i].persistence()) > minimumRelevantPersistence) { + switch(d[i].dim) { + case 0: + minCoords.push_back(birthCoords); + minScalar.push_back(d[i].birth.sfValue); + mapMin.push_back(i); + + sad_1Coords.push_back(deathCoords); + sad_1Scalar.push_back(d[i].death.sfValue); + mapSad_1.push_back(i); + break; + case 1: + sad_2Coords.push_back(birthCoords); + sad_2Scalar.push_back(d[i].birth.sfValue); + mapSad_2.push_back(i); + + maxCoords.push_back(deathCoords); + maxScalar.push_back(d[i].death.sfValue); + mapMax.push_back(i); + break; + case 2: + sad_1Coords.push_back(birthCoords); + sad_1Scalar.push_back(d[i].birth.sfValue); + mapSad_1.push_back(i); + + sad_2Coords.push_back(deathCoords); + sad_2Scalar.push_back(d[i].death.sfValue); + mapSad_2.push_back(i); + break; + } + } + } +} + +void ttk::TrackingFromCriticalPoints::buildCostMatrix( + const std::vector> coords_1, + const std::vector sfValues_1, + const std::vector> coords_2, + const std::vector sfValues_2, + std::vector> &matrix, + float costDeathBirth) { + int size_1 = coords_1.size(); + int size_2 = coords_2.size(); + int matrix_size = (size_1 > 0 && size_2 > 0) ? size_1 + size_2 : 0; + for(int i = 0; i < size_1; i++) { + for(int j = 0; j < size_2; j++) { + matrix[i][j] = criticalPointDistance( + coords_1[i], sfValues_1[i], coords_2[j], sfValues_2[j]); + } + } + if(!adaptiveDeathBirthCost) { + for(int i = size_1; i < matrix_size; i++) { + for(int j = 0; j < size_2; j++) { + matrix[i][j] = costDeathBirth; + } + } + for(int i = 0; i < size_1; i++) { + for(int j = size_2; j < matrix_size; j++) { + matrix[i][j] = costDeathBirth; + } + } + } else if(adaptiveDeathBirthCost) { + for(int j = 0; j < size_2; j++) { + double c = matrix[0][j]; + for(int i = 1; i < size_1; i++) { + c = matrix[i][j] < c ? matrix[i][j] : c; + } + for(int i = size_1; i < matrix_size; i++) { + matrix[i][j] = c / (epsilonAdapt); + } + } + + for(int i = 0; i < size_1; i++) { + double d = matrix[i][0]; + for(int j = 1; j < size_2; j++) { + d = matrix[i][j] < d ? matrix[i][j] : d; + } + for(int j = size_2; j < matrix_size; j++) { + matrix[i][j] = d / (epsilonAdapt); + } + } + } +} + +void ttk::TrackingFromCriticalPoints::performMatchings( + const std::vector &persistenceDiagrams, + std::vector> &maximaMatchings, + std::vector> &sad_1_Matchings, + std::vector> &sad_2_Matchings, + std::vector> &minimaMatchings, + std::vector> &maxMap, + std::vector> &sad_1Map, + std::vector> &sad_2Map, + std::vector> &minMap) { + + int fieldNumber = persistenceDiagrams.size(); + + std::vector>> maxCoords(fieldNumber); + std::vector>> sad_1Coords(fieldNumber); + std::vector>> sad_2Coords(fieldNumber); + std::vector>> minCoords(fieldNumber); + + std::vector> maxScalar(fieldNumber); + std::vector> sad_1Scalar(fieldNumber); + std::vector> sad_2Scalar(fieldNumber); + std::vector> minScalar(fieldNumber); + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(int i = 0; i < fieldNumber; i++) { + + double minimumRelevantPersistence{}; + if(i < fieldNumber - 1) { + minimumRelevantPersistence + = ttk::TrackingFromCriticalPoints::computeRelevantPersistence( + persistenceDiagrams[i], persistenceDiagrams[i + 1]); + } + if(i == fieldNumber - 1) { + minimumRelevantPersistence + = ttk::TrackingFromCriticalPoints::computeRelevantPersistence( + persistenceDiagrams[i - 1], persistenceDiagrams[i]); + } + sortCriticalPoint(persistenceDiagrams[i], minimumRelevantPersistence, + maxCoords[i], sad_1Coords[i], sad_2Coords[i], + minCoords[i], maxScalar[i], sad_1Scalar[i], + sad_2Scalar[i], minScalar[i], maxMap[i], sad_1Map[i], + sad_2Map[i], minMap[i]); + } + +#ifdef TTK_ENABLE_OPENMP +#pragma omp parallel for num_threads(threadNumber_) +#endif // TTK_ENABLE_OPENMP + for(int i = 0; i < fieldNumber - 1; i++) { + + float costDeathBirth + = epsilonConstant + * computeBoundingBoxRadius( + persistenceDiagrams[i], persistenceDiagrams[i + 1]); + int maxSize = (maxCoords[i].size() > 0 && maxCoords[i + 1].size() > 0) + ? maxCoords[i].size() + maxCoords[i + 1].size() + : 0; + int sad_1Size = (sad_1Coords[i].size() > 0 && sad_1Coords[i + 1].size() > 0) + ? sad_1Coords[i].size() + sad_1Coords[i + 1].size() + : 0; + int sad_2Size = (sad_2Coords[i].size() > 0 && sad_2Coords[i + 1].size() > 0) + ? sad_2Coords[i].size() + sad_2Coords[i + 1].size() + : 0; + int minSize = (minCoords[i].size() > 0 && minCoords[i + 1].size() > 0) + ? minCoords[i].size() + minCoords[i + 1].size() + : 0; + + std::vector> maxMatrix( + maxSize, std::vector(maxSize, 0)); + std::vector> sad_1Matrix( + sad_1Size, std::vector(sad_1Size, 0)); + std::vector> sad_2Matrix( + sad_2Size, std::vector(sad_2Size, 0)); + std::vector> minMatrix( + minSize, std::vector(minSize, 0)); + + buildCostMatrix(maxCoords[i], maxScalar[i], maxCoords[i + 1], + maxScalar[i + 1], maxMatrix, costDeathBirth); + buildCostMatrix(sad_1Coords[i], sad_1Scalar[i], sad_1Coords[i + 1], + sad_1Scalar[i + 1], sad_1Matrix, costDeathBirth); + buildCostMatrix(sad_2Coords[i], sad_2Scalar[i], sad_2Coords[i + 1], + sad_2Scalar[i + 1], sad_2Matrix, costDeathBirth); + buildCostMatrix(minCoords[i], minScalar[i], minCoords[i + 1], + minScalar[i + 1], minMatrix, costDeathBirth); + + assignmentSolver(maxMatrix, maximaMatchings[i]); + assignmentSolver(sad_1Matrix, sad_1_Matchings[i]); + assignmentSolver(sad_2Matrix, sad_2_Matchings[i]); + assignmentSolver(minMatrix, minimaMatchings[i]); + } +} + +int ttk::TrackingFromCriticalPoints::computeGlobalId( + const DiagramType &persistenceDiagram, + const CriticalType &type, + const SimplexId &id) { + + switch(persistenceDiagram[id].dim) { + case 0: + return (type == CriticalType::Local_minimum + ? persistenceDiagram[id].birth.id + : persistenceDiagram[id].death.id); + case 1: + return (type == CriticalType::Saddle2 ? persistenceDiagram[id].birth.id + : persistenceDiagram[id].death.id); + case 2: + return (type == CriticalType::Saddle1 ? persistenceDiagram[id].birth.id + : persistenceDiagram[id].death.id); + } + return -1; +} + +void ttk::TrackingFromCriticalPoints::performTrackingForOneType( + const std::vector &persistenceDiagrams, + const std::vector> &matchings, + const std::vector> &map, + const CriticalType ¤tType, + std::vector &trackings, + std::vector> &trackingCosts, + std::vector &trackingPersistence) { + int fieldNumber = matchings.size() + 1; + + SimplexId deathLimitId = map[0].size(); + + std::vector previousStepMap(deathLimitId, -1); + std::vector sw; + + for(unsigned int i = 0; i < matchings[0].size(); i++) { + if(std::get<1>(matchings[0][i]) < deathLimitId) { + SimplexId startLocalId = std::get<0>(matchings[0][i]); + SimplexId endLocalId = std::get<1>(matchings[0][i]); + + SimplexId startId = computeGlobalId( + persistenceDiagrams[0], currentType, map[0][startLocalId]); + SimplexId endId = computeGlobalId( + persistenceDiagrams[1], currentType, map[1][endLocalId]); + std::vector chain = {startId, endId}; + std::tuple> tt + = std::make_tuple(0, 1, chain); + trackings.push_back(tt); + trackingPersistence.push_back( + persistenceDiagrams[0][map[0][startLocalId]].persistence() + + persistenceDiagrams[1][map[1][endLocalId]].persistence()); + std::vector newCostEntry; + newCostEntry.push_back(std::get<2>(matchings[0][i])); + + trackingCosts.push_back(newCostEntry); + + previousStepMap[std::get<1>(matchings[0][i])] = trackings.size() - 1; + } + } + for(int i = 1; i < fieldNumber - 1; i++) { + SimplexId birthLimitId = map[i].size(); + deathLimitId = map[i + 1].size(); + sw.resize(deathLimitId); + for(unsigned int j = 0; j < matchings[i].size(); j++) { + SimplexId startLocalId = std::get<0>(matchings[i][j]); + SimplexId endLocalId = std::get<1>(matchings[i][j]); + + bool wasPreviouslyMatched + = (startLocalId < birthLimitId ? previousStepMap[startLocalId] == -1 + : false); + bool existingTracking = endLocalId < deathLimitId && wasPreviouslyMatched; + + if(existingTracking) { + int trackingId = previousStepMap[startLocalId]; + SimplexId endId = computeGlobalId( + persistenceDiagrams[i + 1], currentType, map[i + 1][endLocalId]); + std::get<1>(trackings[trackingId])++; + std::get<2>(trackings[trackingId]).push_back(endId); + trackingCosts[trackingId].push_back(std::get<2>(matchings[i][j])); + trackingPersistence[trackingId] + += persistenceDiagrams[i + 1][map[i + 1][endLocalId]].persistence(); + sw[endLocalId] = trackingId; + } + + else if(startLocalId < birthLimitId && endLocalId < deathLimitId + && !wasPreviouslyMatched) { + + SimplexId startId = computeGlobalId( + persistenceDiagrams[i], currentType, map[i][startLocalId]); + SimplexId endId = computeGlobalId( + persistenceDiagrams[i + 1], currentType, map[i + 1][endLocalId]); + std::vector chain = {startId, endId}; + std::tuple> tt + = std::make_tuple(i, i + 1, chain); + trackings.push_back(tt); + std::vector newCostEntry; + newCostEntry.push_back(std::get<2>(matchings[i][j])); + trackingCosts.push_back(newCostEntry); + trackingPersistence.push_back( + persistenceDiagrams[i][map[i][startLocalId]].persistence() + + persistenceDiagrams[i + 1][map[i + 1][endLocalId]].persistence()); + sw[endLocalId] = trackings.size() - 1; + } + } + previousStepMap = sw; + sw.clear(); + } + for(unsigned int i = 0; i < trackings.size(); i++) { + trackingPersistence[i] + /= (std::get<1>(trackings[i]) - std::get<0>(trackings[i]) + 1); + } +} + +void ttk::TrackingFromCriticalPoints::performTrackings( + const std::vector &persistenceDiagrams, + const std::vector> &maximaMatchings, + const std::vector> &sad_1_Matchings, + const std::vector> &sad_2_Matchings, + const std::vector> &minimaMatchings, + const std::vector> &maxMap, + const std::vector> &sad_1Map, + const std::vector> &sad_2Map, + const std::vector> &minMap, + std::vector &allTrackings, + std::vector> &allTrackingsCosts, + std::vector &allTrackingsMeanPersistences, + unsigned int (&typesArrayLimits)[]) { + + std::vector trackingsMax; + std::vector trackingsSad_1; + std::vector trackingsSad_2; + std::vector trackingsMin; + + std::vector> maxTrackingCost; + std::vector> sad_1_TrackingCost; + std::vector> sad_2_TrackingCost; + std::vector> minTrackingCost; + + std::vector trackingsPersistenceMax; + std::vector trackingsPersistenceSad_1; + std::vector trackingsPersistenceSad_2; + std::vector trackingsPersistenceMin; + + performTrackingForOneType(persistenceDiagrams, maximaMatchings, maxMap, + CriticalType::Local_maximum, trackingsMax, + maxTrackingCost, trackingsPersistenceMax); + allTrackings.insert( + allTrackings.end(), trackingsMax.begin(), trackingsMax.end()); + allTrackingsMeanPersistences.insert(allTrackingsMeanPersistences.end(), + trackingsPersistenceMax.begin(), + trackingsPersistenceMax.end()); + allTrackingsCosts.insert( + allTrackingsCosts.end(), maxTrackingCost.begin(), maxTrackingCost.end()); + typesArrayLimits[0] = allTrackings.size(); + + performTrackingForOneType(persistenceDiagrams, sad_1_Matchings, sad_1Map, + CriticalType::Saddle1, trackingsSad_1, + sad_1_TrackingCost, trackingsPersistenceSad_1); + allTrackings.insert( + allTrackings.end(), trackingsSad_1.begin(), trackingsSad_1.end()); + allTrackingsMeanPersistences.insert(allTrackingsMeanPersistences.end(), + trackingsPersistenceSad_1.begin(), + trackingsPersistenceSad_1.end()); + allTrackingsCosts.insert(allTrackingsCosts.end(), sad_1_TrackingCost.begin(), + sad_1_TrackingCost.end()); + typesArrayLimits[1] = allTrackings.size(); + + performTrackingForOneType(persistenceDiagrams, sad_2_Matchings, sad_2Map, + CriticalType::Saddle2, trackingsSad_2, + sad_2_TrackingCost, trackingsPersistenceSad_2); + allTrackings.insert( + allTrackings.end(), trackingsSad_2.begin(), trackingsSad_2.end()); + allTrackingsMeanPersistences.insert(allTrackingsMeanPersistences.end(), + trackingsPersistenceSad_2.begin(), + trackingsPersistenceSad_2.end()); + allTrackingsCosts.insert(allTrackingsCosts.end(), sad_2_TrackingCost.begin(), + sad_2_TrackingCost.end()); + typesArrayLimits[2] = allTrackings.size(); + + performTrackingForOneType(persistenceDiagrams, minimaMatchings, minMap, + CriticalType::Local_minimum, trackingsMin, + minTrackingCost, trackingsPersistenceMin); + allTrackings.insert( + allTrackings.end(), trackingsMin.begin(), trackingsMin.end()); + allTrackingsMeanPersistences.insert(allTrackingsMeanPersistences.end(), + trackingsPersistenceMin.begin(), + trackingsPersistenceMin.end()); + allTrackingsCosts.insert( + allTrackingsCosts.end(), minTrackingCost.begin(), minTrackingCost.end()); +} + +void ttk::TrackingFromCriticalPoints::assignmentSolver( + std::vector> &costMatrix, + std::vector &matching) { + if(costMatrix.size() > 0) { + if(assignmentMethod == 0) { + ttk::AssignmentAuction solver; + solver.setInput(costMatrix); + solver.run(matching); + solver.clearMatrix(); + } else if(assignmentMethod == 1) { + ttk::AssignmentMunkres solver; + solver.setInput(costMatrix); + solver.run(matching); + solver.clearMatrix(); + } + } +} diff --git a/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.h b/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.h new file mode 100644 index 0000000000..8a6a65c35d --- /dev/null +++ b/core/base/trackingFromCriticalPoints/TrackingFromCriticalPoints.h @@ -0,0 +1,198 @@ +/// \ingroup base +/// \class ttk::TrackingFromPersistenceDiagrams +/// \author Thomas Daniel +/// \date January 2025. +/// +/// \b Online \b examples: \n +/// - Time +/// tracking example + +#pragma once + +// base code includes +#include +#include +#include + +namespace ttk { + + using trackingTuple = std::tuple>; + + class TrackingFromCriticalPoints : virtual public Debug { + + private: + double epsilonConstant{10e-1}; + double epsilonAdapt{0.5}; + double meshDiameter{1}; + double tolerance{10e-3}; + int assignmentMethod{0}; + double xWeight{1}; + double yWeight{1}; + double zWeight{1}; + double fWeight{1}; + bool adaptiveDeathBirthCost{false}; + + public: + TrackingFromCriticalPoints() { + } + + void setMeshDiameter(double r) { + meshDiameter = r; + } + + void setEpsilon(double e) { + epsilonConstant = e; + } + + void setEpsilonAdapt(double e) { + epsilonAdapt = e; + } + + void setTolerance(double t) { + tolerance = t; + } + + void setAssignmentMethod(int a) { + if(a == 0 || a == 1) { + assignmentMethod = a; + } + } + + void setAdaptDeathBirthCost(bool b) { + adaptiveDeathBirthCost = b; + } + + void setWeights(double PX, double PY, double PZ, double PF) { + xWeight = PX; + yWeight = PY; + zWeight = PZ; + fWeight = PF; + } + + double computeBoundingBoxRadius(const DiagramType &d1, + const DiagramType &d2) { + double maxScalar = d1[0].birth.sfValue; + double minScalar = d1[0].birth.sfValue; + + for(unsigned int i = 0; i < d1.size(); i++) { + maxScalar = std::max(maxScalar, d1[i].birth.sfValue); + maxScalar = std::max(maxScalar, d1[i].death.sfValue); + minScalar = std::min(minScalar, d1[i].birth.sfValue); + minScalar = std::min(minScalar, d1[i].death.sfValue); + } + + for(unsigned int i = 0; i < d2.size(); i++) { + maxScalar = std::max(maxScalar, d2[i].birth.sfValue); + maxScalar = std::max(maxScalar, d2[i].death.sfValue); + minScalar = std::min(minScalar, d2[i].birth.sfValue); + minScalar = std::min(minScalar, d2[i].death.sfValue); + } + + return std::sqrt(std::pow(meshDiameter, 2) + + std::pow(maxScalar - minScalar, 2)); + } + + void + performMatchings(const std::vector &persistenceDiagrams, + std::vector> &maximaMatchings, + std::vector> &sad_1_Matchings, + std::vector> &sad_2_Matchings, + std::vector> &minimaMatchings, + std::vector> &maxMap, + std::vector> &sad_1Map, + std::vector> &sad_2Map, + std::vector> &minMap); + void performTrackings( + const std::vector &persistenceDiagrams, + const std::vector> &maximaMatchings, + const std::vector> &sad_1_Matchings, + const std::vector> &sad_2_Matchings, + const std::vector> &minimaMatchings, + const std::vector> &maxMap, + const std::vector> &sad_1Map, + const std::vector> &sad_2Map, + const std::vector> &minMap, + std::vector &allTrackings, + std::vector> &allTrackingsCost, + std::vector &allTrackingsMeanPersistences, + unsigned int (&typesArrayLimits)[]); + + private: + double computeRelevantPersistence(const DiagramType &d1, + const DiagramType &d2) { + const auto sp = this->tolerance; + const double s = sp > 0.0 && sp < 100.0 ? sp / 100.0 : 0; + + std::vector toSort(d1.size() + d2.size()); + for(size_t i = 0; i < d1.size(); ++i) { + const auto &t = d1[i]; + toSort[i] = std::abs(t.persistence()); + } + for(size_t i = 0; i < d2.size(); ++i) { + const auto &t = d2[i]; + toSort[d1.size() + i] = std::abs(t.persistence()); + } + + const auto minVal = *std::min_element(toSort.begin(), toSort.end()); + const auto maxVal = *std::max_element(toSort.begin(), toSort.end()); + return s * (maxVal - minVal); + } + + // Compute L_p distance betweem (p,f(p)) and (q,f(q)) where p and q are + // critical points + + double criticalPointDistance(const std::array &coords_p1, + const double &sfValue_p1, + const std::array &coords_p2, + const double &sfValue_p2, + const int &p); + + // Sort the critical points by types + + void sortCriticalPoint(const DiagramType &d, + const double minimumRelevantPersistence, + std::vector> &maxCoords, + std::vector> &sad_1Coords, + std::vector> &sad_2Coords, + std::vector> &minCoords, + std::vector &maxScalar, + std::vector &sad1Scalar, + std::vector &sad_2Scalar, + std::vector &minScalar, + std::vector &mapMax, + std::vector &mapSad_1, + std::vector &mapSad_2, + std::vector &mapMin); + + void buildCostMatrix(const std::vector> coords_1, + const std::vector sfValues_1, + const std::vector> coords_2, + const std::vector sfValues_2, + std::vector> &matrix, + float costDeathBirth); + + void localToGlobalMatching(const std::vector &startMap, + const std::vector &endMap, + const std::vector &startPersistence, + const std::vector &endPersistence, + std::vector &matchings, + std::vector &matchingsPersistence); + + void assignmentSolver(std::vector> &costMatrix, + std::vector &matching); + + int computeGlobalId(const DiagramType &persistenceDiagram, + const CriticalType &type, + const SimplexId &id); + + void performTrackingForOneType( + const std::vector &persistenceDiagrams, + const std::vector> &matching, + const std::vector> &map, + const CriticalType ¤tType, + std::vector &tracking, + std::vector> &trackingCosts, + std::vector &trackingPersistence); + }; +} // namespace ttk diff --git a/core/base/trackingFromFields/CMakeLists.txt b/core/base/trackingFromFields/CMakeLists.txt index c22a0553b6..e97a9c353b 100644 --- a/core/base/trackingFromFields/CMakeLists.txt +++ b/core/base/trackingFromFields/CMakeLists.txt @@ -7,5 +7,6 @@ ttk_add_base_library(trackingFromFields persistenceDiagram bottleneckDistance trackingFromPersistenceDiagrams + trackingFromCriticalPoints triangulation ) diff --git a/core/base/trackingFromFields/TrackingFromFields.h b/core/base/trackingFromFields/TrackingFromFields.h index e1a7fd35eb..e35516c947 100644 --- a/core/base/trackingFromFields/TrackingFromFields.h +++ b/core/base/trackingFromFields/TrackingFromFields.h @@ -14,7 +14,7 @@ #include #include #include - +#include namespace ttk { class TrackingFromFields : virtual public Debug { @@ -91,26 +91,26 @@ int ttk::TrackingFromFields::performDiagramComputation( int fieldNumber, std::vector &persistenceDiagrams, const triangulationType *triangulation) { - + #ifdef TTK_ENABLE_OPENMP #pragma omp parallel for num_threads(threadNumber_) #endif // TTK_ENABLE_OPENMP for(int i = 0; i < fieldNumber; ++i) { ttk::PersistenceDiagram persistenceDiagram; persistenceDiagram.setThreadNumber(1); + persistenceDiagram.setBackend(ttk::PersistenceDiagram::BACKEND::DISCRETE_MORSE_SANDWICH); persistenceDiagram.execute(persistenceDiagrams[i], (dataType *)(inputData_[i]), 0, inputOffsets_[i], triangulation); - // Augment diagram. - for(auto &pair : persistenceDiagrams[i]) { - triangulation->getVertexPoint(pair.birth.id, pair.birth.coords[0], - pair.birth.coords[1], pair.birth.coords[2]); - triangulation->getVertexPoint(pair.death.id, pair.death.coords[0], - pair.death.coords[1], pair.death.coords[2]); - pair.birth.sfValue = static_cast(inputData_[i])[pair.birth.id]; - pair.death.sfValue = static_cast(inputData_[i])[pair.death.id]; - } + //for(auto &pair : persistenceDiagrams[i]) { + // triangulation->getVertexPoint(pair.birth.id, pair.birth.coords[0], + // pair.birth.coords[1], pair.birth.coords[2]); + // triangulation->getVertexPoint(pair.death.id, pair.death.coords[0], + // pair.death.coords[1], pair.death.coords[2]); + // pair.birth.sfValue = static_cast(inputData_[i])[pair.birth.id]; + // pair.death.sfValue = static_cast(inputData_[i])[pair.death.id]; + //} } return 0; diff --git a/core/base/trackingFromPersistenceDiagrams/TrackingFromPersistenceDiagrams.cpp b/core/base/trackingFromPersistenceDiagrams/TrackingFromPersistenceDiagrams.cpp index 8dd1c1b924..c44d552fc0 100644 --- a/core/base/trackingFromPersistenceDiagrams/TrackingFromPersistenceDiagrams.cpp +++ b/core/base/trackingFromPersistenceDiagrams/TrackingFromPersistenceDiagrams.cpp @@ -88,21 +88,24 @@ int ttk::TrackingFromPersistenceDiagrams::performTracking( std::vector &allDiagrams, std::vector> &allMatchings, std::vector &trackings) { + auto numPersistenceDiagramsInput = (int)allDiagrams.size(); - for(int in = 1; in < numPersistenceDiagramsInput - 1; ++in) { + + std::vector matchings1 = allMatchings[in - 1]; std::vector matchings2 = allMatchings[in]; auto matchingsSize1 = (int)matchings1.size(); auto matchingsSize2 = (int)matchings2.size(); + int const endIndex = numPersistenceDiagramsInput - 2; for(int i = 0; i < matchingsSize1; ++i) { const auto m1ai0 = std::get<0>(matchings1[i]); const auto m1ai1 = std::get<1>(matchings1[i]); - for(int j = 0; j < matchingsSize2; ++j) { + const auto m2aj0 = std::get<0>(matchings2[j]); const auto m2aj1 = std::get<1>(matchings2[j]); diff --git a/core/vtk/ttkAlgorithm/ttkAlgorithm.cpp b/core/vtk/ttkAlgorithm/ttkAlgorithm.cpp index d816394a51..6e69ee217e 100644 --- a/core/vtk/ttkAlgorithm/ttkAlgorithm.cpp +++ b/core/vtk/ttkAlgorithm/ttkAlgorithm.cpp @@ -104,6 +104,7 @@ vtkDataArray * vtkDataArray *oldOrderArray, ttk::Triangulation *triangulation) { + vtkSmartPointer newOrderArray; auto nVertices = scalarArray->GetNumberOfTuples(); if(oldOrderArray != nullptr && getGlobalOrder) { @@ -148,12 +149,14 @@ vtkDataArray * ->AddArray(newOrderArray); } #else + switch(scalarArray->GetDataType()) { vtkTemplateMacro(ttk::preconditionOrderArray( nVertices, static_cast(ttkUtils::GetVoidPointer(scalarArray)), static_cast(ttkUtils::GetVoidPointer(newOrderArray)), this->threadNumber_)); } + inputData ->GetAttributesAsFieldData( this->GetInputArrayAssociation(scalarArrayIdx, inputData)) diff --git a/core/vtk/ttkTrackingFromFields/ttk.module b/core/vtk/ttkTrackingFromFields/ttk.module index 08d1e2a0e9..6dbaffe55a 100644 --- a/core/vtk/ttkTrackingFromFields/ttk.module +++ b/core/vtk/ttkTrackingFromFields/ttk.module @@ -7,4 +7,5 @@ HEADERS DEPENDS trackingFromFields trackingFromPersistenceDiagrams + trackingFromCriticalPoints ttkAlgorithm diff --git a/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.cpp b/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.cpp index a0ef30e14b..28f824975b 100644 --- a/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.cpp +++ b/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.cpp @@ -38,15 +38,11 @@ int ttkTrackingFromFields::trackWithPersistenceMatching( unsigned long fieldNumber, const triangulationType *triangulation) { - using trackingTuple = ttk::trackingTuple; - - // 1. get persistence diagrams. std::vector persistenceDiagrams(fieldNumber); this->performDiagramComputation( (int)fieldNumber, persistenceDiagrams, triangulation); - // 2. call feature tracking with threshold. std::vector> outputMatchings(fieldNumber - 1); double const spacing = Spacing; @@ -82,8 +78,7 @@ int ttkTrackingFromFields::trackWithPersistenceMatching( componentIds->SetName("ConnectedComponentId"); pointTypeScalars->SetName("CriticalType"); - // (+ vertex id) - std::vector trackingsBase; + std::vector trackingsBase; tfp.performTracking(persistenceDiagrams, outputMatchings, trackingsBase); std::vector> trackingTupleToMerged(trackingsBase.size()); @@ -107,13 +102,116 @@ int ttkTrackingFromFields::trackWithPersistenceMatching( return 1; } +template +int ttkTrackingFromFields::trackWithCriticalPointMatching( + vtkUnstructuredGrid *output, + unsigned long fieldNumber, + const triangulationType *triangulation) { + + float x, y, z; + float maxX, minX, maxY, minY, maxZ, minZ; + triangulation->getVertexPoint(0, minX, minY, minZ); + triangulation->getVertexPoint(0, maxX, maxY, maxZ); + + for(int i = 0; i < triangulation->getNumberOfVertices(); i++) { + triangulation->getVertexPoint(i, x, y, z); + maxX = std::max(x, maxX); + maxX = std::min(x, minX); + maxY = std::max(y, maxX); + minY = std::min(y, minY); + maxZ = std::max(z, maxZ); + minZ = std::min(z, minZ); + } + + double const costDeathBirth = CostDeathBirth; + double const tolerance = (double)Tolerance; + float meshDiameter + = std::sqrt(std::pow(maxX - minX, 2) + std::pow(maxY - minY, 2) + + std::pow(maxZ - minZ, 2)); + int assignmentMethod = AssignmentMethod; + bool adaptDeathBirthCost = AdaptDeathBirthCost; + double epsilonAdapt = EpsilonAdapt; + + ttk::TrackingFromCriticalPoints tracker; + tracker.setMeshDiameter(meshDiameter); + tracker.setTolerance(tolerance); + tracker.setEpsilon(costDeathBirth); + tracker.setAdaptDeathBirthCost(adaptDeathBirthCost); + tracker.setAssignmentMethod(assignmentMethod); + tracker.setEpsilonAdapt(epsilonAdapt); + tracker.setWeights(PX, PY, PZ, PF); + + tracker.setThreadNumber(this->threadNumber_); + + std::vector persistenceDiagrams(fieldNumber); + this->performDiagramComputation( + (int)fieldNumber, persistenceDiagrams, triangulation); + + std::vector> maximaMatchings(fieldNumber - 1); + std::vector> sad_1_Matchings(fieldNumber - 1); + std::vector> sad_2_Matchings(fieldNumber - 1); + std::vector> minimaMatchings(fieldNumber - 1); + + std::vector> maxMap(fieldNumber); + std::vector> sad_1Map(fieldNumber); + std::vector> sad_2Map(fieldNumber); + std::vector> minMap(fieldNumber); + + tracker.performMatchings(persistenceDiagrams, maximaMatchings, + sad_1_Matchings, sad_2_Matchings, minimaMatchings, + maxMap, sad_1Map, sad_2Map, minMap); + + vtkNew const points{}; + vtkNew const outputMesh{}; + + vtkNew costs{}; + vtkNew averagePersistences{}; + vtkNew valueScalars{}; + vtkNew globalVertexIds{}; + vtkNew lengthScalars{}; + vtkNew timeScalars{}; + vtkNew connectedComponentIds{}; + vtkNew pointsCriticalType{}; + + costs->SetName("Costs"); + averagePersistences->SetName("AveragePersistence"); + valueScalars->SetName("Scalar"); + globalVertexIds->SetName("VertexGlobalId"); + lengthScalars->SetName("ComponentLength"); + timeScalars->SetName("TimeStep"); + connectedComponentIds->SetName("ConnectedComponentId"); + pointsCriticalType->SetName("CriticalType"); + + std::vector allTrackings; + std::vector allTrackingsMeanPersistence; + std::vector> allTrackingsCosts; + + unsigned int typesArrayLimits[3] = {}; + tracker.performTrackings( + persistenceDiagrams, maximaMatchings, sad_1_Matchings, sad_2_Matchings, + minimaMatchings, maxMap, sad_1Map, sad_2Map, minMap, allTrackings, + allTrackingsCosts, allTrackingsMeanPersistence, typesArrayLimits); + + double const spacing = Spacing; + bool const useGeometricSpacing = UseGeometricSpacing; + + ttkTrackingFromPersistenceDiagrams::buildMesh( + triangulation, allTrackings, allTrackingsCosts, allTrackingsMeanPersistence, + useGeometricSpacing, spacing, points, outputMesh, pointsCriticalType, + timeScalars, lengthScalars, globalVertexIds, connectedComponentIds, costs, + averagePersistences, typesArrayLimits); + + output->ShallowCopy(outputMesh); + + return 1; +} + int ttkTrackingFromFields::RequestData(vtkInformation *ttkNotUsed(request), vtkInformationVector **inputVector, vtkInformationVector *outputVector) { auto input = vtkDataSet::GetData(inputVector[0]); auto output = vtkUnstructuredGrid::GetData(outputVector); - ttk::Triangulation *triangulation = ttkAlgorithm::GetTriangulation(input); if(!triangulation) return 0; @@ -177,6 +275,7 @@ int ttkTrackingFromFields::RequestData(vtkInformation *ttkNotUsed(request), std::string const algorithm = DistanceAlgorithm; int const pvalg = PVAlgorithm; bool useTTKMethod = false; + bool trackWithCriticalPoints = (pvalg == 2); if(pvalg >= 0) { switch(pvalg) { @@ -217,7 +316,7 @@ int ttkTrackingFromFields::RequestData(vtkInformation *ttkNotUsed(request), // 0. get data int const fieldNumber = inputScalarFields.size(); std::vector inputFields(fieldNumber); - for(int i = 0; i < fieldNumber; ++i) { + for(int i = 0; i < fieldNumber; i++) { inputFields[i] = ttkUtils::GetVoidPointer(inputScalarFields[i]); } this->setInputScalars(inputFields); @@ -234,11 +333,16 @@ int ttkTrackingFromFields::RequestData(vtkInformation *ttkNotUsed(request), this->setInputOffsets(inputOrders); int status = 0; - if(useTTKMethod) { + if(useTTKMethod && !trackWithCriticalPoints) { ttkVtkTemplateMacro( inputScalarFields[0]->GetDataType(), triangulation->getType(), (status = this->trackWithPersistenceMatching( output, fieldNumber, (TTK_TT *)triangulation->getData()))); + } else if(useTTKMethod && trackWithCriticalPoints) { + ttkVtkTemplateMacro( + inputScalarFields[0]->GetDataType(), triangulation->getType(), + (status = this->trackWithCriticalPointMatching( + output, fieldNumber, (TTK_TT *)triangulation->getData()))); } else { this->printMsg("The specified matching method is not supported."); } diff --git a/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.h b/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.h index 36f497c9cf..bb20e3abf7 100644 --- a/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.h +++ b/core/vtk/ttkTrackingFromFields/ttkTrackingFromFields.h @@ -30,10 +30,12 @@ #pragma once +#include #include #include // VTK Module +#include #include #include #include @@ -73,8 +75,17 @@ class TTKTRACKINGFROMFIELDS_EXPORT ttkTrackingFromFields /// @{ vtkSetMacro(Tolerance, double); vtkGetMacro(Tolerance, double); + + vtkSetMacro(CostDeathBirth, double); + vtkGetMacro(CostDeathBirth, double); /// @} + vtkSetMacro(AdaptDeathBirthCost, bool); + vtkGetMacro(AdaptDeathBirthCost, bool); + + vtkSetMacro(EpsilonAdapt, double); + vtkGetMacro(EpsilonAdapt, double); + /// @brief Importance weight for the X component of the extremum. /// @{ vtkSetMacro(PX, double); @@ -105,6 +116,12 @@ class TTKTRACKINGFROMFIELDS_EXPORT ttkTrackingFromFields vtkGetMacro(PS, double); /// @} + /// @brief Importance weight for function values. + /// @{ + vtkSetMacro(PF, double); + vtkGetMacro(PF, double); + /// @} + /// @brief Value of the parameter p for the Wp (p-th Wasserstein) distance /// computation (type "inf" for the Bottleneck distance). /// @{ @@ -121,6 +138,9 @@ class TTKTRACKINGFROMFIELDS_EXPORT ttkTrackingFromFields /// @{ vtkSetMacro(PVAlgorithm, int); vtkGetMacro(PVAlgorithm, int); + + vtkSetMacro(AssignmentMethod, int); + vtkGetMacro(AssignmentMethod, int); /// @} /// @brief For the translation of the second set of critical points even the @@ -170,9 +190,15 @@ class TTKTRACKINGFROMFIELDS_EXPORT ttkTrackingFromFields double Tolerance{1}; double PX{1}; double PY{1}; - double PZ{0}; + double PZ{1}; double PE{0}; double PS{0}; + double PF{1}; + + double CostDeathBirth{0.1}; + double EpsilonAdapt{0.5}; + int AssignmentMethod{0}; + bool AdaptDeathBirthCost{false}; // Bottleneck config. bool UseGeometricSpacing{false}; @@ -187,4 +213,9 @@ class TTKTRACKINGFROMFIELDS_EXPORT ttkTrackingFromFields int trackWithPersistenceMatching(vtkUnstructuredGrid *output, unsigned long fieldNumber, const triangulationType *triangulation); + + template + int trackWithCriticalPointMatching(vtkUnstructuredGrid *output, + unsigned long fieldNumber, + const triangulationType *triangulation); }; diff --git a/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.cpp b/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.cpp index 97a2e1e2d0..ff0c20a24d 100644 --- a/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.cpp +++ b/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.cpp @@ -2,15 +2,6 @@ #include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include vtkStandardNewMacro(ttkTrackingFromPersistenceDiagrams); @@ -55,7 +46,6 @@ int ttkTrackingFromPersistenceDiagrams::buildMesh( vtkIntArray *componentIds, vtkIntArray *pointTypeScalars, const ttk::Debug &dbg) { - using ttk::CriticalType; int currentVertex = 0; for(size_t k = 0; k < trackings.size(); ++k) { @@ -71,6 +61,8 @@ int ttkTrackingFromPersistenceDiagrams::buildMesh( return 0; } + + for(int c = 0; c < chainLength - 1; ++c) { const auto &matchings1 = outputMatchings[numStart + c]; int const d1id = numStart + c; @@ -135,6 +127,7 @@ int ttkTrackingFromPersistenceDiagrams::buildMesh( z1 += spacing * (numStart + c); } + // Postproc component ids. int cid = k; bool hasMergedFirst = false; @@ -271,6 +264,7 @@ int ttkTrackingFromPersistenceDiagrams::buildMesh( lengthScalars->InsertTuple1(currentVertex, chainLength); currentVertex++; + } } @@ -286,6 +280,7 @@ int ttkTrackingFromPersistenceDiagrams::buildMesh( return 0; } + int ttkTrackingFromPersistenceDiagrams::RequestData( vtkInformation *ttkNotUsed(request), vtkInformationVector **inputVector, diff --git a/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.h b/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.h index b18738372d..f9f5210c85 100644 --- a/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.h +++ b/core/vtk/ttkTrackingFromPersistenceDiagrams/ttkTrackingFromPersistenceDiagrams.h @@ -5,6 +5,15 @@ // VTK Module #include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include class vtkUnstructuredGrid; @@ -77,6 +86,77 @@ class TTKTRACKINGFROMPERSISTENCEDIAGRAMS_EXPORT vtkIntArray *pointTypeScalars, const ttk::Debug &dbg); +template + static int buildMesh( + const triangulationType *triangulation, + const std::vector &trackings, + const std::vector> &allTrackingsCosts, + const std::vector &allTrackingsPersistence, + const bool useGeometricSpacing, + const double spacing, + vtkPoints *points, + vtkUnstructuredGrid *outputMesh, + vtkIntArray *pointsCriticalType, + vtkIntArray * timeScalars, + vtkIntArray *lengthScalars, + vtkIntArray *globalVertexIds, + vtkIntArray *connectedComponentIds, + vtkDoubleArray *costs, + vtkDoubleArray *averagePersistence, + unsigned int *sizes){ + + int pointCpt = 0; + int edgeCpt=0; + for (unsigned int i = 0 ; i < trackings.size() ; i++){ + ttk::CriticalType currentType=ttk::CriticalType::Local_minimum; + if(i < sizes[0])currentType=ttk::CriticalType::Local_maximum; + else if(i < sizes[1] && i >=sizes[0])currentType = ttk::CriticalType::Saddle1; + else if(i < sizes[2] && i >=sizes[1])currentType = ttk::CriticalType::Saddle2; + int startTime = std::get<0>(trackings[i]); + std::vector chain = std::get<2>(trackings[i]); + + float x=0; + float y=0; + float z=0; + triangulation->getVertexPoint(chain[0], x, y, z); + if(useGeometricSpacing)z+=startTime *spacing; + points->InsertNextPoint(x,y,z); + globalVertexIds->InsertTuple1(pointCpt, (int)chain[0]); + pointsCriticalType->InsertTuple1(pointCpt, (int)currentType); + timeScalars->InsertTuple1(pointCpt, startTime); + vtkIdType edge[2]; + for (unsigned int j = 1 ; j < chain.size() ; j++){ + triangulation->getVertexPoint(chain[j], x, y, z); + if(useGeometricSpacing)z+=(j+startTime)*spacing; + edge[0]=pointCpt; + pointCpt++; + edge[1]=pointCpt; + points->InsertNextPoint(x, y, z); + globalVertexIds->InsertTuple1(pointCpt, (int)chain[j]); + outputMesh->InsertNextCell(VTK_LINE, 2, edge); + pointsCriticalType->InsertTuple1(pointCpt, (int)currentType); + timeScalars->InsertTuple1(pointCpt, startTime+j); + lengthScalars->InsertTuple1(edgeCpt, chain.size()-1); + connectedComponentIds->InsertTuple1(edgeCpt, i); + costs->InsertTuple1(edgeCpt, allTrackingsCosts[i][j-1]); + averagePersistence->InsertTuple1(edgeCpt, allTrackingsPersistence[i]); + edgeCpt++; + } + pointCpt++; + } + + outputMesh->SetPoints(points); + outputMesh->GetCellData()->AddArray(lengthScalars); + outputMesh->GetCellData()->AddArray(connectedComponentIds); + outputMesh->GetCellData()->AddArray(averagePersistence); + outputMesh->GetCellData()->AddArray(costs); + outputMesh->GetPointData()->AddArray(pointsCriticalType); + outputMesh->GetPointData()->AddArray(timeScalars); + outputMesh->GetPointData()->AddArray(globalVertexIds); + + return 0; + } + protected: ttkTrackingFromPersistenceDiagrams(); diff --git a/paraview/patch/headless.yml b/paraview/patch/headless.yml index 0692628233..39186a8b99 100644 --- a/paraview/patch/headless.yml +++ b/paraview/patch/headless.yml @@ -72,7 +72,7 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-12, macos-14] + os: [macos-13, macos-14, macos-15] env: CCACHE_DIR: /Users/runner/work/ttk/.ccache CCACHE_MAXSIZE: 200M @@ -80,15 +80,13 @@ jobs: - uses: actions/checkout@v4 name: Checkout TTK-ParaView source code - - name: Remove hosted Python - run: | - sudo rm -rf /Library/Frameworks/Python.framework/ - sudo rm -rf /usr/local/Frameworks/Python.framework/ + - uses: actions/setup-python@v5 + with: + python-version: '3.12' - name: Install macOS dependencies run: | # ParaView dependencies - brew reinstall python brew install --cask xquartz brew install ninja open-mpi @@ -97,7 +95,7 @@ jobs: mkdir build && cd build cmake \ -DCMAKE_BUILD_TYPE=Release \ - -DPARAVIEW_USE_MPI=ON \ + -DPARAVIEW_USE_MPI=OFF \ -DPARAVIEW_USE_QT=OFF \ -GNinja \ $GITHUB_WORKSPACE @@ -112,6 +110,22 @@ jobs: cd build cpack -G TGZ + - name: Install ParaView + run: | + cd build + sudo cmake --build . --target install + # pvpython does not embed the correct PYTHONPATH + echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV + + - name: Test Python imports + run: | + python3 -m vtk + python3 -m paraview.simple + pvpython -m vtk + pvpython -m paraview.simple + env: + DYLD_LIBRARY_PATH: /usr/local/lib + - name: Upload compressed binaries uses: actions/upload-artifact@v4 with: @@ -240,22 +254,29 @@ jobs: file: ttk-paraview-headless-ubuntu-24.04/ttk-paraview.deb asset_name: ttk-paraview-headless-ubuntu-24.04.deb + - name: Upload MacOS 15 .tar.gz as Release Asset + uses: svenstaro/upload-release-action@v2 + with: + repo_token: ${{ secrets.GITHUB_TOKEN }} + tag: ${{ github.ref }} + file: ttk-paraview-headless-macos-15/ttk-paraview.tar.gz + asset_name: ttk-paraview-headless-macos-15.tar.gz - - name: Upload MacOS 14 (arm64) .tar.gz as Release Asset + - name: Upload MacOS 14 .tar.gz as Release Asset uses: svenstaro/upload-release-action@v2 with: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} file: ttk-paraview-headless-macos-14/ttk-paraview.tar.gz - asset_name: ttk-paraview-headless-macos-14-arm64.tar.gz + asset_name: ttk-paraview-headless-macos-14.tar.gz - - name: Upload MacOS 12 .tar.gz as Release Asset + - name: Upload MacOS 13 .tar.gz as Release Asset uses: svenstaro/upload-release-action@v2 with: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} - file: ttk-paraview-headless-macos-12/ttk-paraview.tar.gz - asset_name: ttk-paraview-headless-macos-12.tar.gz + file: ttk-paraview-headless-macos-13/ttk-paraview.tar.gz + asset_name: ttk-paraview-headless-macos-13.tar.gz - name: Upload .exe as Release Asset diff --git a/paraview/patch/package.yml b/paraview/patch/package.yml index 52f150e365..846a69df7b 100644 --- a/paraview/patch/package.yml +++ b/paraview/patch/package.yml @@ -75,20 +75,21 @@ jobs: runs-on: ${{ matrix.os }} strategy: matrix: - os: [macos-12, macos-14] + os: [macos-13, macos-14, macos-15] + env: + CCACHE_DIR: /Users/runner/work/ttk/.ccache + CCACHE_MAXSIZE: 200M steps: - uses: actions/checkout@v4 name: Checkout TTK-ParaView source code - - name: Remove hosted Python - run: | - sudo rm -rf /Library/Frameworks/Python.framework/ - sudo rm -rf /usr/local/Frameworks/Python.framework + - uses: actions/setup-python@v5 + with: + python-version: '3.12' - name: Install macOS dependencies run: | # ParaView dependencies - brew reinstall python brew install --cask xquartz brew install mesa glew qt@5 ninja @@ -98,7 +99,6 @@ jobs: cmake \ -DCMAKE_BUILD_TYPE=Release \ -DQt5_DIR=$(brew --prefix qt@5)/lib/cmake/Qt5 \ - -DPython3_ROOT_DIR=$(brew --prefix python) \ -GNinja \ $GITHUB_WORKSPACE @@ -112,6 +112,22 @@ jobs: cd build cpack -G TGZ + - name: Install ParaView + run: | + cd build + sudo cmake --build . --target install + # pvpython does not embed the correct PYTHONPATH + echo "PYTHONPATH=/usr/local/lib/python3.12/site-packages:$PYTHONPATH" >> $GITHUB_ENV + + - name: Test Python imports + run: | + python3 -m vtk + python3 -m paraview.simple + pvpython -m vtk + pvpython -m paraview.simple + env: + DYLD_LIBRARY_PATH: /usr/local/lib + - name: Upload compressed binaries uses: actions/upload-artifact@v4 with: @@ -239,21 +255,29 @@ jobs: file: ttk-paraview-ubuntu-24.04/ttk-paraview.deb asset_name: ttk-paraview-$tag-ubuntu-24.04.deb - - name: Upload MacOS 14 (arm64) .tar.gz as Release Asset + - name: Upload MacOS 15 .tar.gz as Release Asset + uses: svenstaro/upload-release-action@v2 + with: + repo_token: ${{ secrets.GITHUB_TOKEN }} + tag: ${{ github.ref }} + file: ttk-paraview-macos-15/ttk-paraview.tar.gz + asset_name: ttk-paraview-$tag-macos-15.tar.gz + + - name: Upload MacOS 14 .tar.gz as Release Asset uses: svenstaro/upload-release-action@v2 with: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} file: ttk-paraview-macos-14/ttk-paraview.tar.gz - asset_name: ttk-paraview-$tag-macos-14-arm64.tar.gz + asset_name: ttk-paraview-$tag-macos-14.tar.gz - - name: Upload MacOS 12 .tar.gz as Release Asset + - name: Upload MacOS 13 .tar.gz as Release Asset uses: svenstaro/upload-release-action@v2 with: repo_token: ${{ secrets.GITHUB_TOKEN }} tag: ${{ github.ref }} - file: ttk-paraview-macos-12/ttk-paraview.tar.gz - asset_name: ttk-paraview-$tag-macos-12.tar.gz + file: ttk-paraview-macos-13/ttk-paraview.tar.gz + asset_name: ttk-paraview-$tag-macos-13.tar.gz - name: Upload .exe as Release Asset uses: svenstaro/upload-release-action@v2 diff --git a/paraview/xmls/TrackingFromFields.xml b/paraview/xmls/TrackingFromFields.xml index a29e2ebab1..9d78f85480 100644 --- a/paraview/xmls/TrackingFromFields.xml +++ b/paraview/xmls/TrackingFromFields.xml @@ -94,25 +94,45 @@ by a list of scalar fields) and which computes a tracking mesh." - - - - + number_of_elements="1" + default_values="0" + > - - + + + + Method for computing matchings. - + + + + + + + + + + Assignment algorithm used betweem each time step. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Importance weight for extrema. @@ -145,6 +245,20 @@ by a list of scalar fields) and which computes a tracking mesh." default_values="1" panel_visibility="advanced"> + + + + + + + + Importance weight for saddles. @@ -153,7 +267,7 @@ by a list of scalar fields) and which computes a tracking mesh." label="X weight" command="SetPX" number_of_elements="1" - default_values="0" + default_values="1" panel_visibility="advanced"> Importance weight for the X component of @@ -165,7 +279,7 @@ the extremum label="Y weight" command="SetPY" number_of_elements="1" - default_values="0" + default_values="1" panel_visibility="advanced"> Importance weight for the Y component of @@ -177,13 +291,30 @@ the extremum label="Z weight" command="SetPZ" number_of_elements="1" - default_values="0" + default_values="1" panel_visibility="advanced"> Importance weight for the Z component of the extremum + + + + + + Importance weight for the function value of the extremum + +