DGtal 1.3.0
Searching...
No Matches
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#include "ConfigExamples.h"
41
42#include <DGtal/io/boards/Board2D.h>
44
45using namespace DGtal;
46
48// Using standard 3D digital space.
51// The following typedefs are useful
55typedef SurfMesh::Face Face;
60typedef PC::Solver Solver;
61typedef PC::Vector Vector;
62typedef PC::Triplet Triplet;
63typedef std::vector<Vertex> chain;
65typedef Conditions::IntegerVector IntegerVector;
66
67
68//Polyscope global
69polyscope::SurfaceMesh *psMesh;
70polyscope::SurfaceMesh *psParam;
71
72//DEC
73PC *calculus;
74
75//surface mesh global
76SurfMesh surfmesh;
77std::vector<std::vector<size_t>> faces;
78std::vector<RealPoint> positions;
79
80
87std::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
114void 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
130 for (size_t v = 0;v<n;v++)
131 {
132 pos[v] = RealPoint{scale*UV(v,0),scale*UV(v,1),0.} + avg;
133 }
134 polyscope::registerSurfaceMesh("On circle parametrization", pos, faces)->setEnabled(false);
135
136 Board2D board;
137 auto pairs = surfmesh.allEdgeVertices();
138 for(auto p: pairs)
139 {
140 if ((std::pow(UV(p.first,0)-UV(p.second,0),2.0)+std::pow(UV(p.first,1)- UV(p.second,1),2.0)) > 0.0)
141 board.drawLine(UV(p.first,0), UV(p.first,1), UV(p.second,0), UV(p.second,1));
142 }
143 board.saveSVG("layout.svg");
144
145}
146
153DenseMatrix HarmonicParametrization()
154{
155 auto n = surfmesh.nbVertices();
156 std::cout<<"Nb boundary edges = "<< surfmesh.computeManifoldBoundaryEdges().size()<<std::endl;
157 std::vector<chain> chains = surfmesh.computeManifoldBoundaryChains();
158 //choose longest chain as boundary of the parametrization
159 std::cout<<"Nb boundaries = "<< chains.size() << std::endl;
160
161 //choose longest chain as boundary of the parametrization
162 auto B = *std::max_element(chains.begin(),chains.end(),[] (const chain& A,const chain& B) {return A.size() < B.size();});
163
164 IntegerVector boundary = IntegerVector::Zero(n);
165 for (Vertex v : B)
166 boundary(v) = 1;
167
168 std::pair<Vector,Vector> uv_b = FixBoundaryParametrization(B);//maps boundary to circle
169
170 calculus = new PC(surfmesh);
171
172 //Impose dirichlet boundary condition to laplace problem
173 Vector Z = Vector::Zero(n);
174 SparseMatrix L = calculus->globalLaplaceBeltrami();
175 SparseMatrix L_d = Conditions::dirichletOperator( L, boundary );
176
177 PC::Solver solver;
178 solver.compute(L_d);
179
180 Vector b_u = Conditions::dirichletVector( L, Z,boundary, uv_b.first );
181 Vector b_v = Conditions::dirichletVector( L, Z,boundary, uv_b.second );
182
183 Vector rslt_u_d = solver.solve(b_u);
184 Vector rslt_v_d = solver.solve(b_v);
185
186 Vector rslt_u = Conditions::dirichletSolution(rslt_u_d,boundary,uv_b.first);
187 Vector rslt_v = Conditions::dirichletSolution(rslt_v_d,boundary,uv_b.second);
188
189 DenseMatrix uv(n,2);
190 uv.col(0) = rslt_u;
191 uv.col(1) = rslt_v;
192
193 return uv;
194}
195
196int main(int argc, char **argv)
197{
198 //Import Voxel Model
200 params("surfaceComponents", "All");
201
203 std::string inputFilename(examplesPath + "samples/lucy-64.vol" );
204
205 auto binary_image = SH3::makeBinaryImage(inputFilename, params );
206
207 //offset K space to create boundary
208 auto K = SH3::getKSpace( binary_image, params );
209 K.init(K.lowerBound()+SH3::Point(5,0,0),K.upperBound(),true);
210
211 auto surface = SH3::makeDigitalSurface( binary_image, K, params );
212 auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
213
214
215 //Need to convert the faces
216 for(size_t face= 0 ; face < primalSurface->nbFaces(); ++face)
217 faces.push_back(primalSurface->incidentVertices( face ));
218
219 //Recasting to vector of vertices
220 positions = primalSurface->positions();
221
222 surfmesh = SurfMesh(positions.begin(),
223 positions.end(),
224 faces.begin(),
225 faces.end());
226
227
228 // Initialize polyscope
229 polyscope::init();
230
231 //compute parametrization
232 DenseMatrix UV = HarmonicParametrization();
233
234 //visualize parametrization on 2D circle
235 VisualizeParametrizationOnCircle(UV);
236
237 psMesh = polyscope::registerSurfaceMesh("Digital Surface", positions, faces);
239
240 //set correct view for voxel mesh
241 polyscope::view::upDir = polyscope::view::UpDir::XUp;
242
243 polyscope::show();
244 return EXIT_SUCCESS;
245}
Aim: This class specializes a 'Board' class so as to display DGtal objects more naturally (with <<)....
Definition: Board2D.h:71
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)
const Point & lowerBound() const
Return the lower bound for digital points in this space.
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 & upperBound() const
Return the upper bound for digital points in this space.
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
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...
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
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 drawLine(double x1, double y1, double x2, double y2, int depthValue=-1)
Definition: Board.cpp:368
void saveSVG(const char *filename, PageSize size=Board::BoundingBox, double margin=10.0) const
Definition: Board.cpp:1012
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
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
Scalar distance(const Vertex i, const Vertex j) const
Definition: SurfaceMesh.h:692
Edges computeManifoldBoundaryEdges() const
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
Size nbVertices() const
Definition: SurfaceMesh.h:280
const std::vector< VertexPair > & allEdgeVertices() const
Definition: SurfaceMesh.h:381
std::vector< Vertices > computeManifoldBoundaryChains() const
Definition: SurfaceMesh.h:473
int main()
Definition: testBits.cpp:56
FreemanChain< int >::Vector Vector
KSpace K
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
EigenLinearAlgebraBackend::DenseMatrix DenseMatrix
TriMesh::Face Face
TriMesh::Vertex Vertex