This class implements [41] on polygonal surfaces (using Discrete differential calculus on polygonal surfaces). More...

#include <DGtal/dec/GeodesicsInHeat.h>

Public Types
typedef TPolygonalCalculus	PolygonalCalculus
typedef PolygonalCalculus::SparseMatrix	SparseMatrix
typedef PolygonalCalculus::DenseMatrix	DenseMatrix
typedef PolygonalCalculus::Solver	Solver
typedef PolygonalCalculus::Vector	Vector
typedef PolygonalCalculus::Vertex	Vertex
typedef PolygonalCalculus::LinAlg	LinAlgBackend
typedef DirichletConditions< LinAlgBackend >	Conditions
typedef Conditions::IntegerVector	IntegerVector

Public Member Functions
	GeodesicsInHeat ()=delete
	GeodesicsInHeat (ConstAlias< PolygonalCalculus > calculus)
	~GeodesicsInHeat ()=default
	GeodesicsInHeat (const GeodesicsInHeat &other)=delete
	GeodesicsInHeat (GeodesicsInHeat &&other)=delete
GeodesicsInHeat &	operator= (const GeodesicsInHeat &other)=delete
GeodesicsInHeat &	operator= (GeodesicsInHeat &&other)=delete
void	init (double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
void	addSource (const Vertex aV)
void	clearSource ()
Vector	source () const
Vector	compute () const
bool	isValid () const

Private Attributes
const PolygonalCalculus *	myCalculus
	The underlying PolygonalCalculus instance.
SparseMatrix	myHeatOpe
	The operator for heat diffusion.
Solver	myPoissonSolver
	Poisson solver.
Solver	myHeatSolver
	Heat solver.
Vector	mySource
	Source vector.
Vertex	myLastSourceIndex
	Vertex index to the last source point (to shift the distances)
bool	myIsInit
	Validitate flag.
double	myLambda
	Lambda parameter.
bool	myManageBoundary
IntegerVector	myBoundary
	The boundary characteristic vector.
Solver	myHeatDirichletSolver
	Heat solver with Dirichlet boundary conditions.

Detailed Description

template<typename TPolygonalCalculus>
class DGtal::GeodesicsInHeat< TPolygonalCalculus >

This class implements [41] on polygonal surfaces (using Discrete differential calculus on polygonal surfaces).

Description of template class 'GeodesicsInHeat'

see Geodesics In Heat using discrete polygonal calculus for details and examples.

Template Parameters

a	model of PolygonalCalculus.

Definition at line 62 of file GeodesicsInHeat.h.

Member Typedef Documentation

◆ Conditions

template<typename TPolygonalCalculus>

typedef DirichletConditions< LinAlgBackend > DGtal::GeodesicsInHeat< TPolygonalCalculus >::Conditions

Definition at line 74 of file GeodesicsInHeat.h.

◆ DenseMatrix

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::DenseMatrix DGtal::GeodesicsInHeat< TPolygonalCalculus >::DenseMatrix

Definition at line 69 of file GeodesicsInHeat.h.

◆ IntegerVector

template<typename TPolygonalCalculus>

typedef Conditions::IntegerVector DGtal::GeodesicsInHeat< TPolygonalCalculus >::IntegerVector

Definition at line 75 of file GeodesicsInHeat.h.

◆ LinAlgBackend

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::LinAlg DGtal::GeodesicsInHeat< TPolygonalCalculus >::LinAlgBackend

Definition at line 73 of file GeodesicsInHeat.h.

◆ PolygonalCalculus

template<typename TPolygonalCalculus>

typedef TPolygonalCalculus DGtal::GeodesicsInHeat< TPolygonalCalculus >::PolygonalCalculus

Definition at line 67 of file GeodesicsInHeat.h.

◆ Solver

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::Solver DGtal::GeodesicsInHeat< TPolygonalCalculus >::Solver

Definition at line 70 of file GeodesicsInHeat.h.

◆ SparseMatrix

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::SparseMatrix DGtal::GeodesicsInHeat< TPolygonalCalculus >::SparseMatrix

Definition at line 68 of file GeodesicsInHeat.h.

◆ Vector

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::Vector DGtal::GeodesicsInHeat< TPolygonalCalculus >::Vector

Definition at line 71 of file GeodesicsInHeat.h.

◆ Vertex

template<typename TPolygonalCalculus>

typedef PolygonalCalculus::Vertex DGtal::GeodesicsInHeat< TPolygonalCalculus >::Vertex

Definition at line 72 of file GeodesicsInHeat.h.

Constructor & Destructor Documentation

◆ GeodesicsInHeat() [1/4]

template<typename TPolygonalCalculus>

DGtal::GeodesicsInHeat< TPolygonalCalculus >::GeodesicsInHeat ( )

delete

Default constructor.

Referenced by GeodesicsInHeat(), GeodesicsInHeat(), operator=(), and operator=().

◆ GeodesicsInHeat() [2/4]

template<typename TPolygonalCalculus>

DGtal::GeodesicsInHeat< TPolygonalCalculus >::GeodesicsInHeat ( ConstAlias< PolygonalCalculus > calculus )

inline

Constructor from an existing polygonal calculus. T

Parameters

calculus a instance of PolygonalCalculus

Definition at line 84 of file GeodesicsInHeat.h.

                                                           : myCalculus(&calculus)
    {
      myIsInit=false;
    }

References calculus, myCalculus, and myIsInit.

◆ ~GeodesicsInHeat()

template<typename TPolygonalCalculus>

DGtal::GeodesicsInHeat< TPolygonalCalculus >::~GeodesicsInHeat ( )

default

Destructor.

◆ GeodesicsInHeat() [3/4]

template<typename TPolygonalCalculus>

DGtal::GeodesicsInHeat< TPolygonalCalculus >::GeodesicsInHeat ( const GeodesicsInHeat< TPolygonalCalculus > & other )

delete

Copy constructor.

Parameters

other the object to clone.

References GeodesicsInHeat().

◆ GeodesicsInHeat() [4/4]

template<typename TPolygonalCalculus>

DGtal::GeodesicsInHeat< TPolygonalCalculus >::GeodesicsInHeat ( GeodesicsInHeat< TPolygonalCalculus > && other )

delete

Move constructor.

Parameters

other the object to move.

References GeodesicsInHeat().

Member Function Documentation

◆ addSource()

template<typename TPolygonalCalculus>

void DGtal::GeodesicsInHeat< TPolygonalCalculus >::addSource ( const Vertex aV )

inline

Adds a source point at a vertex aV

Parameters

aV	the Vertex

Definition at line 183 of file GeodesicsInHeat.h.

    {
      ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
      myLastSourceIndex = aV;
      mySource( aV ) = 1.0;
    }

References myLastSourceIndex, and mySource.

◆ clearSource()

template<typename TPolygonalCalculus>

void DGtal::GeodesicsInHeat< TPolygonalCalculus >::clearSource ( )

inline

Removes all source Diracs.

Definition at line 192 of file GeodesicsInHeat.h.

    {
      mySource = Vector::Zero(myCalculus->nbVertices());
    }

References myCalculus, and mySource.

◆ compute()

template<typename TPolygonalCalculus>

Vector DGtal::GeodesicsInHeat< TPolygonalCalculus >::compute ( ) const

inline

Main computation of the Geodesic In Heat

Returns: the estimated geodesic distances from the sources.

Definition at line 209 of file GeodesicsInHeat.h.

    {
      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());
    }

References DGtal::DirichletConditions< LinAlgBackend >::dirichletSolution(), DGtal::DirichletConditions< LinAlgBackend >::dirichletVector(), myBoundary, myCalculus, myHeatDirichletSolver, myHeatOpe, myHeatSolver, myIsInit, myLastSourceIndex, myManageBoundary, myPoissonSolver, mySource, and surfmesh.

◆ init()

template<typename TPolygonalCalculus>

void DGtal::GeodesicsInHeat< TPolygonalCalculus >::init	(	double	dt,
		double	lambda = 1.0,
		bool	boundary_with_mixed_solution = false )

inline

Initialize the solvers with dt as timestep for the heat diffusion and lambda parameter for the polygonal calculus, which guarantee definiteness for positive lambda.

Parameters

dt	the timestep
lambda	timestep
boundary_with_mixed_solution	when 'true' and when the surface has boundaries, mix two solutions of the heat diffusion operation (Neumann and Dirichlet null conditions on boundary).

Definition at line 135 of file GeodesicsInHeat.h.

    {
      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 );
    }

References DGtal::DirichletConditions< LinAlgBackend >::dirichletOperator(), dt, laplacian(), myBoundary, myCalculus, myHeatDirichletSolver, myHeatOpe, myHeatSolver, myIsInit, myLambda, myManageBoundary, myPoissonSolver, mySource, and surfmesh.

◆ isValid()

template<typename TPolygonalCalculus>

bool DGtal::GeodesicsInHeat< TPolygonalCalculus >::isValid ( ) const

inline

Returns: true if the calculus is valid.

Definition at line 269 of file GeodesicsInHeat.h.

    {
      return myIsInit && myCalculus->isValid();
    }

References myCalculus, and myIsInit.

◆ operator=() [1/2]

template<typename TPolygonalCalculus>

GeodesicsInHeat & DGtal::GeodesicsInHeat< TPolygonalCalculus >::operator= ( const GeodesicsInHeat< TPolygonalCalculus > & other )

delete

Copy assignment operator.

Parameters

other the object to copy.

Returns: a reference on 'this'.

References GeodesicsInHeat().

◆ operator=() [2/2]

template<typename TPolygonalCalculus>

GeodesicsInHeat & DGtal::GeodesicsInHeat< TPolygonalCalculus >::operator= ( GeodesicsInHeat< TPolygonalCalculus > && other )

delete

Move assignment operator.

Parameters

other the object to move.

Returns: a reference on 'this'.

References GeodesicsInHeat().

◆ source()

template<typename TPolygonalCalculus>

Vector DGtal::GeodesicsInHeat< TPolygonalCalculus >::source ( ) const

inline

Returns: the source point vector.

Definition at line 200 of file GeodesicsInHeat.h.

    {
      FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
      return mySource;
    }

References myIsInit, and mySource.

Field Documentation

◆ myBoundary

template<typename TPolygonalCalculus>

IntegerVector DGtal::GeodesicsInHeat< TPolygonalCalculus >::myBoundary

private

The boundary characteristic vector.

Definition at line 307 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ myCalculus

template<typename TPolygonalCalculus>

const PolygonalCalculus* DGtal::GeodesicsInHeat< TPolygonalCalculus >::myCalculus

private

The underlying PolygonalCalculus instance.

Definition at line 279 of file GeodesicsInHeat.h.

Referenced by clearSource(), compute(), GeodesicsInHeat(), init(), and isValid().

◆ myHeatDirichletSolver

template<typename TPolygonalCalculus>

Solver DGtal::GeodesicsInHeat< TPolygonalCalculus >::myHeatDirichletSolver

private

Heat solver with Dirichlet boundary conditions.

Definition at line 310 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ myHeatOpe

template<typename TPolygonalCalculus>

SparseMatrix DGtal::GeodesicsInHeat< TPolygonalCalculus >::myHeatOpe

private

The operator for heat diffusion.

Definition at line 282 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ myHeatSolver

template<typename TPolygonalCalculus>

Solver DGtal::GeodesicsInHeat< TPolygonalCalculus >::myHeatSolver

private

Heat solver.

Definition at line 288 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ myIsInit

template<typename TPolygonalCalculus>

bool DGtal::GeodesicsInHeat< TPolygonalCalculus >::myIsInit

private

Validitate flag.

Definition at line 297 of file GeodesicsInHeat.h.

Referenced by compute(), GeodesicsInHeat(), init(), isValid(), and source().

◆ myLambda

template<typename TPolygonalCalculus>

double DGtal::GeodesicsInHeat< TPolygonalCalculus >::myLambda

private

Lambda parameter.

Definition at line 300 of file GeodesicsInHeat.h.

Referenced by init().

◆ myLastSourceIndex

template<typename TPolygonalCalculus>

Vertex DGtal::GeodesicsInHeat< TPolygonalCalculus >::myLastSourceIndex

private

Vertex index to the last source point (to shift the distances)

Definition at line 294 of file GeodesicsInHeat.h.

Referenced by addSource(), and compute().

◆ myManageBoundary

template<typename TPolygonalCalculus>

bool DGtal::GeodesicsInHeat< TPolygonalCalculus >::myManageBoundary

private

When 'true', manage boundaries with a mixed solution of Neumann and Dirichlet conditions.

Definition at line 304 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ myPoissonSolver

template<typename TPolygonalCalculus>

Solver DGtal::GeodesicsInHeat< TPolygonalCalculus >::myPoissonSolver

private

Poisson solver.

Definition at line 285 of file GeodesicsInHeat.h.

Referenced by compute(), and init().

◆ mySource

template<typename TPolygonalCalculus>

Vector DGtal::GeodesicsInHeat< TPolygonalCalculus >::mySource

private

Source vector.

Definition at line 291 of file GeodesicsInHeat.h.

Referenced by addSource(), clearSource(), compute(), init(), and source().

The documentation for this class was generated from the following file:

GeodesicsInHeat.h

Public Types

Public Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ Conditions

◆ DenseMatrix

◆ IntegerVector

◆ LinAlgBackend

◆ PolygonalCalculus

◆ Solver

◆ SparseMatrix

◆ Vector

◆ Vertex

Constructor & Destructor Documentation

◆ GeodesicsInHeat() [1/4]

◆ GeodesicsInHeat() [2/4]

◆ ~GeodesicsInHeat()

◆ GeodesicsInHeat() [3/4]

◆ GeodesicsInHeat() [4/4]

Member Function Documentation

◆ addSource()

◆ clearSource()

◆ compute()

◆ init()

◆ isValid()

◆ operator=() [1/2]

◆ operator=() [2/2]

◆ source()

Field Documentation

◆ myBoundary

◆ myCalculus

◆ myHeatDirichletSolver

◆ myHeatOpe

◆ myHeatSolver

◆ myIsInit

◆ myLambda

◆ myLastSourceIndex

◆ myManageBoundary

◆ myPoissonSolver

◆ mySource