DGtal 1.4.2
Loading...
Searching...
No Matches
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
49namespace 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
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");
187 mySource( aV ) = 1.0;
188 }
189
193 {
194 mySource = Vector::Zero(myCalculus->nbVertices());
195 }
196
201 {
202 FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
203 return mySource;
204 }
205
206
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
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 crane2013 on polygonal surfaces (using Discrete differential calculus on polygo...
GeodesicsInHeat(GeodesicsInHeat &&other)=delete
PolygonalCalculus::DenseMatrix DenseMatrix
double myLambda
Lambda parameter.
Vertex myLastSourceIndex
Vertex index to the last source point (to shift the distances)
TPolygonalCalculus PolygonalCalculus
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
DirichletConditions< LinAlgBackend > Conditions
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.
void addSource(const Vertex aV)
IntegerVector myBoundary
The boundary characteristic vector.
Solver myPoissonSolver
Poisson solver.
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
PolygonalCalculus::LinAlg LinAlgBackend
GeodesicsInHeat & operator=(const GeodesicsInHeat &other)=delete
PolygonalCalculus::Vertex Vertex
PolygonalCalculus::SparseMatrix SparseMatrix
PolygonalCalculus::Solver Solver
GeodesicsInHeat(ConstAlias< PolygonalCalculus > calculus)
Vector mySource
Source vector.
PolygonalCalculus::Vector Vector
Conditions::IntegerVector IntegerVector
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
MySurfaceMesh::Vertex Vertex
Vertex type.
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
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....
std::size_t Index
The type used for numbering vertices and faces.