DiscreteVectorField.h

Go to the documentation of this file.

 
#pragma once
 
// base code includes
#include <DiscreteGradient.h>
#include <Geometry.h>
#include <SurfaceGeometrySmoother.h>
#include <Triangulation.h>
#include <VisitedMask.h>
 
#include <algorithm>
#include <array>
#include <functional>
#include <queue>
#include <set>
#include <utility>
 
using ttk::dcg::Cell;
 
namespace ttk {
 
  namespace dcvf {
 
    enum connectionValue { NULL_CONNECTION = -1, GHOST_CONNECTION = -2 };

    struct CellOutExt : Cell {
      explicit CellOutExt(const int dim, const SimplexId id) : Cell{dim, id} {
      }
      explicit CellOutExt(const int dim,
                          const SimplexId id,
                          const std::array<SimplexId, 3> &lowVerts,
                          const std::array<float, 3> &lowVertWeights,
                          const std::array<uint8_t, 3> &faces)
        : Cell{dim, id}, lowVerts_{lowVerts},
          lowVertWeights_{lowVertWeights}, faces_{faces} {
      }
 
      // ID values for Outward vertices in current Outward star
      // (1 for edges, 2 for triangles, 3 for tetras)
      const std::array<SimplexId, 3> lowVerts_{};
      // (weightValue for) Outward vertices in current Outward star, align with
      // lowVerts_ (1 for edges, 2 for triangles, 3 for tetras)
      const std::array<float, 3> lowVertWeights_{};
      // indices of faces (cells of dimensions dim_ - 1) in Outward 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 Outward star
      bool paired_{false};
    };

    using vectorValue = SurfaceGeometrySmoother::Point;

    class DiscreteVectorField : virtual public Debug {
 
    public:
      DiscreteVectorField() {
        this->setDebugMsgPrefix("DiscreteVectorField");
#ifdef TTK_ENABLE_MPI
        hasMPISupport_ = true;
#endif
      }

      void setReverseFullOrbit(bool data) {
        reverseFullOrbit = data;
      }

      template <typename dataType, typename triangulationType>
      int buildField(const triangulationType &triangulation);

      template <typename dataType, typename triangulationType>
      bool compare(const triangulationType &triangulation,
                   SimplexId vertexA,
                   SimplexId vertexB,
                   float &weightValue) const;

      inline void setInputVectorField(const void *const data,
                                      const size_t mTime) {
        inputVectorField_ = std::make_pair(data, mTime);
      }

      template <typename dataType>
      vectorValue getVectorValueAt(SimplexId &vertex) const {
#ifndef TTK_ENABLE_KAMIKAZE
        if(inputVectorField_.first == nullptr) {
          this->printErr(
            "Attempting to access vector before setting input vector field");
          return {};
        }
#endif
        vectorValue result;
        const dataType *dataArray
          = reinterpret_cast<const dataType *>(inputVectorField_.first);
        result[0] = static_cast<float>(dataArray[3 * vertex + 0]);
        result[1] = static_cast<float>(dataArray[3 * vertex + 1]);
        result[2] = static_cast<float>(dataArray[3 * vertex + 2]);
 
        return result;
      }

      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();
 
          data->preconditionEdgeTriangles();
          data->preconditionTriangles();
          if(dim >= 2) {
            data->preconditionBoundaryEdges();
          }
          if(dim == 2) {
            data->preconditionCellEdges();
          } else if(dim == 3) {
            data->preconditionBoundaryTriangles();
            data->preconditionVertexTriangles();
            data->preconditionTriangleEdges();
            data->preconditionTriangleStars();
            data->preconditionCellTriangles();
          }
        }
      }

      int getDimensionality() const;

      int getNumberOfDimensions() const;

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

      template <typename dataType, 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 dataType, typename triangulationType>
      int getAscendingPath(const Cell &cell,
                           std::vector<Cell> &vpath,
                           const triangulationType &triangulation,
                           const bool stopOnCycle) const;

      template <typename dataType, typename triangulationType>
      int getAscendingPathRecursive(const Cell &cell,
                                    std::vector<Cell> &vpath,
                                    const triangulationType &triangulation,
                                    std::vector<char> &previousDescPaths,
                                    std::vector<char> &previousAscPaths) const;

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

      template <typename dataType, typename triangulationType>
      int getDescendingPathRecursive(const Cell &cell,
                                     std::vector<Cell> &vpath,
                                     const triangulationType &triangulation,
                                     std::vector<char> &previousDescPaths,
                                     std::vector<char> &previousAscPaths) 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>
      void getAscendingPathThroughWall(
        const Cell &saddle1,
        const Cell &saddle2,
        const std::vector<bool> &isVisited,
        std::vector<Cell> *const vpath,
        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 dataType, typename triangulationType>
      SimplexId
        getCellGreaterVertex(const Cell c,
                             const triangulationType &triangulation) const;

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

      template <typename dataType, 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 dataType, 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 setVectorGlyphs(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 outwardStarType = std::array<std::vector<CellOutExt>, 4>;

      template <typename dataType, typename triangulationType>
      inline void OutwardStar(outwardStarType &os,
                              const SimplexId a,
                              const triangulationType &triangulation) const;

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

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

      template <typename dataType, typename triangulationType>
      int processOutwardStars(const triangulationType &triangulation);

      void initMemory(const AbstractTriangulation &triangulation);
 
    public:
      template <typename dataType, typename triangulationType>
      float getPersistence(const std::vector<Cell> &vpath,
                           const triangulationType &triangulation) 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 reverseAlternatingPath(const std::vector<Cell> &vpath,
                                 const triangulationType &triangulation) const;

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

      template <typename triangulationType>
      int reverseDescendingPathOnWall(
        const std::vector<Cell> &vpath,
        const triangulationType &triangulation) const;
 
    protected:
      int dimensionality_{-1};
      SimplexId numberOfVertices_{};
      bool reverseFullOrbit{true};
 
      // spare storage for discrete vectors internal structure
      std::array<std::vector<SimplexId>, 6> localVectors_{};
      // former cache key (vector field pointer + timestamp)
      std::pair<const void *, size_t> inputVectorField_{};
      // pointer to either cache entry corresponding to inputVectorField_ or
      // localVectors_ (because cache is not implemented)
      std::array<std::vector<SimplexId>, 6> *vectors_{};
    };
 
  } // namespace dcvf
} // namespace ttk
 
#include <DiscreteVectorField_Template.h>