DGtal  1.4.beta
testPolygonalCalculus.cpp
Go to the documentation of this file.
1 
32 #include <iostream>
34 #include "DGtal/base/Common.h"
35 #include "ConfigTest.h"
36 #include "DGtalCatch.h"
37 #include "DGtal/helpers/StdDefs.h"
38 
39 #include "DGtal/dec/PolygonalCalculus.h"
40 #include "DGtal/shapes/SurfaceMesh.h"
41 #include "DGtal/shapes/MeshHelpers.h"
42 #include "DGtal/helpers/Shortcuts.h"
43 #include "DGtal/math/linalg/DirichletConditions.h"
45 
46 using namespace std;
47 using namespace DGtal;
48 using namespace Z3i;
49 
51 // Functions for testing class PolygonalCalculus.
53 
55 {
56  return RealPoint(v(0),v(1),v(2));
57 }
58 
59 TEST_CASE( "Testing PolygonalCalculus" )
60 {
62  std::vector<RealPoint> positions = { RealPoint( 0, 0, 0 ) ,
63  RealPoint( 1, 0, 0 ) ,
64  RealPoint( 0, 1, 0 ) ,
65  RealPoint( 1, 1, 0 ) ,
66  RealPoint( 0, 0, 1 ) ,
67  RealPoint( 1, 0, 1 ) ,
68  RealPoint( 0, 1, 1 ) ,
69  RealPoint( 1, 1, 1 ) ,
70  RealPoint( 1, 0, 2 ) ,
71  RealPoint( 0, 0, 2 ) };
72  std::vector<Mesh::Vertices> faces = { { 1, 0, 2, 3 },
73  { 0, 1, 5, 4 } ,
74  { 1, 3, 7, 5 } ,
75  { 3, 2, 6, 7 } ,
76  { 2, 0, 4, 6 } ,
77  { 4, 5, 8, 9 } };
78 
79  Mesh box(positions.cbegin(), positions.cend(),
80  faces.cbegin(), faces.cend());
81 
83 
84  SECTION("Construction and basic operators")
85  {
86  REQUIRE( boxCalculus.isValid() );
87  REQUIRE( boxCalculus.nbVertices() == positions.size() );
88 
90  auto x = boxCalculus.X(f);
91  auto d = boxCalculus.D(f);
92  auto a = boxCalculus.A(f);
93 
94  //Checking X
96  REQUIRE( vecToRealPoint(vec ) == positions[1]);
97  vec = x.row(1);
98  REQUIRE( vecToRealPoint(vec ) == positions[0]);
99  vec = x.row(2);
100  REQUIRE( vecToRealPoint(vec ) == positions[2]);
101  vec = x.row(3);
102  REQUIRE( vecToRealPoint(vec ) == positions[3]);
103 
104  trace.info()<< boxCalculus <<std::endl;
105 
106  //Some D and A
107  REQUIRE( d(1,1) == -1 );
108  REQUIRE( d(0,1) == 1 );
109  REQUIRE( d(0,2) == 0 );
110 
111  REQUIRE( a(1,1) == 0.5 );
112  REQUIRE( a(0,1) == 0.5 );
113  REQUIRE( a(0,2) == 0 );
114 
115  auto vectorArea = boxCalculus.vectorArea(f);
116 
117  //Without correction, this should match
118  for(auto f=0; f<6; ++f)
119  REQUIRE( boxCalculus.faceArea(f) == box.faceArea(f) );
120 
121  box.computeFaceNormalsFromPositions();
122  for(auto f=0; f<6; ++f)
123  {
124  auto cn = boxCalculus.faceNormalAsDGtalVector(f);
125  auto n = box.faceNormal(f);
126  REQUIRE( cn == n );
127  }
128 
129  RealPoint c = boxCalculus.centroidAsDGtalPoint(f);
130  RealPoint expected(0.5,0.5,0.0);
131  REQUIRE(c == expected);
132  }
133 
134  SECTION("Derivatives")
135  {
137  auto d = boxCalculus.D(f);
138 
139  auto nf = boxCalculus.faceDegree(f);
141  phi << 1.0, 3.0, 2.0, 6.0;
142  expected << 2,-1,4,-5;
143  auto dphi = d*phi; // n_f x 1 matrix
144  REQUIRE(dphi == expected);
145 
146  }
147 
148  SECTION("Structural propertes")
149  {
151  auto nf = boxCalculus.faceDegree(f);
153  phi << 1.0, 3.0, 2.0, 6.0;
154 
155  auto G = boxCalculus.gradient(f);
156  auto gphi = G*phi;
157  auto coG = boxCalculus.coGradient(f);
158  auto cogphi = coG*phi;
159 
160  // grad . cograd == 0
161  REQUIRE( gphi.dot(cogphi) == 0.0);
162 
163  // Gf = Uf Df
164  REQUIRE( G == boxCalculus.sharp(f)*boxCalculus.D(f));
165 
166  // UV = I - nn^t (lemma4)
168  REQUIRE( boxCalculus.sharp(f)*boxCalculus.flat(f) == PolygonalCalculus< RealPoint,RealVector >::DenseMatrix::Identity(3,3) - n*n.transpose() );
169 
170  // P^2 = P (lemma6)
171  auto P = boxCalculus.P(f);
172  REQUIRE( P*P == P);
173 
174  // PV=0 (lemma5)
175  REQUIRE( (P*boxCalculus.flat(f)).norm() == 0.0);
176  }
177 
178  SECTION("Div / Curl")
179  {
181  auto curl = boxCalculus.curl(f);
182  //Not a great test BTW
183  REQUIRE(curl.norm() == 2.0);
184  }
185 
186  SECTION("Local Laplace-Beltrami")
187  {
189  auto nf = box.incidentVertices(f).size();
190 
191  auto L = boxCalculus.laplaceBeltrami(f);
193  phi << 1.0, 1.0, 1.0, 1.0;
194  expected << 0,0,0,0;
195  auto lphi = L*phi;
196  REQUIRE( lphi == expected);
197  }
198 
199  SECTION("Local Connection-Laplace-Beltrami")
200  {
202  auto nf = box.incidentVertices(f).size();
203 
204  auto L = boxCalculus.connectionLaplacian(f);
206  phi << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0;
207  //since connection laplacian transports the phi vectors to the face,
208  //it's not expected to have 0 since these vectors aren't actually the same
209  auto lphi = L*phi;
210  //but we can still check that L is semi-definite
211  double det = L.determinant()+1.;
212  REQUIRE( det == Approx(1.0));
213  REQUIRE( lphi[2] == Approx(-3.683));
214  }
215 
216  SECTION("Covariant Operators")
217  {
219  auto nf = box.incidentVertices(f).size();
220 
222  phi << 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0;
223  auto CG = boxCalculus.covariantGradient(f,phi);
224  auto CP = boxCalculus.covariantProjection(f,phi);
225 
226  //check sizes
227  REQUIRE( CG.rows() == 2);
228  REQUIRE( CG.cols() == 2);
229  REQUIRE( CP.rows() == faces[f].size());
230  REQUIRE( CP.cols() == 2);
231 
232  REQUIRE( CG(0,0) == Approx(0.707106));
233  REQUIRE( CP(0,0) == Approx(1.224744));
234  }
235  SECTION("Check lumped mass matrix")
236  {
238  double a=0.0;
239  for(auto v=0; v < box.nbVertices(); ++v )
240  a += M.coeffRef(v,v);
241 
242  double fa=0.0;
243  for(auto f=0; f < box.nbFaces(); ++f )
244  fa += box.faceArea(f);
245  REQUIRE( a == fa );
246  }
247 
248  SECTION("Checking cache")
249  {
250  auto cacheU = boxCalculus.getOperatorCacheMatrix( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.sharp(f);} );
251  REQUIRE( cacheU.size() == 6 );
252 
253  auto cacheC = boxCalculus.getOperatorCacheVector( [&](const PolygonalCalculus< RealPoint,RealVector >::Face f){ return boxCalculus.centroid(f);} );
254  REQUIRE( cacheC.size() == 6 );
255  }
256 
257  SECTION("Internal cache")
258  {
259  PolygonalCalculus< RealPoint,RealVector > boxCalculusCached(box,true);
260  trace.info()<< boxCalculusCached <<std::endl;
261 
262  trace.beginBlock("Without cache");
263  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix L(box.nbVertices(),box.nbVertices());
264  for(auto i=0; i < 1000 ; ++i)
265  L += i*boxCalculus.globalLaplaceBeltrami();
266  auto tps = trace.endBlock();
267 
268  trace.beginBlock("With cache");
269  PolygonalCalculus< RealPoint,RealVector >::SparseMatrix LC(box.nbVertices(),box.nbVertices());
270  for(auto i=0; i < 1000 ; ++i)
271  LC += i*boxCalculusCached.globalLaplaceBeltrami();
272  auto tpsC = trace.endBlock();
273  REQUIRE(tpsC < tps);
274  REQUIRE(L.norm() == Approx(LC.norm()));
275 
276  }
277 
278 }
279 
280 TEST_CASE( "Testing PolygonalCalculus and DirichletConditions" )
281 {
282  typedef Shortcuts< KSpace > SH3;
284  typedef Mesh::Index Index;
287 
288  // Build a more complex surface.
289  auto params = SH3::defaultParameters();
290 
291  params( "polynomial", "0.1*y*y -0.1*x*x - 2.0*z" )( "gridstep", 2.0 );
292  auto implicit_shape = SH3::makeImplicitShape3D ( params );
293  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
294  auto K = SH3::getKSpace( params );
295  auto binary_image = SH3::makeBinaryImage( digitized_shape, params );
296  auto surface = SH3::makeDigitalSurface( binary_image, K, params );
297  auto primalSurface = SH3::makePrimalSurfaceMesh(surface);
298 
299  std::vector<std::vector< Index > > faces;
300  std::vector<RealPoint> positions = primalSurface->positions();
301  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
302  faces.push_back(primalSurface->incidentVertices( face ));
303 
304  Mesh surfmesh = Mesh( positions.begin(), positions.end(),
305  faces.begin(), faces.end() );
306  auto boundaryEdges = surfmesh.computeManifoldBoundaryEdges();
307 
308  SECTION("Check surface")
309  {
310  REQUIRE( surfmesh.nbVertices() == 1364 );
311  REQUIRE( surfmesh.nbFaces() == 1279 );
312  REQUIRE( boundaryEdges.size() == 168 );
313  }
314 
315  // Builds calculus and solve a Poisson problem with Dirichlet boundary conditions
316  PolyDEC calculus( surfmesh );
317  // Laplace opeartor
319  // value on boundary
320  PolyDEC::Vector g = calculus.form0();
321  // characteristic set of boundary
322  DC::IntegerVector b = DC::IntegerVector::Zero( g.rows() );
323 
324  SECTION("Solve Poisson problem with boundary Dirichlet conditions")
325  {
326  for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
327  {
328  std::cout << "scale=" << scale << std::endl;
329  auto phi = [&]( Index v)
330  {
331  return cos(scale*(surfmesh.position(v)[0]))
332  * (scale*surfmesh.position(v)[1]);
333  };
334 
335  for(auto &e: boundaryEdges)
336  {
337  auto adjVertices = surfmesh.edgeVertices(e);
338  auto v1 = adjVertices.first;
339  auto v2 = adjVertices.second;
340  g(v1) = phi(v1);
341  g(v2) = phi(v2);
342  b(v1) = 1;
343  b(v2) = 1;
344  }
345  // Solve Δu=0 with g as boundary conditions
346  PolyDEC::Solver solver;
347  PolyDEC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, b );
348  solver.compute( L_dirichlet );
349  REQUIRE( solver.info() == Eigen::Success );
350  PolyDEC::Vector g_dirichlet = DC::dirichletVector( L, g, b, g );
351  PolyDEC::Vector x_dirichlet = solver.solve( g_dirichlet );
352  REQUIRE( solver.info() == Eigen::Success );
353  PolyDEC::Vector u = DC::dirichletSolution( x_dirichlet, b, g );
354  double min_phi = 0.0;
355  double max_phi = 0.0;
356  double min_u = 0.0;
357  double max_u = 0.0;
358  double min_i_u = 0.0;
359  double max_i_u = 0.0;
360  double exp_u = 0.0;
361  double exp_u2 = 0.0;
362  for ( auto v = 0; v < surfmesh.nbVertices(); ++v )
363  {
364  min_phi = std::min( min_phi, phi( v ) );
365  max_phi = std::max( max_phi, phi( v ) );
366  min_u = std::min( min_u , u ( v ) );
367  max_u = std::max( max_u , u ( v ) );
368  if ( b( v ) == 0.0 )
369  {
370  min_i_u = std::min( min_i_u, u ( v ) );
371  max_i_u = std::max( max_i_u, u ( v ) );
372  }
373  }
374  REQUIRE( min_phi <= min_u );
375  REQUIRE( max_phi >= max_u );
376  REQUIRE( min_phi < min_i_u );
377  REQUIRE( max_phi > max_i_u );
378  } // for ( double scale = 0.1; scale < 2.0; scale *= 2.0 )
379  }
380 };
381 
DGtal::PolygonalCalculus::Face
MySurfaceMesh::Face Face
Face type.
Definition: PolygonalCalculus.h:87
vecToRealPoint
RealPoint vecToRealPoint(PolygonalCalculus< RealPoint, RealPoint >::Vector &v)
Definition: testPolygonalCalculus.cpp:54
DGtal::PolygonalCalculus::getOperatorCacheMatrix
std::vector< DenseMatrix > getOperatorCacheMatrix(const std::function< DenseMatrix(Face)> &perFaceOperator) const
Definition: PolygonalCalculus.h:942
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
max
int max(int a, int b)
Definition: testArithmeticalDSS.cpp:1108
DGtal::PolygonalCalculus::globalLaplaceBeltrami
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Definition: PolygonalCalculus.h:798
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:102
DGtal::PolygonalCalculus::centroidAsDGtalPoint
Real3dPoint centroidAsDGtalPoint(const Face f) const
Definition: PolygonalCalculus.h:436
Index
SMesh::Index Index
Definition: fullConvexitySphereGeodesics.cpp:117
DGtal::Mesh::Index
VertexStorage::size_type Index
Definition: Mesh.h:120
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::SurfaceMesh::computeManifoldBoundaryEdges
Edges computeManifoldBoundaryEdges() const
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
REQUIRE
REQUIRE(domain.isInside(aPoint))
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:70
DGtal::PolygonalCalculus::covariantProjection
DenseMatrix covariantProjection(const Face f, const Vector &uf)
Definition: PolygonalCalculus.h:715
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
DGtal::PolygonalCalculus::form0
Form form0() const
Definition: PolygonalCalculus.h:771
DGtal::PolygonalCalculus::laplaceBeltrami
DenseMatrix laplaceBeltrami(const Face f, const double lambda=1.0) const
Definition: PolygonalCalculus.h:545
DGtal::SurfaceMesh::nbVertices
Size nbVertices() const
Definition: SurfaceMesh.h:280
SparseMatrix
EigenLinearAlgebraBackend::SparseMatrix SparseMatrix
Definition: testHeatLaplace.cpp:50
DGtal::PolygonalCalculus::P
DenseMatrix P(const Face f) const
Definition: PolygonalCalculus.h:461
DGtal::Trace::info
std::ostream & info()
DGtal::PolygonalCalculus::D
DenseMatrix D(const Face f) const
Definition: PolygonalCalculus.h:261
DGtal::Shortcuts::defaultParameters
static Parameters defaultParameters()
Definition: Shortcuts.h:203
DGtal::PolygonalCalculus::sharp
DenseMatrix sharp(const Face f) const
Definition: PolygonalCalculus.h:445
DGtal::Mesh
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Definition: Mesh.h:91
DGtal::PolygonalCalculus::getOperatorCacheVector
std::vector< Vector > getOperatorCacheVector(const std::function< Vector(Face)> &perFaceVectorOperator) const
Definition: PolygonalCalculus.h:966
DGtal::PolygonalCalculus::gradient
DenseMatrix gradient(const Face f) const
Definition: PolygonalCalculus.h:390
DGtal::Shortcuts::makeImplicitShape3D
static CountedPtr< ImplicitShape3D > makeImplicitShape3D(const Parameters &params=parametersImplicitShape3D())
Definition: Shortcuts.h:282
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::PolygonalCalculus::faceDegree
size_t faceDegree(Face f) const
Definition: PolygonalCalculus.h:1006
DGtal::PolygonalCalculus::centroid
Vector centroid(const Face f) const
Definition: PolygonalCalculus.h:428
DGtal::PolygonalCalculus::covariantGradient
DenseMatrix covariantGradient(const Face f, const Vector &uf)
Definition: PolygonalCalculus.h:699
DGtal::SurfaceMesh::edgeVertices
const VertexPair & edgeVertices(Edge e) const
Definition: SurfaceMesh.h:329
DGtal::PolygonalCalculus::isValid
bool isValid() const
Definition: PolygonalCalculus.h:1054
DGtal::Mesh::nbFaces
Size nbFaces() const
DGtal::PolygonalCalculus::curl
DenseMatrix curl(const Face f) const
Definition: PolygonalCalculus.h:519
DGtal::Shortcuts::makeDigitizedImplicitShape3D
static CountedPtr< DigitizedImplicitShape3D > makeDigitizedImplicitShape3D(CountedPtr< ImplicitShape3D > shape, Parameters params=parametersDigitizedImplicitShape3D())
Definition: Shortcuts.h:523
Vector
FreemanChain< int >::Vector Vector
Definition: testCombinDSS.cpp:60
DGtal::PolygonalCalculus::nbVertices
size_t nbVertices() const
Definition: PolygonalCalculus.h:1012
TEST_CASE
TEST_CASE("Testing PolygonalCalculus")
Definition: testPolygonalCalculus.cpp:59
DGtal::PolygonalCalculus::X
DenseMatrix X(const Face f) const
Definition: PolygonalCalculus.h:236
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::PointVector
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:165
SECTION
SECTION("Testing constant forward iterators")
Definition: testSimpleRandomAccessRangeFromPoint.cpp:66
DGtal::PolygonalCalculus::faceNormal
Vector faceNormal(const Face f) const
Definition: PolygonalCalculus.h:348
DGtal::PolygonalCalculus::flat
DenseMatrix flat(const Face f) const
Definition: PolygonalCalculus.h:404
DGtal::PolygonalCalculus::globalLumpedMassMatrix
SparseMatrix globalLumpedMassMatrix() const
Definition: PolygonalCalculus.h:825
DGtal::PolygonalCalculus::connectionLaplacian
DenseMatrix connectionLaplacian(const Face &f, double lambda=1.0) const
Definition: PolygonalCalculus.h:753
DGtal::Shortcuts::makeBinaryImage
static CountedPtr< BinaryImage > makeBinaryImage(Domain shapeDomain)
Definition: Shortcuts.h:561
RealPoint
Z2i::RealPoint RealPoint
Definition: testAstroid2D.cpp:46
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288
DGtal::PolygonalCalculus::A
DenseMatrix A(const Face f) const
Definition: PolygonalCalculus.h:295