DiscreteGradient.h

Go to the documentation of this file.

 
#pragma once
 
// base code includes
#include <FTMTree.h>
#include <Geometry.h>
#include <Triangulation.h>
#include <VisitedMask.h>
#include <algorithm>
#include <array>
#include <functional>
#include <queue>
#include <set>
#include <utility>
 
namespace ttk {
 
  namespace dcg {
    struct Cell {
      explicit Cell() = default;
 
      explicit Cell(const int dim, const SimplexId id) : dim_{dim}, id_{id} {
      }
 
      inline bool operator==(const Cell &other) const {
        return std::tie(this->dim_, this->id_)
               == std::tie(other.dim_, other.id_);
      }
 
      inline std::string to_string() const {
        return '{' + std::to_string(this->dim_) + ' '
               + std::to_string(this->id_) + '}';
      }
 
      int dim_{-1};
      SimplexId id_{-1};
    };
 
    enum gradientValue { NULL_GRADIENT = -1, GHOST_GRADIENT = -2 };

    struct CellExt : Cell {
      explicit CellExt(const int dim, const SimplexId id) : Cell{dim, id} {
      }
      explicit CellExt(const int dim,
                       const SimplexId id,
                       const std::array<SimplexId, 3> &lowVerts,
                       const std::array<uint8_t, 3> &faces)
        : Cell{dim, id}, lowVerts_{lowVerts}, faces_{faces} {
      }
 
      // (order field value on) lower vertices in current lower star
      // (1 for edges, 2 for triangles, 3 for tetras)
      const std::array<SimplexId, 3> lowVerts_{};
      // indices of faces (cells of dimensions dim_ - 1) in lower star
      // structure, only applicable for triangles (2 edge faces)  and tetras (3
      // triangle faces)
      const std::array<uint8_t, 3> faces_{};
      // if cell has been paired with another in current lower star
      bool paired_{false};
    };

    struct VPath {
      explicit VPath() = default;
 
      explicit VPath(const bool isValid,
                     const SimplexId segmentId,
                     const SimplexId source,
                     const SimplexId destination,
                     const SimplexId sourceSlot,
                     const SimplexId destinationSlot,
                     const double persistence)
        : isValid_{isValid}, states_{1}, segments_{segmentId}, source_{source},
          destination_{destination}, sourceSlot_{sourceSlot},
          destinationSlot_{destinationSlot}, persistence_{persistence} {
      }
 
      explicit VPath(const bool isValid,
                     const std::vector<char> &states,
                     const std::vector<SimplexId> &segments,
                     const SimplexId source,
                     const SimplexId destination,
                     const SimplexId sourceSlot,
                     const SimplexId destinationSlot,
                     const double persistence)
        : isValid_{isValid}, states_{states}, segments_{segments},
          source_{source}, destination_{destination}, sourceSlot_{sourceSlot},
          destinationSlot_{destinationSlot}, persistence_{persistence} {
      }
 
      explicit VPath(const bool isValid,
                     std::vector<char> &&states,
                     std::vector<SimplexId> &&segments,
                     const SimplexId source,
                     const SimplexId destination,
                     const SimplexId sourceSlot,
                     const SimplexId destinationSlot,
                     const double persistence)
        : isValid_{isValid}, states_{states}, segments_{segments},
          source_{source}, destination_{destination}, sourceSlot_{sourceSlot},
          destinationSlot_{destinationSlot}, persistence_{persistence} {
      }

      int invalidate() {
        isValid_ = false;
        source_ = -1;
        destination_ = -1;
        persistence_ = -1;
        clear();
 
        return 0;
      }

      int clear() {
        states_.clear();
        segments_.clear();
 
        return 0;
      }
 
      bool isValid_{};
      std::vector<char> states_{};
      std::vector<SimplexId> segments_{};
      SimplexId source_{-1};
      SimplexId destination_{-1};
      SimplexId sourceSlot_{-1};
      SimplexId destinationSlot_{-1};
      double persistence_{};
    };

    struct CriticalPoint {
      explicit CriticalPoint() = default;
 
      explicit CriticalPoint(const Cell &cell) : cell_{cell}, numberOfSlots_{} {
      }
 
      explicit CriticalPoint(const Cell &cell,
                             const std::vector<SimplexId> &vpaths)
        : cell_{cell}, vpaths_{vpaths}, numberOfSlots_{} {
      }
 
      explicit CriticalPoint(const Cell &cell, std::vector<SimplexId> &&vpaths)
        : cell_{cell}, vpaths_{vpaths}, numberOfSlots_{} {
      }

      SimplexId omp_addSlot() {
        SimplexId numberOfSlots = 0;
 
#ifdef TTK_ENABLE_OPENMP
#pragma omp atomic capture
#endif
        numberOfSlots = (numberOfSlots_++);
 
        return numberOfSlots;
      }

      SimplexId addSlot() {
        return (numberOfSlots_++);
      }

      int clear() {
        vpaths_.clear();
 
        return 0;
      }
 
      Cell cell_{};
      std::vector<SimplexId> vpaths_{};
      SimplexId numberOfSlots_{};
    };

    template <typename dataType>
    struct SaddleSaddleVPathComparator {
      bool
        operator()(const std::tuple<dataType, SimplexId, SimplexId> &v1,
                   const std::tuple<dataType, SimplexId, SimplexId> &v2) const {
        const dataType persistence1 = std::get<0>(v1);
        const dataType persistence2 = std::get<0>(v2);
 
        const SimplexId vpathId1 = std::get<1>(v1);
        const SimplexId vpathId2 = std::get<1>(v2);
 
        const SimplexId saddleId1 = std::get<2>(v1);
        const SimplexId saddleId2 = std::get<2>(v2);
 
        if(persistence1 != persistence2) {
          return (persistence1 < persistence2);
        }
 
        if(saddleId1 != saddleId2) {
          return (saddleId1 < saddleId2);
        }
 
        return (vpathId1 < vpathId2);
      }
    };

    class DiscreteGradient : virtual public Debug {
 
    public:
      enum class BACKEND {
        CLASSIC_BACKEND = 0,
        STOCHASTIC_BACKEND = 1,
      };
 
      DiscreteGradient() {
        this->setDebugMsgPrefix("DiscreteGradient");
#ifdef TTK_ENABLE_MPI
        hasMPISupport_ = true;
#endif
      }

      template <typename dataType, typename triangulationType>
      int buildGradient(const triangulationType &triangulation,
                        bool bypassCache = false,
                        const std::vector<bool> *updateMask = nullptr);
 
      template <typename triangulationType>
      int buildGradient(const triangulationType &triangulation,
                        bool bypassCache = false,
                        const std::vector<bool> *updateMask = nullptr);
 
#ifdef TTK_ENABLE_MPI
      template <typename triangulationType>
      int getSimplexRankWithGradientType(const triangulationType &triangulation,
                                         const ttk::SimplexId &lid,
                                         const int gradientType);
#endif
      inline void setInputScalarField(const void *const data,
                                      const size_t mTime) {
        inputScalarField_ = std::make_pair(data, mTime);
      }
 
      inline void setBackend(const BACKEND newBackend) {
        OldBackEnd = BackEnd;
        BackEnd = newBackend;
      }
 
      inline void setSeed(const unsigned int &newSeed) {
        OldSeed = Seed;
        Seed = newSeed;
      }
 
      inline void
        setReturnSaddleConnectors(const bool &returnSaddleConnectors) {
        OldReturnSaddleConnectors = ReturnSaddleConnectors;
        ReturnSaddleConnectors = returnSaddleConnectors;
      }
 
      inline bool newParameters() {
        return OldSeed != Seed || OldBackEnd != BackEnd
               || OldReturnSaddleConnectors != ReturnSaddleConnectors;
      }
 
      inline void setSaddleConnectorsPersistenceThreshold(double threshold) {
        SaddleConnectorsPersistenceThreshold = threshold;
      }

      inline void preconditionTriangulation(AbstractTriangulation *const data) {
        if(data != nullptr) {
          const auto dim{data->getDimensionality()};
 
          data->preconditionBoundaryVertices();
          data->preconditionVertexNeighbors();
          data->preconditionVertexEdges();
          data->preconditionVertexStars();
          data->preconditionEdges();
          data->preconditionEdgeStars();
          if(dim >= 2) {
            data->preconditionBoundaryEdges();
          }
          if(dim == 2) {
            data->preconditionCellEdges();
          } else if(dim == 3) {
            data->preconditionBoundaryTriangles();
            data->preconditionVertexTriangles();
            data->preconditionEdgeTriangles();
            data->preconditionTriangles();
            data->preconditionTriangleEdges();
            data->preconditionTriangleStars();
            data->preconditionCellTriangles();
          }
        }
      }
 
      static inline void
        clearCache(const AbstractTriangulation &triangulation) {
        triangulation.gradientCache_.clear();
      }

      inline void setLocalGradient() {
        this->gradient_ = &this->localGradient_;
      }

      inline void setInputOffsets(const SimplexId *const data) {
        inputOffsets_ = data;
      }

      int getDimensionality() const;

      int getNumberOfDimensions() const;

      template <typename triangulationType>
      SimplexId getNumberOfCells(const int dimension,
                                 const triangulationType &triangulation) const;

      template <typename triangulationType>
      bool isBoundary(const Cell &cell,
                      const triangulationType &triangulation) const;

      bool isCellCritical(const int cellDim, const SimplexId cellId) const;
      bool isCellCritical(const Cell &cell) const;

      template <typename triangulationType>
      SimplexId getPairedCell(const Cell &cell,
                              const triangulationType &triangulation,
                              bool isReverse = false) const;

      template <typename triangulationType>
      int getAscendingPath(const Cell &cell,
                           std::vector<Cell> &vpath,
                           const triangulationType &triangulation,
                           const bool enableCycleDetector = false) const;

      template <typename triangulationType>
      int getDescendingPath(const Cell &cell,
                            std::vector<Cell> &vpath,
                            const triangulationType &triangulation) const;

      template <typename triangulationType>
      bool getDescendingPathThroughWall(const Cell &saddle2,
                                        const Cell &saddle1,
                                        const std::vector<bool> &isVisited,
                                        std::vector<Cell> *const vpath,
                                        const triangulationType &triangulation,
                                        const bool stopIfMultiConnected = false,
                                        const bool enableCycleDetector
                                        = false) const;

      template <typename triangulationType>
      bool getAscendingPathThroughWall(const Cell &saddle1,
                                       const Cell &saddle2,
                                       const std::vector<bool> &isVisited,
                                       std::vector<Cell> *const vpath,
                                       const triangulationType &triangulation,
                                       const bool stopIfMultiConnected = false,
                                       const bool enableCycleDetector = false,
                                       bool *const cycleFound = nullptr) const;

      template <typename triangulationType>
      bool detectGradientCycle(const Cell &cell,
                               const triangulationType &triangulation) const;

      template <typename triangulationType>
      int getDescendingWall(const Cell &cell,
                            VisitedMask &mask,
                            const triangulationType &triangulation,
                            std::vector<Cell> *const wall = nullptr,
                            std::vector<SimplexId> *const saddles
                            = nullptr) const;

      template <typename triangulationType>
      int getAscendingWall(const Cell &cell,
                           VisitedMask &mask,
                           const triangulationType &triangulation,
                           std::vector<Cell> *const wall = nullptr,
                           std::vector<SimplexId> *const saddles
                           = nullptr) const;

      template <typename triangulationType>
      SimplexId
        getCellGreaterVertex(const Cell c,
                             const triangulationType &triangulation) const;

      template <typename triangulationType>
      SimplexId
        getCellLowerVertex(const Cell c,
                           const triangulationType &triangulation) const;

      template <typename triangulationType>
      int setCriticalPoints(
        const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
        std::vector<std::array<float, 3>> &points,
        std::vector<char> &cellDimensions,
        std::vector<SimplexId> &cellIds,
        std::vector<char> &isOnBoundary,
        std::vector<SimplexId> &PLVertexIdentifiers,
        const triangulationType &triangulation) const;

      template <typename triangulationType>
      int setCriticalPoints(std::vector<std::array<float, 3>> &points,
                            std::vector<char> &cellDimensions,
                            std::vector<SimplexId> &cellIds,
                            std::vector<char> &isOnBoundary,
                            std::vector<SimplexId> &PLVertexIdentifiers,
                            const triangulationType &triangulation) const;

      template <typename triangulationType>
      int getCriticalPoints(
        std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
        const triangulationType &triangulation) const;

      template <typename triangulationType>
      int getCriticalPoints(std::vector<Cell> &criticalPoints,
                            const triangulationType &triangulation) const;
 
#ifdef TTK_ENABLE_MPI
      void setCellToGhost(const int cellDim, const SimplexId cellId);
 
#endif
      int setManifoldSize(
        const std::array<std::vector<SimplexId>, 4> &criticalCellsByDim,
        const SimplexId *const ascendingManifold,
        const SimplexId *const descendingManifold,
        std::vector<SimplexId> &manifoldSize) const;

      template <typename triangulationType>
      int setGradientGlyphs(std::vector<std::array<float, 3>> &points,
                            std::vector<char> &points_pairOrigins,
                            std::vector<char> &cells_pairTypes,
                            std::vector<SimplexId> &cellsIds,
                            std::vector<char> &cellsDimensions,
                            const triangulationType &triangulation) const;

      CriticalType criticalTypeFromCellDimension(const int dim) const;
 
    private:
      using lowerStarType = std::array<std::vector<CellExt>, 4>;

      template <typename triangulationType>
      inline void lowerStar(lowerStarType &ls,
                            const SimplexId a,
                            const SimplexId *const offsets,
                            const triangulationType &triangulation) const;

      template <typename triangulationType>
      inline void lowerStarWithMask(lowerStarType &ls,
                                    const SimplexId a,
                                    const SimplexId *const offsets,
                                    const triangulationType &triangulation,
                                    const std::vector<bool> *updateMask
                                    = nullptr) const;

      std::pair<size_t, SimplexId>
        numUnpairedFaces(const CellExt &c, const lowerStarType &ls) const;
      std::pair<size_t, SimplexId>
        numUnpairedFacesTriangle(const CellExt &c,
                                 const lowerStarType &ls) const;
      std::pair<size_t, SimplexId>
        numUnpairedFacesTetra(const CellExt &c, const lowerStarType &ls) const;

      template <typename triangulationType>
      inline void pairCells(CellExt &alpha,
                            CellExt &beta,
                            const triangulationType &triangulation);

      template <typename triangulationType>
      int processLowerStars(const SimplexId *const offsets,
                            const triangulationType &triangulation);

      template <typename triangulationType>
      int processLowerStarsWithMask(const SimplexId *const offsets,
                                    const triangulationType &triangulation,
                                    const std::vector<bool> *updateMask
                                    = nullptr);
 
      template <typename dataType, typename triangulationType>
      int processLowerStarsStochastic(const SimplexId *const offsets,
                                      const triangulationType &triangulation);
 
      template <typename triangulationType>
      void buildStencil(const SimplexId &x,
                        const std::array<float, 3> xCoords,
                        const triangulationType &triangulation,
                        std::vector<SimplexId> &stencilIds,
                        std::vector<float> &stencilLength);

      void computeDerivatives(
        const SimplexId &x,
        const std::vector<SimplexId> &stencilIds,
        const std::array<float, 3> &xCoords,
        const std::vector<std::array<float, 3>> &stencilCoords,
        double (&grad)[3]);

      void initMemory(const AbstractTriangulation &triangulation);

 
      template <typename dataType, typename triangulationType>
      int getRemovableSaddles1(
        const std::vector<std::pair<SimplexId, char>> &criticalPoints,
        const bool allowBoundary,
        std::vector<char> &isRemovableSaddle,
        std::vector<SimplexId> &pl2dmt_saddle,
        const triangulationType &triangulation);

 
      template <typename dataType, typename triangulationType>
      int getRemovableSaddles2(
        const std::vector<std::pair<SimplexId, char>> &criticalPoints,
        const bool allowBoundary,
        std::vector<char> &isRemovableSaddle,
        std::vector<SimplexId> &pl2dmt_saddle,
        const triangulationType &triangulation);

      template <typename dataType, typename triangulationType>
      int initializeSaddleSaddleConnections1(
        const std::vector<char> &isRemovableSaddle1,
        const std::vector<char> &isRemovableSaddle2,
        const bool allowBruteForce,
        std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::vector<SimplexId> &saddle1Index,
        std::vector<SimplexId> &saddle2Index,
        const triangulationType &triangulation) const;

 
      template <typename dataType, typename triangulationType>
      int initializeSaddleSaddleConnections2(
        const std::vector<char> &isRemovableSaddle1,
        const std::vector<char> &isRemovableSaddle2,
        const bool allowBruteForce,
        std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::vector<SimplexId> &saddle1Index,
        std::vector<SimplexId> &saddle2Index,
        const triangulationType &triangulation) const;

 
      template <typename dataType>
      int orderSaddleSaddleConnections1(
        const std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::set<std::tuple<dataType, SimplexId, SimplexId>,
                 SaddleSaddleVPathComparator<dataType>> &S);

 
      template <typename dataType>
      int orderSaddleSaddleConnections2(
        const std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::set<std::tuple<dataType, SimplexId, SimplexId>,
                 SaddleSaddleVPathComparator<dataType>> &S);

      template <typename dataType, typename triangulationType>
      int processSaddleSaddleConnections1(
        const int iterationThreshold,
        const std::vector<char> &isPL,
        const bool allowBoundary,
        const bool allowBruteForce,
        const bool returnSaddleConnectors,
        std::set<std::tuple<dataType, SimplexId, SimplexId>,
                 SaddleSaddleVPathComparator<dataType>> &S,
        std::vector<SimplexId> &pl2dmt_saddle1,
        std::vector<SimplexId> &pl2dmt_saddle2,
        std::vector<char> &isRemovableSaddle1,
        std::vector<char> &isRemovableSaddle2,
        std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::vector<SimplexId> &saddle1Index,
        std::vector<SimplexId> &saddle2Index,
        const triangulationType &triangulation);

 
      template <typename dataType, typename triangulationType>
      int processSaddleSaddleConnections2(
        const int iterationThreshold,
        const std::vector<char> &isPL,
        const bool allowBoundary,
        const bool allowBruteForce,
        const bool returnSaddleConnectors,
        std::set<std::tuple<dataType, SimplexId, SimplexId>,
                 SaddleSaddleVPathComparator<dataType>> &S,
        std::vector<SimplexId> &pl2dmt_saddle1,
        std::vector<SimplexId> &pl2dmt_saddle2,
        std::vector<char> &isRemovableSaddle1,
        std::vector<char> &isRemovableSaddle2,
        std::vector<VPath> &vpaths,
        std::vector<CriticalPoint> &criticalPoints,
        std::vector<SimplexId> &saddle1Index,
        std::vector<SimplexId> &saddle2Index,
        const triangulationType &triangulation);
 
    public:
      template <typename dataType, typename triangulationType>
      dataType getPersistence(const Cell &up,
                              const Cell &down,
                              const dataType *const scalars,
                              const triangulationType &triangulation) const;

      bool isMinimum(const Cell &cell) const;

      bool isSaddle1(const Cell &cell) const;

      bool isSaddle2(const Cell &cell) const;

      bool isMaximum(const Cell &cell) const;

      template <typename triangulationType>
      int reverseAscendingPath(const std::vector<Cell> &vpath,
                               const triangulationType &triangulation) const;

      template <typename triangulationType>
      int reverseDescendingPath(const std::vector<Cell> &vpath,
                                const triangulationType &triangulation) const;

      template <typename triangulationType>
      int reverseAscendingPathOnWall(const std::vector<Cell> &vpath,
                                     const triangulationType &triangulation,
                                     bool cancelReversal = false) const;

      template <typename triangulationType>
      int reverseDescendingPathOnWall(
        const std::vector<Cell> &vpath,
        const triangulationType &triangulation) const;

      int getCriticalPointMap(
        const std::vector<std::pair<SimplexId, char>> &criticalPoints,
        std::vector<char> &isPL);

      template <typename dataType, typename triangulationType>
      int simplifySaddleSaddleConnections1(
        const std::vector<std::pair<SimplexId, char>> &criticalPoints,
        const std::vector<char> &isPL,
        const int iterationThreshold,
        const bool allowBoundary,
        const bool allowBruteForce,
        const bool returnSaddleConnectors,
        const triangulationType &triangulation);

      template <typename dataType, typename triangulationType>
      int simplifySaddleSaddleConnections2(
        const std::vector<std::pair<SimplexId, char>> &criticalPoints,
        const std::vector<char> &isPL,
        const int iterationThreshold,
        const bool allowBoundary,
        const bool allowBruteForce,
        const bool returnSaddleConnectors,
        const triangulationType &triangulation);
 
      template <typename dataType, typename triangulationType>
      int filterSaddleConnectors(const bool allowBoundary,
                                 const triangulationType &triangulation);
 
    protected:
      int dimensionality_{-1};
      int IterationThreshold{50};
      double SaddleConnectorsPersistenceThreshold{};
      SimplexId numberOfVertices_{};
      BACKEND BackEnd{BACKEND::CLASSIC_BACKEND};
      BACKEND OldBackEnd{BACKEND::CLASSIC_BACKEND};
      unsigned int Seed{};
      unsigned int OldSeed{};
      bool ReturnSaddleConnectors{};
      bool OldReturnSaddleConnectors{};
      // spare storage (bypass cache) for gradient internal structure
      AbstractTriangulation::gradientType localGradient_{};
      // cache key (scalar field pointer + timestamp)
      AbstractTriangulation::gradientKeyType inputScalarField_{};
      // pointer to either cache entry corresponding to inputScalarField_ or
      // localGradient_ (if cache is bypassed)
      AbstractTriangulation::gradientType *gradient_{};
      const SimplexId *inputOffsets_{};
    };
 
  } // namespace dcg
} // namespace ttk
 
#include <DiscreteGradient_Template.h>