DiscreteGradient.h

Go to the documentation of this file.

 
 
#pragma once
 
// base code includes
#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};
    };
 
    class DiscreteGradient : virtual public Debug {
 
    public:
      DiscreteGradient() {
        this->setDebugMsgPrefix("DiscreteGradient");
#ifdef TTK_ENABLE_MPI
        hasMPISupport_ = true;
#endif
      }
 
      template <typename triangulationType>
      int buildGradient(const triangulationType &triangulation,
                        bool bypassCache = false);
 
      inline void setInputScalarField(const void *const data,
                                      const size_t mTime) {
        inputScalarField_ = std::make_pair(data, mTime);
      }
 
      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;
 
#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;
 
    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;
 
      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;
 
      CriticalType criticalTypeFromCellDimension(const int dim) 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);
 
      void initMemory(const AbstractTriangulation &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;
 
    protected:
      int dimensionality_{-1};
      SimplexId numberOfVertices_{};
 
      // 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>