DGtal  1.5.beta
GeodesicsInHeat.h
1 
17 #pragma once
18 
31 #if defined(GeodesicsInHeat_RECURSES)
32 #error Recursive header files inclusion detected in GeodesicsInHeat.h
33 #else // defined(GeodesicsInHeat_RECURSES)
35 #define GeodesicsInHeat_RECURSES
36 
37 #if !defined GeodesicsInHeat_h
39 #define GeodesicsInHeat_h
40 
42 // Inclusions
43 #include <iostream>
44 #include "DGtal/base/Common.h"
45 #include "DGtal/base/ConstAlias.h"
46 #include "DGtal/math/linalg/DirichletConditions.h"
48 
49 namespace DGtal
50 {
52  // template class GeodesicsInHeat
61  template <typename TPolygonalCalculus>
63  {
64  // ----------------------- Standard services ------------------------------
65  public:
66 
67  typedef TPolygonalCalculus PolygonalCalculus;
76 
80  GeodesicsInHeat() = delete;
81 
85  {
86  myIsInit=false;
87  }
88 
92  ~GeodesicsInHeat() = default;
93 
98  GeodesicsInHeat ( const GeodesicsInHeat & other ) = delete;
99 
104  GeodesicsInHeat ( GeodesicsInHeat && other ) = delete;
105 
111  GeodesicsInHeat & operator= ( const GeodesicsInHeat & other ) = delete;
112 
119 
120 
121 
122  // ----------------------- Interface --------------------------------------
123 
135  void init( double dt, double lambda = 1.0,
136  bool boundary_with_mixed_solution = false )
137  {
138  myIsInit = true;
139  myLambda = lambda;
140 
141  SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
142  SparseMatrix mass = myCalculus->globalLumpedMassMatrix();
143  myHeatOpe = mass - dt*laplacian;
144 
145  // from https://geometry-central.net
146  // NOTE: In theory, it should not be necessary to shift the Laplacian: the Polydec Laplace is always PSD. However, when the
147  // matrix is only positive SEMIdefinite, some solvers may not work (ie Eigen's Cholesky solver doesn't work, but
148  // Suitesparse does).
149  SparseMatrix Id = SparseMatrix(myCalculus->nbVertices(),myCalculus->nbVertices());
150  Id.setIdentity();
151  laplacian += 1e-6 * Id;
152 
153  //Prefactorizing
154  myPoissonSolver.compute( laplacian );
155  myHeatSolver.compute ( myHeatOpe );
156 
157  //empty source
158  mySource = Vector::Zero(myCalculus->nbVertices());
159 
160  // Manage boundaries
161  myManageBoundary = false;
162  if ( ! boundary_with_mixed_solution ) return;
163  myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
164  const auto surfmesh = myCalculus->getSurfaceMeshPtr();
165  const auto edges = surfmesh->computeManifoldBoundaryEdges();
166  for ( auto e : edges )
167  {
168  const auto vtcs = surfmesh->edgeVertices( e );
169  myBoundary[ vtcs.first ] = 1;
170  myBoundary[ vtcs.second ] = 1;
171  }
172  myManageBoundary = ! edges.empty();
173  if ( ! myManageBoundary ) return;
174  // Prepare solver for a problem with Dirichlet conditions.
176  // Prefactoring
177  myHeatDirichletSolver.compute( heatOpe_d );
178  }
179 
183  void addSource(const Vertex aV)
184  {
185  ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
186  myLastSourceIndex = aV;
187  mySource( aV ) = 1.0;
188  }
189 
192  void clearSource()
193  {
194  mySource = Vector::Zero(myCalculus->nbVertices());
195  }
196 
200  Vector source() const
201  {
202  FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
203  return mySource;
204  }
205 
206 
209  Vector compute() const
210  {
211  FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
212  //Heat diffusion
213  Vector heatDiffusion = myHeatSolver.solve(mySource);
214  ASSERT(myHeatSolver.info()==Eigen::Success);
215 
216  // Take care of boundaries
217  if ( myManageBoundary )
218  {
219  Vector bValues = Vector::Zero( myCalculus->nbVertices() );
221  myBoundary, bValues );
222  Vector bSol = myHeatDirichletSolver.solve( bSources );
223  Vector heatDiffusionDirichlet
224  = Conditions::dirichletSolution( bSol, myBoundary, bValues );
225  heatDiffusion = 0.5 * ( heatDiffusion + heatDiffusionDirichlet );
226  }
227  Vector divergence = Vector::Zero(myCalculus->nbVertices());
228  auto cpt=0;
229  auto surfmesh = myCalculus->getSurfaceMeshPtr();
230 
231  // Heat, normalization and divergence per face
232  for(typename PolygonalCalculus::MySurfaceMesh::Index f=0; f< myCalculus->nbFaces(); ++f)
233  {
234  Vector faceHeat( myCalculus->degree(f));
235  cpt=0;
236  auto vertices = surfmesh->incidentVertices(f);
237  for(auto v: vertices)
238  {
239  faceHeat(cpt) = heatDiffusion( v );
240  ++cpt;
241  }
242  // ∇heat / ∣∣∇heat∣∣
243  Vector grad = -myCalculus->gradient(f) * faceHeat;
244  grad.normalize();
245 
246  // div
247  DenseMatrix oneForm = myCalculus->flat(f)*grad;
248  Vector divergenceFace = myCalculus->divergence( f ) * oneForm;
249  cpt=0;
250  for(auto v: vertices)
251  {
252  divergence(v) += divergenceFace(cpt);
253  ++cpt;
254  }
255  }
256 
257  // Last Poisson solve
258  Vector distVec = myPoissonSolver.solve(divergence);
259  ASSERT(myPoissonSolver.info()==Eigen::Success);
260 
261  //Source val
262  auto sourceval = distVec(myLastSourceIndex);
263  //shifting the distances to get 0 at sources
264  return distVec - sourceval*Vector::Ones(myCalculus->nbVertices());
265  }
266 
267 
269  bool isValid() const
270  {
271  return myIsInit && myCalculus->isValid();
272  }
273 
274  // ----------------------- Private --------------------------------------
275 
276  private:
277 
280 
283 
286 
289 
292 
295 
297  bool myIsInit;
298 
300  double myLambda;
301 
305 
308 
311 
312 
313  }; // end of class GeodesicsInHeat
314 } // namespace DGtal
315 
316 // //
318 
319 #endif // !defined GeodesicsInHeat_h
320 
321 #undef GeodesicsInHeat_RECURSES
322 #endif // else defined(GeodesicsInHeat_RECURSES)
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
Aim: A helper class to solve a system with Dirichlet boundary conditions.
LinearAlgebraBackend::IntegerVector IntegerVector
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
GeodesicsInHeat(GeodesicsInHeat &&other)=delete
DirichletConditions< LinAlgBackend > Conditions
PolygonalCalculus::LinAlg LinAlgBackend
double myLambda
Lambda parameter.
PolygonalCalculus::SparseMatrix SparseMatrix
Vertex myLastSourceIndex
Vertex index to the last source point (to shift the distances)
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
GeodesicsInHeat(const GeodesicsInHeat &other)=delete
void init(double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
Solver myHeatSolver
Heat solver.
bool myIsInit
Validitate flag.
SparseMatrix myHeatOpe
The operator for heat diffusion.
PolygonalCalculus::Vertex Vertex
void addSource(const Vertex aV)
IntegerVector myBoundary
The boundary characteristic vector.
Solver myPoissonSolver
Poisson solver.
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
PolygonalCalculus::Vector Vector
TPolygonalCalculus PolygonalCalculus
PolygonalCalculus::Solver Solver
GeodesicsInHeat(ConstAlias< PolygonalCalculus > calculus)
PolygonalCalculus::DenseMatrix DenseMatrix
Vector mySource
Source vector.
GeodesicsInHeat & operator=(const GeodesicsInHeat &other)=delete
Conditions::IntegerVector IntegerVector
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
MySurfaceMesh::Vertex Vertex
Vertex type.
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
LinAlg::DenseVector Vector
Type of Vector.
SurfMesh surfmesh
PolyCalculus * calculus
DGtal is the top-level namespace which contains all DGtal functions and types.
void laplacian(Shape &shape, const Options &options, std::function< double(const RealPoint3D &)> input_function, std::function< double(const RealPoint3D &)> target_function, int argc, char **argv)
Aim: Provide linear algebra backend using Eigen dense and sparse matrix as well as dense vector....
Definition: EigenSupport.h:96
std::size_t Index
The type used for numbering vertices and faces.
Definition: SurfaceMesh.h:105