DGtal 1.3.0
Searching...
No Matches
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
44using namespace DGtal;
45using namespace Z3i;
46
47// Using standard 3D digital space.
50// The following typedefs are useful
52typedef SurfMesh::Face Face;
54
55//Polyscope global
56polyscope::SurfaceMesh *psMesh;
57SurfMesh surfmesh;
58float scale = 0.1;
60
61//Restriction of a scalar function to vertices
62double 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
82void initPhi()
83{
84 phiEigen.resize(surfmesh.nbVertices());
85 for(auto i = 0; i < surfmesh.nbVertices(); ++i)
86 phiEigen(i) = phiVertex(i);
88}
89
90void 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 {
107
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
125}
126
127
128void 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 {
145
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
163}
164
165
166void 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
204}
205
206void 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
223int main()
224{
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}
Aim: A helper class to solve a system with Dirichlet boundary conditions.
Implements differential operators on polygonal surfaces from .
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Real3dVector faceNormalAsDGtalVector(const Face f) const
Vector vectorArea(const Face f) const
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
double faceArea(const Face f) const
LinAlg::DenseVector Vector
Type of Vector.
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
static Parameters parametersGeometryEstimation()
static Parameters defaultParameters()
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
static KSpace getKSpace(const Point &low, const Point &up, Parameters params=parametersKSpace())
Definition: Shortcuts.h:332
std::map< Cell, IdxVertex > Cell2Index
Definition: Shortcuts.h:189
static CountedPtr< DigitalSurface > makeDigitalSurface(CountedPtr< TPointPredicate > bimage, const KSpace &K, const Parameters &params=parametersDigitalSurface())
Definition: Shortcuts.h:1209
static Parameters defaultParameters()
Definition: Shortcuts.h:203
static CountedPtr< SurfaceMesh > makePrimalSurfaceMesh(Cell2Index &c2i, CountedPtr< ::DGtal::DigitalSurface< TContainer > > aSurface)
Definition: Shortcuts.h:2372
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
void beginBlock(const std::string &keyword="")
double endBlock()
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:154
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
const Vertices & incidentVertices(Face f) const
Definition: SurfaceMesh.h:307
Size nbFaces() const
Definition: SurfaceMesh.h:288
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
Size nbVertices() const
Definition: SurfaceMesh.h:280
Edges computeNonManifoldEdges() const
int main(int argc, char **argv)
KSpace K
TriMesh::Face Face
TriMesh::Vertex Vertex