DGtal  1.4.beta
exampleHarmonicParametrization.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/math/linalg/DirichletConditions.h"
37 
38 #include <polyscope/polyscope.h>
39 #include <polyscope/surface_mesh.h>
40 
41 #include "ConfigExamples.h"
42 
44 
45 using namespace DGtal;
46 
48 // Using standard 3D digital space.
51 // The following typedefs are useful
55 typedef SurfMesh::Face Face;
56 typedef SurfMesh::Vertex Vertex;
60 typedef PC::Solver Solver;
61 typedef PC::Vector Vector;
62 typedef PC::Triplet Triplet;
63 typedef std::vector<Vertex> chain;
65 typedef Conditions::IntegerVector IntegerVector;
66 
67 
68 //Polyscope global
69 polyscope::SurfaceMesh *psMesh;
70 polyscope::SurfaceMesh *psParam;
71 
72 //DEC
73 PC *calculus;
74 
75 //surface mesh global
76 SurfMesh surfmesh;
77 std::vector<std::vector<size_t>> faces;
78 std::vector<RealPoint> positions;
79 
80 
87 std::pair<Vector,Vector> FixBoundaryParametrization(const std::vector<Vertex>& boundary)
88 {
89  auto nb = boundary.size();
90  auto n = surfmesh.nbVertices();
91  Vector u = Vector::Zero(n),v = Vector::Zero(n);
92  double totalBoundaryLength = 0;
93  for (Vertex i = 0;i<nb;i++)
94  totalBoundaryLength += surfmesh.distance(boundary[(i+1)%nb],boundary[i]);
95 
96  double partialSum = 0;
97  for (Vertex i = 0;i<nb;i++)
98  {
99  double th = 2*M_PI*partialSum/totalBoundaryLength;
100  auto vi = boundary[i];
101  auto vj = boundary[(i+1)%nb];
102  u(vi) = std::cos(th);
103  v(vj) = std::sin(th);
104  partialSum += surfmesh.distance(vi,vj);
105  }
106  return {u,v};
107 }
108 
114 void VisualizeParametrizationOnCircle(const DenseMatrix& UV)
115 {
116  auto n= surfmesh.nbVertices();
117  std::vector<RealPoint> pos(n);
118  double scale = 0;
119  RealPoint avg = {0.,0.,0.};
120  for (size_t v = 0;v<n;v++)
121  {
122  auto p = surfmesh.position(v);
123  avg += p;
124  if (p.norm() > scale)
125  scale = p.norm();
126  }
127  avg /= surfmesh.nbVertices();
128  scale /= 2;
129  for (size_t v = 0;v<n;v++)
130  pos[v] = RealPoint{scale*UV(v,0),scale*UV(v,1),0.} + avg;
131  polyscope::registerSurfaceMesh("On circle parametrization", pos, faces)->setEnabled(false);
132 }
133 
140 DenseMatrix HarmonicParametrization()
141 {
142  auto n = surfmesh.nbVertices();
143  std::cout<<"Nb boundary edges = "<< surfmesh.computeManifoldBoundaryEdges().size()<<std::endl;
144  std::vector<chain> chains = surfmesh.computeManifoldBoundaryChains();
145  //choose longest chain as boundary of the parametrization
146  std::cout<<"Nb boundaries = "<< chains.size() << std::endl;
147 
148  //choose longest chain as boundary of the parametrization
149  auto B = *std::max_element(chains.begin(),chains.end(),[] (const chain& A,const chain& B) {return A.size() < B.size();});
150 
151  IntegerVector boundary = IntegerVector::Zero(n);
152  for (Vertex v : B)
153  boundary(v) = 1;
154 
155  std::pair<Vector,Vector> uv_b = FixBoundaryParametrization(B);//maps boundary to circle
156 
157  calculus = new PC(surfmesh);
158 
159  //Impose dirichlet boundary condition to laplace problem
160  Vector Z = Vector::Zero(n);
161  SparseMatrix L = calculus->globalLaplaceBeltrami();
162  SparseMatrix L_d = Conditions::dirichletOperator( L, boundary );
163 
164  PC::Solver solver;
165  solver.compute(L_d);
166 
167  Vector b_u = Conditions::dirichletVector( L, Z,boundary, uv_b.first );
168  Vector b_v = Conditions::dirichletVector( L, Z,boundary, uv_b.second );
169 
170  Vector rslt_u_d = solver.solve(b_u);
171  Vector rslt_v_d = solver.solve(b_v);
172 
173  Vector rslt_u = Conditions::dirichletSolution(rslt_u_d,boundary,uv_b.first);
174  Vector rslt_v = Conditions::dirichletSolution(rslt_v_d,boundary,uv_b.second);
175 
176  DenseMatrix uv(n,2);
177  uv.col(0) = rslt_u;
178  uv.col(1) = rslt_v;
179 
180  return uv;
181 }
182 
183 int main(int argc, char **argv)
184 {
185  //Import Voxel Model
186  auto params = SH3::defaultParameters() | SHG3::defaultParameters() | SHG3::parametersGeometryEstimation();
187  params("surfaceComponents", "All");
188 
189  //load .vol
190  std::string inputFilename(examplesPath + "samples/bunny-64.vol" );
191 
192  auto binary_image = SH3::makeBinaryImage(inputFilename, params );
193 
194  //offset K space to create boundary
195  auto K = SH3::getKSpace( binary_image, params );
196  K.init(K.lowerBound()+SH3::Point(5,0,0),K.upperBound(),true);
197 
198  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
199  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
200 
201 
202  //Need to convert the faces
203  for(size_t face= 0 ; face < primalSurface->nbFaces(); ++face)
204  faces.push_back(primalSurface->incidentVertices( face ));
205 
206  //Recasting to vector of vertices
207  positions = primalSurface->positions();
208 
209  surfmesh = SurfMesh(positions.begin(),
210  positions.end(),
211  faces.begin(),
212  faces.end());
213 
214 
215  // Initialize polyscope
216  polyscope::init();
217 
218  //compute parametrization
219  DenseMatrix UV = HarmonicParametrization();
220 
221  //visualize parametrization on 2D circle
222  VisualizeParametrizationOnCircle(UV);
223 
224  psMesh = polyscope::registerSurfaceMesh("Digital Surface", positions, faces);
225  psMesh->addVertexParameterizationQuantity("Harmonic parametrization",UV)->setEnabled(true);
226 
227  //set correct view for voxel mesh
228  polyscope::view::upDir = polyscope::view::UpDir::XUp;
229 
230  polyscope::show();
231  return EXIT_SUCCESS;
232 }
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)
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
const Point & lowerBound() const
Return the lower bound for digital points in this space.
const Point & upperBound() const
Return the upper bound for digital points in this space.
double norm(const NormType type=L_2) const
Implements differential operators on polygonal surfaces from .
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
LinAlg::Triplet Triplet
Type of sparse matrix triplet.
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
DigitalPlane::Point Vector
SMesh::Vertices Vertices
DGtal is the top-level namespace which contains all DGtal functions and types.
Aim: Represents an embedded mesh as faces and a list of vertices. Vertices may be shared among faces ...
Definition: SurfaceMesh.h:92
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:647
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
std::vector< Vertices > computeManifoldBoundaryChains() const
Definition: SurfaceMesh.h:483
Scalar distance(const Vertex i, const Vertex j) const
Definition: SurfaceMesh.h:702
Edges computeManifoldBoundaryEdges() const
Size nbVertices() const
Definition: SurfaceMesh.h:288
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
EigenLinearAlgebraBackend::DenseMatrix DenseMatrix
ShortcutsGeometry< Z3i::KSpace > SHG3
TriMesh::Face Face
PointVector< 3, double > RealPoint
TriMesh::Vertex Vertex