GeodesicsInHeat.h

 
#pragma once
 
#if defined(GeodesicsInHeat_RECURSES)
#error Recursive header files inclusion detected in GeodesicsInHeat.h
#else // defined(GeodesicsInHeat_RECURSES)
#define GeodesicsInHeat_RECURSES
 
#if !defined GeodesicsInHeat_h
#define GeodesicsInHeat_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/base/ConstAlias.h"
#include "DGtal/math/linalg/DirichletConditions.h"
 
namespace DGtal
{
  // template class GeodesicsInHeat
  template <typename TPolygonalCalculus>
  class GeodesicsInHeat
  {
    // ----------------------- Standard services ------------------------------
  public:
 
    typedef TPolygonalCalculus PolygonalCalculus;
    typedef typename PolygonalCalculus::SparseMatrix SparseMatrix;
    typedef typename PolygonalCalculus::DenseMatrix DenseMatrix;
    typedef typename PolygonalCalculus::Solver Solver;
    typedef typename PolygonalCalculus::Vector Vector;
    typedef typename PolygonalCalculus::Vertex Vertex;
    typedef typename PolygonalCalculus::LinAlg LinAlgBackend;
    typedef DirichletConditions< LinAlgBackend > Conditions;
    typedef typename Conditions::IntegerVector IntegerVector;
    
    GeodesicsInHeat() = delete;
    
    GeodesicsInHeat(ConstAlias<PolygonalCalculus> calculus): myCalculus(&calculus)
    {
      myIsInit=false;
    }
    
    ~GeodesicsInHeat() = default;
 
    GeodesicsInHeat ( const GeodesicsInHeat & other ) = delete;
 
    GeodesicsInHeat ( GeodesicsInHeat && other ) = delete;
 
    GeodesicsInHeat & operator= ( const GeodesicsInHeat & other ) = delete;
 
    GeodesicsInHeat & operator= ( GeodesicsInHeat && other ) = delete;
 
 
    
    // ----------------------- Interface --------------------------------------
   
    void init( double dt, double lambda = 1.0,
               bool boundary_with_mixed_solution = false  )
    {
      myIsInit = true;
      myLambda = lambda;
 
      SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
      SparseMatrix mass      = myCalculus->globalLumpedMassMatrix();
      myHeatOpe              = mass - dt*laplacian;
    
      // from https://geometry-central.net
      // NOTE: In theory, it should not be necessary to shift the Laplacian: the Polydec Laplace is always PSD. However, when the
      // matrix is only positive SEMIdefinite, some solvers may not work (ie Eigen's Cholesky solver doesn't work, but
      // Suitesparse does).
      SparseMatrix Id = SparseMatrix(myCalculus->nbVertices(),myCalculus->nbVertices());
      Id.setIdentity();
      laplacian += 1e-6 * Id;
      
      //Prefactorizing
      myPoissonSolver.compute( laplacian );
      myHeatSolver.compute   ( myHeatOpe );
      
      //empty source
      mySource    = Vector::Zero(myCalculus->nbVertices());
 
      // Manage boundaries
      myManageBoundary = false;
      if ( ! boundary_with_mixed_solution ) return;
      myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
      const auto surfmesh = myCalculus->getSurfaceMeshPtr();
      const auto edges    = surfmesh->computeManifoldBoundaryEdges();
      for ( auto e : edges )
        {
          const auto vtcs = surfmesh->edgeVertices( e );
          myBoundary[ vtcs.first  ] = 1;
          myBoundary[ vtcs.second ] = 1;
        }
      myManageBoundary = ! edges.empty();
      if ( ! myManageBoundary ) return;
      // Prepare solver for a problem with Dirichlet conditions.
      SparseMatrix heatOpe_d = Conditions::dirichletOperator( myHeatOpe, myBoundary );
      // Prefactoring
      myHeatDirichletSolver.compute( heatOpe_d );
    }
    
    void addSource(const Vertex aV)
    {
      ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
      myLastSourceIndex = aV;
      mySource( aV ) = 1.0;
    }
    
    void clearSource()
    {
      mySource = Vector::Zero(myCalculus->nbVertices());
    }
    
    Vector source() const
    {
      FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
      return mySource;
    }
    
    
    Vector compute() const
    {
      FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
      //Heat diffusion
      Vector heatDiffusion = myHeatSolver.solve(mySource);
      ASSERT(myHeatSolver.info()==Eigen::Success);
 
      // Take care of boundaries
      if ( myManageBoundary )
        {
          Vector bValues  = Vector::Zero( myCalculus->nbVertices() );
          Vector bSources = Conditions::dirichletVector( myHeatOpe, mySource,
                                                         myBoundary, bValues );
          Vector bSol     = myHeatDirichletSolver.solve( bSources );
          Vector heatDiffusionDirichlet
                          = Conditions::dirichletSolution( bSol, myBoundary, bValues );
          heatDiffusion = 0.5 * ( heatDiffusion + heatDiffusionDirichlet );
        }
      Vector divergence    = Vector::Zero(myCalculus->nbVertices());
      auto cpt=0;
      auto surfmesh = myCalculus->getSurfaceMeshPtr();
      
      // Heat, normalization and divergence per face
      for(typename PolygonalCalculus::MySurfaceMesh::Index f=0; f< myCalculus->nbFaces(); ++f)
        {
          Vector faceHeat( myCalculus->degree(f));
          cpt=0;
          auto vertices = surfmesh->incidentVertices(f);
          for(auto v: vertices)
            {
              faceHeat(cpt) = heatDiffusion( v );
              ++cpt;
            }
          // ∇heat / ∣∣∇heat∣∣
          Vector grad = -myCalculus->gradient(f) * faceHeat;
          grad.normalize();
      
          // div
          DenseMatrix   oneForm = myCalculus->flat(f)*grad;
          Vector divergenceFace = myCalculus->divergence( f ) * oneForm;
          cpt=0;
          for(auto v: vertices)
            {
              divergence(v) += divergenceFace(cpt);
              ++cpt;
            }
        }
      
      // Last Poisson solve
      Vector distVec = myPoissonSolver.solve(divergence);
      ASSERT(myPoissonSolver.info()==Eigen::Success);
 
      //Source val
      auto sourceval = distVec(myLastSourceIndex);
      //shifting the distances to get 0 at sources
      return distVec - sourceval*Vector::Ones(myCalculus->nbVertices());
    }
    
    
    bool isValid() const
    {
      return myIsInit && myCalculus->isValid();
    }
    
    // ----------------------- Private --------------------------------------
 
  private:
    
    const PolygonalCalculus *myCalculus;
 
    SparseMatrix myHeatOpe;
    
    Solver myPoissonSolver;
 
    Solver myHeatSolver;
 
    Vector mySource;
 
    Vertex myLastSourceIndex;
  
    bool myIsInit;
 
    double myLambda;
 
    bool myManageBoundary;
 
    IntegerVector myBoundary;
    
    Solver myHeatDirichletSolver;
  
  
  }; // end of class GeodesicsInHeat
} // namespace DGtal
 
//                                                                           //
 
#endif // !defined GeodesicsInHeat_h
 
#undef GeodesicsInHeat_RECURSES
#endif // else defined(GeodesicsInHeat_RECURSES)