DGtal  1.5.beta
exampleVectorHeatMethod.cpp
1 
27 
28 #include <iostream>
29 
30 #include <DGtal/base/Common.h>
31 #include <DGtal/helpers/StdDefs.h>
32 #include <DGtal/helpers/Shortcuts.h>
33 #include <DGtal/helpers/ShortcutsGeometry.h>
34 #include <DGtal/shapes/SurfaceMesh.h>
35 #include <DGtal/dec/PolygonalCalculus.h>
36 #include <DGtal/dec/VectorsInHeat.h>
37 
38 #include <polyscope/polyscope.h>
39 #include <polyscope/surface_mesh.h>
40 
42 
43 using namespace DGtal;
44 
45 // Using standard 3D digital space.
48 // The following typedefs are useful
52 typedef SurfMesh::Face Face;
53 typedef SurfMesh::Vertex Vertex;
55 typedef PC::Vector Vector;
57 
58 //Polyscope global
59 polyscope::SurfaceMesh *psMesh;
60 
62 
63 //Polygonal Calculus and VectorsInHeat solvers
64 PC *calculus;
66 
67 //sources
68 std::vector<Vector> X_0;
69 std::vector<Vertex> idX_0;
70 
71 //rotation matrix
73 
74 bool noSources = true;
75 bool toggle=false;
76 
81 void addRandomSource()
82 {
83  size_t id = rand()%surfmesh.nbVertices();
84  VHM->addSource(id,Eigen::Vector3d::Random(3).normalized());
85 
86  idX_0.push_back(id);
87  X_0[id] = VHM->extrinsicVectorSourceAtVertex(id);
88  psMesh->addVertexVectorQuantity("X_0",X_0);
89  noSources = false;
90 }
91 
96 void diffuse()
97 {
98  if (noSources)
99  addRandomSource();
100  psMesh->addVertexVectorQuantity("VHM field",VHM->compute());
101 }
102 
106 void precompute()
107 {
108  auto nv = surfmesh.nbVertices();
109  auto ael = surfmesh.averageEdgeLength();
110  VHM->init(ael);//init vector heat method solvers (should be ael^2 but smoother results that way)
111 
112  X_0.resize(nv,Vector::Zero(3));//extrinsic Source vectors
113 
114  psMesh->addVertexVectorQuantity("X_0",X_0);
115 }
116 
117 void clearSources()
118 {
119  VHM->clearSource();
120  noSources=true;
121  idX_0.clear();
122  X_0.clear();
123  X_0.resize(surfmesh.nbVertices(),Vector::Zero(3));
124  //cleanup the visualization
125  psMesh->addVertexVectorQuantity("X_0",X_0);
126  psMesh->addVertexVectorQuantity("VHM field",X_0);
127 }
128 
129 void rotate()
130 {
131  VHM->clearSource();
132  for(const auto id: idX_0)
133  {
134  Vector x = R*X_0[id];
135  VHM->addSource(id, x);
136  X_0[id] = VHM->extrinsicVectorSourceAtVertex(id);
137  }
138  psMesh->addVertexVectorQuantity("X_0",X_0);
139 }
140 
141 
142 void myCallback()
143 {
144  if(ImGui::Button("Compute Vector Field"))
145  {
146  diffuse();
147  }
148  if(ImGui::Button("Add random source"))
149  {
150  addRandomSource();
151  }
152  if(ImGui::Button("Clear sources"))
153  {
154  clearSources();
155  }
156 
157  if(ImGui::Button("Start/stop rotating sources"))
158  toggle = !toggle;
159 
160 
161  if (toggle)
162  {
163  rotate();
164  diffuse();
165  }
166 }
167 
168 int main(int argc, char **argv)
169 {
170  std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
171  std::vector<RealPoint> positions;
172 
173  if (argc <= 1)
174  {
175  trace.error()<<"Missing vol file. Usage: exampleVectorHeatMethod bunny.vol"<<std::endl;
176  exit(2);
177  }
178 
179  //load voxel model
180  auto params = SH3::defaultParameters() | SHG3::defaultParameters() | SHG3::parametersGeometryEstimation();
181  params("surfaceComponents", "All");
182  auto binary_image = SH3::makeBinaryImage(argv[1], params );
183  auto K = SH3::getKSpace( binary_image, params );
184  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
185  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
186 
187  //Need to convert the faces
188  for(size_t face= 0 ; face < primalSurface->nbFaces(); ++face)
189  faces.push_back(primalSurface->incidentVertices( face ));
190 
191  //Recasting to vector of vertices
192  positions = primalSurface->positions();
193 
194  surfmesh = SurfMesh(positions.begin(),
195  positions.end(),
196  faces.begin(),
197  faces.end());
198 
199  //instantiate PolyDEC
200  calculus = new PC(surfmesh);
201 
202  //instantiate VHM
203  VHM = new VectorsInHeat<PC>(calculus);
204 
205  //For the rotation of input VF
206  PC::DenseMatrix Rotx(3,3),Roty(3,3),Rotz(3,3);
207  double theta=0.05;
208  Rotx << 1, 0,0,0,cos(theta),-sin(theta),0,sin(theta),cos(theta);
209  Roty << cos(theta), 0,sin(theta),0,1,0,-sin(theta),0,cos(theta);
210  Rotz << cos(theta), -sin(theta),0,sin(theta),cos(theta),0,0,0,1;
211  R=Rotz*Roty*Rotx;
212 
213  //Initialize polyscope
214  polyscope::init();
215 
216  psMesh = polyscope::registerSurfaceMesh("Digital Surface", positions, faces);
217 
218  //Initialize solvers
219  precompute();
220 
221  polyscope::view::upDir = polyscope::view::UpDir::XUp;
222 
223  polyscope::state::userCallback = myCallback;
224 
225  polyscope::show();
226  return EXIT_SUCCESS;
227 
228 }
Implements differential operators on polygonal surfaces from .
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
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...
Aim: This class is used to simplify shape and surface creation. With it, you can create new shapes an...
Definition: Shortcuts.h:105
std::ostream & error()
This class implements on polygonal surfaces (using Discrete differential calculus on polygonal surfa...
Definition: VectorsInHeat.h:63
SurfaceMesh< RealPoint, RealVector > SurfMesh
void myCallback()
polyscope::SurfaceMesh * psMesh
SurfMesh surfmesh
void precompute()
PolyCalculus * calculus
void clearSources()
CountedPtr< SH3::DigitalSurface > surface
CountedPtr< SH3::BinaryImage > binary_image
DigitalPlane::Point Vector
SMesh::Vertices Vertices
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:153
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
std::vector< Vertex > Vertices
The type that defines a list/range of vertices (e.g. to define faces)
Definition: SurfaceMesh.h:112
TRealPoint RealPoint
Definition: SurfaceMesh.h:93
Size nbVertices() const
Definition: SurfaceMesh.h:288
Scalar averageEdgeLength() const
int main(int argc, char **argv)
Shortcuts< KSpace > SH3
K init(Point(0, 0, 0), Point(512, 512, 512), true)
KSpace K
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
ShortcutsGeometry< Z3i::KSpace > SHG3
TriMesh::Face Face
PointVector< 3, double > RealPoint
TriMesh::Vertex Vertex