DGtal  1.4.beta
dgtalCalculus-bunny.cpp
1 
25 #include <iostream>
26 #include <string>
27 #include <DGtal/base/Common.h>
28 #include <DGtal/helpers/StdDefs.h>
29 #include <DGtal/helpers/Shortcuts.h>
30 #include <DGtal/helpers/ShortcutsGeometry.h>
31 #include <DGtal/shapes/SurfaceMesh.h>
32 #include <DGtal/geometry/surfaces/DigitalSurfaceRegularization.h>
33 #include <DGtal/dec/PolygonalCalculus.h>
34 #include <DGtal/math/linalg/DirichletConditions.h>
35 
36 #include <polyscope/polyscope.h>
37 #include <polyscope/surface_mesh.h>
38 #include <polyscope/point_cloud.h>
39 #include <Eigen/Dense>
40 #include <Eigen/Sparse>
41 
42 #include "ConfigExamples.h"
43 
44 using namespace DGtal;
45 using namespace Z3i;
46 
47 // Using standard 3D digital space.
50 // The following typedefs are useful
52 typedef SurfMesh::Face Face;
53 typedef SurfMesh::Vertex Vertex;
54 
55 //Polyscope global
56 polyscope::SurfaceMesh *psMesh;
57 SurfMesh surfmesh;
58 float scale = 0.1;
60 
61 //Restriction of a scalar function to vertices
62 double phiVertex(const Vertex v)
63 {
64  return cos(scale*(surfmesh.position(v)[0]))*sin(scale*surfmesh.position(v)[1]);
65 }
66 
67 //Restriction of a scalar function to vertices
69 {
70  auto vertices = surfmesh.incidentVertices(f);
71  auto nf = vertices.size();
72  Eigen::VectorXd ph(nf);
73  size_t cpt=0;
74  for(auto v: vertices)
75  {
76  ph(cpt) = phiVertex(v);
77  ++cpt;
78  }
79  return ph;
80 }
81 
82 void initPhi()
83 {
84  phiEigen.resize(surfmesh.nbVertices());
85  for(auto i = 0; i < surfmesh.nbVertices(); ++i)
86  phiEigen(i) = phiVertex(i);
87  psMesh->addVertexScalarQuantity("Phi", phiEigen);
88 }
89 
90 void initQuantities()
91 {
92  trace.beginBlock("Basic quantities");
94 
97  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> normals;
98  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> vectorArea;
99  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> centroids;
100 
101  std::vector<double> faceArea;
102 
103  for(auto f=0; f < surfmesh.nbFaces(); ++f)
104  {
106  gradients.push_back( grad );
107 
109  cogradients.push_back( cograd );
110 
111  normals.push_back(calculus.faceNormalAsDGtalVector(f));
112 
114  vectorArea.push_back({vA(0) , vA(1), vA(2)});
115 
116  faceArea.push_back( calculus.faceArea(f));
117  }
118  trace.endBlock();
119 
120  psMesh->addFaceVectorQuantity("Gradients", gradients);
121  psMesh->addFaceVectorQuantity("co-Gradients", cogradients);
122  psMesh->addFaceVectorQuantity("Normals", normals);
123  psMesh->addFaceScalarQuantity("Face area", faceArea);
124  psMesh->addFaceVectorQuantity("Vector area", vectorArea);
125 }
126 
127 
128 void initQuantitiesCached()
129 {
130  trace.beginBlock("Basic quantities (cached)");
132 
135  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> normals;
136  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> vectorArea;
137  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> centroids;
138 
139  std::vector<double> faceArea;
140 
141  for(auto f=0; f < surfmesh.nbFaces(); ++f)
142  {
144  gradients.push_back( grad );
145 
147  cogradients.push_back( cograd );
148 
149  normals.push_back(calculus.faceNormalAsDGtalVector(f));
150 
152  vectorArea.push_back({vA(0) , vA(1), vA(2)});
153 
154  faceArea.push_back( calculus.faceArea(f));
155  }
156  trace.endBlock();
157 
158  psMesh->addFaceVectorQuantity("Gradients", gradients);
159  psMesh->addFaceVectorQuantity("co-Gradients", cogradients);
160  psMesh->addFaceVectorQuantity("Normals", normals);
161  psMesh->addFaceScalarQuantity("Face area", faceArea);
162  psMesh->addFaceVectorQuantity("Vector area", vectorArea);
163 }
164 
165 
166 void computeLaplace()
167 {
169  trace.beginBlock("Operator construction...");
171  trace.endBlock();
172 
173  const auto nbv = surfmesh.nbVertices();
175 
176  //Setting some random sources
178  DC::IntegerVector p = DC::nullBoundaryVector( g );
179  for(auto cpt=0; cpt< 10;++cpt)
180  {
181  int idx = rand() % nbv;
182  g( idx ) = rand() % 100;
183  p( idx ) = 1.0;
184  }
185 
186  //Solve Δu=0 with g as boundary conditions
188 
189  trace.beginBlock("Prefactorization...");
190  DC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, p );
191  solver.compute( L_dirichlet );
192  ASSERT(solver.info()==Eigen::Success);
193  trace.endBlock();
194 
195  trace.beginBlock("Solve...");
196  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector g_dirichlet = DC::dirichletVector( L, g, p, g );
197  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector x_dirichlet = solver.solve( g_dirichlet );
198  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector u = DC::dirichletSolution( x_dirichlet, p, g );
199  ASSERT(solver.info()==Eigen::Success);
200  trace.endBlock();
201 
202  psMesh->addVertexScalarQuantity("g", g);
203  psMesh->addVertexScalarQuantity("u", u);
204 }
205 
206 void myCallback()
207 {
208  ImGui::SliderFloat("Phi scale", &scale, 0., 1.);
209  if (ImGui::Button("Phi and basic operators"))
210  {
211  initPhi();
212  initQuantities();
213  }
214  if (ImGui::Button("Phi and basic operators (cached)"))
215  {
216  initPhi();
217  initQuantitiesCached();
218  }
219  if(ImGui::Button("Compute Laplace problem"))
220  computeLaplace();
221 }
222 
223 int main()
224 {
225  auto params = SH3::defaultParameters() | SHG3::defaultParameters() | SHG3::parametersGeometryEstimation();
226  params("surfaceComponents", "All");
227 
228  std::string filename = examplesPath + std::string("/samples/bunny-32.vol");
229  auto binary_image = SH3::makeBinaryImage(filename, params );
230  auto K = SH3::getKSpace( binary_image, params );
231  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
232  SH3::Cell2Index c2i;
233  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
234 
235  //Need to convert the faces
236  std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
237  std::vector<RealPoint> positions;
238 
239  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
240  faces.push_back(primalSurface->incidentVertices( face ));
241 
242  //Recasting to vector of vertices
243  positions = primalSurface->positions();
244 
245  surfmesh = SurfMesh(positions.begin(),
246  positions.end(),
247  faces.begin(),
248  faces.end());
249 
250  std::cout<<"number of non-manifold Edges = " << surfmesh.computeNonManifoldEdges().size()<<std::endl;
251  psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
252 
253  // Initialize polyscope
254  polyscope::init();
255 
256  // Set the callback function
257  polyscope::state::userCallback = myCallback;
258  polyscope::show();
259  return EXIT_SUCCESS;
260 }
DGtal::SurfaceMesh::Vertex
Index Vertex
Definition: SurfaceMesh.h:108
DGtal::Shortcuts::Cell2Index
std::map< Cell, IdxVertex > Cell2Index
Definition: Shortcuts.h:189
DGtal::PolygonalCalculus::vectorArea
Vector vectorArea(const Face f) const
Definition: PolygonalCalculus.h:316
DGtal::SurfaceMesh::position
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
DGtal::Trace::endBlock
double endBlock()
SH3
Shortcuts< KSpace > SH3
Definition: testArithmeticalDSSComputerOnSurfels.cpp:49
DGtal::PolygonalCalculus::faceNormalAsDGtalVector
Real3dVector faceNormalAsDGtalVector(const Face f) const
Definition: PolygonalCalculus.h:358
DGtal::PolygonalCalculus::globalLaplaceBeltrami
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Definition: PolygonalCalculus.h:798
Face
TriMesh::Face Face
Definition: testTriangulatedSurface.cpp:56
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:102
DGtal::trace
Trace trace
Definition: Common.h:154
K
KSpace K
Definition: testCubicalComplex.cpp:62
DGtal::SurfaceMesh
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:91
DGtal::Shortcuts::makeDigitalSurface
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
DGtal::PolygonalCalculus::faceArea
double faceArea(const Face f) const
Definition: PolygonalCalculus.h:340
DGtal::Trace::beginBlock
void beginBlock(const std::string &keyword="")
DGtal::Shortcuts::makePrimalSurfaceMesh
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:70
DGtal::Shortcuts
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:104
DGtal::PolygonalCalculus::coGradient
DenseMatrix coGradient(const Face f) const
Definition: PolygonalCalculus.h:367
Vertex
TriMesh::Vertex Vertex
Definition: testTriangulatedSurface.cpp:57
DGtal::SurfaceMesh::nbVertices
Size nbVertices() const
Definition: SurfaceMesh.h:280
boost::vertices
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
SparseMatrix
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: testHeatLaplace.cpp:50
DGtal::SurfaceMesh::computeNonManifoldEdges
Edges computeNonManifoldEdges() const
SHG3
ShortcutsGeometry< Z3i::KSpace > SHG3
Definition: testIntegralInvariantShortcuts.cpp:47
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::PolygonalCalculus::gradient
DenseMatrix gradient(const Face f) const
Definition: PolygonalCalculus.h:390
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::PolygonalCalculus::Solver
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
Definition: PolygonalCalculus.h:107
DGtal::ShortcutsGeometry
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: ShortcutsGeometry.h:76
DGtal::SurfaceMesh::Face
Index Face
Definition: SurfaceMesh.h:106
Vector
FreemanChain< int >::Vector Vector
Definition: testCombinDSS.cpp:60
main
int main(int argc, char **argv)
Definition: testArithmeticDSS-benchmark.cpp:147
DGtal::Shortcuts::getKSpace
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
DGtal::DirichletConditions
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Definition: DirichletConditions.h:97
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:96
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
DGtal::SurfaceMesh::incidentVertices
const Vertices & incidentVertices(Face f) const
Definition: SurfaceMesh.h:307
init
K init(Point(0, 0, 0), Point(512, 512, 512), true)
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288