DGtal  1.3.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)
34 
35 #define GeodesicsInHeat_RECURSES
36 
37 #if !defined GeodesicsInHeat_h
38 
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 
118  GeodesicsInHeat & operator= ( GeodesicsInHeat && other ) = delete;
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 I(myCalculus->nbVertices(),myCalculus->nbVertices());
142  I.setIdentity();
143  SparseMatrix laplacian = myCalculus->globalLaplaceBeltrami( lambda );
144  SparseMatrix mass = myCalculus->globalLumpedMassMatrix();
145  myHeatOpe = mass - dt*laplacian;
146 
147  //Prefactorizing
148  myPoissonSolver.compute( laplacian );
149  myHeatSolver.compute ( myHeatOpe );
150 
151  //empty source
152  mySource = Vector::Zero(myCalculus->nbVertices());
153 
154  // Manage boundaries
155  myManageBoundary = false;
156  if ( ! boundary_with_mixed_solution ) return;
157  myBoundary = IntegerVector::Zero(myCalculus->nbVertices());
158  const auto surfmesh = myCalculus->getSurfaceMeshPtr();
159  const auto edges = surfmesh->computeManifoldBoundaryEdges();
160  for ( auto e : edges )
161  {
162  const auto vtcs = surfmesh->edgeVertices( e );
163  myBoundary[ vtcs.first ] = 1;
164  myBoundary[ vtcs.second ] = 1;
165  }
166  myManageBoundary = ! edges.empty();
167  if ( ! myManageBoundary ) return;
168  // Prepare solver for a problem with Dirichlet conditions.
170  // Prefactoring
171  myHeatDirichletSolver.compute( heatOpe_d );
172  }
173 
177  void addSource(const Vertex aV)
178  {
179  ASSERT_MSG(aV < myCalculus->nbVertices(), "Vertex is not in the surface mesh vertex range");
180  myLastSourceIndex = aV;
181  mySource( aV ) = 1.0;
182  }
183 
187  Vector source() const
188  {
189  FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
190  return mySource;
191  }
192 
193 
196  Vector compute() const
197  {
198  FATAL_ERROR_MSG(myIsInit, "init() method must be called first");
199  //Heat diffusion
200  Vector heatDiffusion = myHeatSolver.solve(mySource);
201  // Take care of boundaries
202  if ( myManageBoundary )
203  {
204  Vector bValues = Vector::Zero( myCalculus->nbVertices() );
206  myBoundary, bValues );
207  Vector bSol = myHeatDirichletSolver.solve( bSources );
208  Vector heatDiffusionDirichlet
209  = Conditions::dirichletSolution( bSol, myBoundary, bValues );
210  heatDiffusion = 0.5 * ( heatDiffusion + heatDiffusionDirichlet );
211  }
212 
213  Vector divergence = Vector::Zero(myCalculus->nbVertices());
214  auto cpt=0;
215  auto surfmesh = myCalculus->getSurfaceMeshPtr();
216 
217  // Heat, normalization and divergence per face
218  for(auto f=0; f< myCalculus->nbFaces(); ++f)
219  {
220  Vector faceHeat( myCalculus->degree(f));
221  cpt=0;
222  auto vertices = surfmesh->incidentVertices(f);
223  for(auto v: vertices)
224  {
225  faceHeat(cpt) = heatDiffusion( v );
226  ++cpt;
227  }
228  // ∇heat / ∣∣∇heat∣∣
229  Vector grad = -myCalculus->gradient(f) * faceHeat;
230  grad.normalize();
231 
232  // div
233  DenseMatrix oneForm = myCalculus->flat(f)*grad;
234  Vector divergenceFace = myCalculus->divergence( f, myLambda ) * oneForm;
235  cpt=0;
236  for(auto v: vertices)
237  {
238  divergence(v) += divergenceFace(cpt);
239  ++cpt;
240  }
241  }
242 
243  // Last Poisson solve
244  Vector distVec = myPoissonSolver.solve(divergence);
245 
246  //Source val
247  auto sourceval = distVec(myLastSourceIndex);
248 
249  //shifting the distances to get 0 at sources
250  return distVec - sourceval*Vector::Ones(myCalculus->nbVertices());
251  }
252 
253 
255  bool isValid() const
256  {
257  return myIsInit && myCalculus->isValid();
258  }
259 
260  // ----------------------- Private --------------------------------------
261 
262  private:
263 
266 
269 
272 
275 
278 
281 
283  bool myIsInit;
284 
286  double myLambda;
287 
291 
294 
297 
298 
299  }; // end of class GeodesicsInHeat
300 } // namespace DGtal
301 
302 // //
304 
305 #endif // !defined GeodesicsInHeat_h
306 
307 #undef GeodesicsInHeat_RECURSES
308 #endif // else defined(GeodesicsInHeat_RECURSES)
DGtal::GeodesicsInHeat::Conditions
DirichletConditions< LinAlgBackend > Conditions
Definition: GeodesicsInHeat.h:74
DGtal::GeodesicsInHeat::myBoundary
IntegerVector myBoundary
The boundary characteristic vector.
Definition: GeodesicsInHeat.h:293
DGtal::GeodesicsInHeat::GeodesicsInHeat
GeodesicsInHeat()=delete
DGtal::ConstAlias
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
DGtal::GeodesicsInHeat::myCalculus
const PolygonalCalculus * myCalculus
The underlying PolygonalCalculus instance.
Definition: GeodesicsInHeat.h:265
DGtal::GeodesicsInHeat::DenseMatrix
PolygonalCalculus::DenseMatrix DenseMatrix
Definition: GeodesicsInHeat.h:69
DGtal::DirichletConditions::dirichletVector
static DenseVector dirichletVector(const SparseMatrix &A, const DenseVector &b, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:170
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:103
DGtal::GeodesicsInHeat::PolygonalCalculus
TPolygonalCalculus PolygonalCalculus
Definition: GeodesicsInHeat.h:67
DGtal::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
DGtal::GeodesicsInHeat::source
Vector source() const
Definition: GeodesicsInHeat.h:187
DGtal::GeodesicsInHeat::myHeatOpe
SparseMatrix myHeatOpe
The operator for heat diffusion.
Definition: GeodesicsInHeat.h:268
DGtal::DirichletConditions::dirichletSolution
static DenseVector dirichletSolution(const DenseVector &xd, const IntegerVector &p, const DenseVector &u)
Definition: DirichletConditions.h:203
DGtal::DirichletConditions::dirichletOperator
static SparseMatrix dirichletOperator(const SparseMatrix &A, const IntegerVector &p)
Definition: DirichletConditions.h:130
DGtal::GeodesicsInHeat::LinAlgBackend
PolygonalCalculus::LinAlg LinAlgBackend
Definition: GeodesicsInHeat.h:73
DGtal::PolygonalCalculus::DenseMatrix
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
Definition: PolygonalCalculus.h:101
DGtal::GeodesicsInHeat::mySource
Vector mySource
Source vector.
Definition: GeodesicsInHeat.h:277
laplacian
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)
Definition: sphereCotangentLaplaceOperator.cpp:89
DGtal::GeodesicsInHeat::GeodesicsInHeat
GeodesicsInHeat(ConstAlias< PolygonalCalculus > calculus)
Definition: GeodesicsInHeat.h:84
DGtal::GeodesicsInHeat::myManageBoundary
bool myManageBoundary
Definition: GeodesicsInHeat.h:290
DGtal::GeodesicsInHeat::Vertex
PolygonalCalculus::Vertex Vertex
Definition: GeodesicsInHeat.h:72
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::EigenLinearAlgebraBackend
Aim: Provide linear algebra backend using Eigen dense and sparse matrix as well as dense vector....
Definition: EigenSupport.h:97
DGtal::SurfaceMesh::edgeVertices
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
DGtal::GeodesicsInHeat::compute
Vector compute() const
Definition: GeodesicsInHeat.h:196
DGtal::PolygonalCalculus::Solver
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
Definition: PolygonalCalculus.h:108
DGtal::GeodesicsInHeat::SparseMatrix
PolygonalCalculus::SparseMatrix SparseMatrix
Definition: GeodesicsInHeat.h:68
DGtal::GeodesicsInHeat
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
Definition: GeodesicsInHeat.h:62
DGtal::DirichletConditions::IntegerVector
LinearAlgebraBackend::IntegerVector IntegerVector
Definition: DirichletConditions.h:105
DGtal::GeodesicsInHeat::addSource
void addSource(const Vertex aV)
Definition: GeodesicsInHeat.h:177
DGtal::GeodesicsInHeat::IntegerVector
Conditions::IntegerVector IntegerVector
Definition: GeodesicsInHeat.h:75
DGtal::GeodesicsInHeat::operator=
GeodesicsInHeat & operator=(const GeodesicsInHeat &other)=delete
DGtal::GeodesicsInHeat::~GeodesicsInHeat
~GeodesicsInHeat()=default
DGtal::PolygonalCalculus::Vertex
MySurfaceMesh::Vertex Vertex
Vertex type.
Definition: PolygonalCalculus.h:85
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::GeodesicsInHeat::myHeatDirichletSolver
Solver myHeatDirichletSolver
Heat solver with Dirichlet boundary conditions.
Definition: GeodesicsInHeat.h:296
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:97
DGtal::GeodesicsInHeat::myHeatSolver
Solver myHeatSolver
Heat solver.
Definition: GeodesicsInHeat.h:274
DGtal::GeodesicsInHeat::isValid
bool isValid() const
Definition: GeodesicsInHeat.h:255
DGtal::GeodesicsInHeat::init
void init(double dt, double lambda=1.0, bool boundary_with_mixed_solution=false)
Definition: GeodesicsInHeat.h:135
DGtal::GeodesicsInHeat::myLambda
double myLambda
Lambda parameter.
Definition: GeodesicsInHeat.h:286
DGtal::GeodesicsInHeat::myLastSourceIndex
Vertex myLastSourceIndex
Vertex index to the last source point (to shift the distances)
Definition: GeodesicsInHeat.h:280
DGtal::SurfaceMesh::incidentVertices
const Vertices & incidentVertices(Face f) const
Definition: SurfaceMesh.h:307
DGtal::GeodesicsInHeat::myIsInit
bool myIsInit
Validitate flag.
Definition: GeodesicsInHeat.h:283
DGtal::GeodesicsInHeat::myPoissonSolver
Solver myPoissonSolver
Poisson solver.
Definition: GeodesicsInHeat.h:271
DGtal::GeodesicsInHeat::Solver
PolygonalCalculus::Solver Solver
Definition: GeodesicsInHeat.h:70
DGtal::GeodesicsInHeat::Vector
PolygonalCalculus::Vector Vector
Definition: GeodesicsInHeat.h:71