DGtal  1.4.beta
PolygonalCalculus.h
1 
17 #pragma once
18 
30 // Inclusions
32 #include <iostream>
33 #include <functional>
34 #include <vector>
35 #include <string>
36 #include <map>
37 #include <unordered_map>
38 #include "DGtal/base/ConstAlias.h"
39 #include "DGtal/base/Common.h"
40 #include "DGtal/shapes/SurfaceMesh.h"
41 #include "DGtal/math/linalg/EigenSupport.h"
43 
44 namespace DGtal
45 {
46 
48 // template class PolygonalCalculus
69 template <typename TRealPoint, typename TRealVector>
71 {
72  // ----------------------- Standard services ------------------------------
73 public:
74 
76  static const Dimension dimension = TRealPoint::dimension;
77  BOOST_STATIC_ASSERT( ( dimension == 3 ) );
78 
81 
85  typedef typename MySurfaceMesh::Vertex Vertex;
87  typedef typename MySurfaceMesh::Face Face;
92 
98  typedef Vector Form;
105 
108 
111 
117  bool globalInternalCacheEnabled = false):
118  mySurfaceMesh(&surf), myGlobalCacheEnabled(globalInternalCacheEnabled)
119  {
120  myEmbedder =[&](Face f,Vertex v){ return mySurfaceMesh->position(v);};
122  init();
123  };
124 
132  const std::function<Real3dPoint(Face,Vertex)> &embedder,
133  bool globalInternalCacheEnabled = false):
134  mySurfaceMesh(&surf), myEmbedder(embedder), myGlobalCacheEnabled(globalInternalCacheEnabled)
135  {
137  init();
138  };
139 
147  const std::function<Real3dVector(Vertex)> &embedder,
148  bool globalInternalCacheEnabled = false):
149  mySurfaceMesh(&surf), myVertexNormalEmbedder(embedder),
150  myGlobalCacheEnabled(globalInternalCacheEnabled)
151  {
152  myEmbedder = [&](Face f,Vertex v){ return mySurfaceMesh->position(v); };
153  init();
154  };
155 
167  const std::function<Real3dPoint(Face,Vertex)> &pos_embedder,
168  const std::function<Vector(Vertex)> &normal_embedder,
169  bool globalInternalCacheEnabled = false) :
170  mySurfaceMesh(&surf), myEmbedder(pos_embedder),
171  myVertexNormalEmbedder(normal_embedder),
172  myGlobalCacheEnabled(globalInternalCacheEnabled)
173  {
174  init();
175  };
176 
180  PolygonalCalculus() = delete;
181 
185  ~PolygonalCalculus() = default;
186 
191  PolygonalCalculus ( const PolygonalCalculus & other ) = delete;
192 
197  PolygonalCalculus ( PolygonalCalculus && other ) = delete;
198 
204  PolygonalCalculus & operator= ( const PolygonalCalculus & other ) = delete;
205 
211  PolygonalCalculus & operator= ( PolygonalCalculus && other ) = delete;
212 
214 
215  // ----------------------- embedding services --------------------------
216  //MARK: Embedding services
219 
222  void setEmbedder(const std::function<Real3dPoint(Face,Vertex)> &externalFunctor)
223  {
224  myEmbedder = externalFunctor;
225  }
226 
227 
228  // ----------------------- Per face operators --------------------------------------
229  //MARK: Per face operator on scalars
232 
236  DenseMatrix X(const Face f) const
237  {
238  if (checkCache(X_,f))
239  return myGlobalCache[X_][f];
240 
241  const auto vertices = mySurfaceMesh->incidentVertices(f);
242  const auto nf = myFaceDegree[f];
243  DenseMatrix Xt(nf,3);
244  size_t cpt=0;
245  for(auto v: vertices)
246  {
247  Xt(cpt,0) = myEmbedder(f,v)[0];
248  Xt(cpt,1) = myEmbedder(f,v)[1];
249  Xt(cpt,2) = myEmbedder(f,v)[2];
250  ++cpt;
251  }
252 
253  setInCache(X_,f,Xt);
254  return Xt;
255  }
256 
257 
261  DenseMatrix D(const Face f) const
262  {
263  if (checkCache(D_,f))
264  return myGlobalCache[D_][f];
265 
266  const auto nf = myFaceDegree[f];
267  DenseMatrix d = DenseMatrix::Zero(nf ,nf);
268  for(auto i=0; i < nf; ++i)
269  {
270  d(i,i) = -1.;
271  d(i, (i+1)%nf) = 1.;
272  }
273 
274  setInCache(D_,f,d);
275  return d;
276  }
277 
281  DenseMatrix E(const Face f) const
282  {
283  if (checkCache(E_,f))
284  return myGlobalCache[E_][f];
285 
286  DenseMatrix op = D(f)*X(f);
287 
288  setInCache(E_,f,op);
289  return op;
290  }
291 
295  DenseMatrix A(const Face f) const
296  {
297  if (checkCache(A_,f))
298  return myGlobalCache[A_][f];
299 
300  const auto nf = myFaceDegree[f];
301  DenseMatrix a = DenseMatrix::Zero(nf ,nf);
302  for(auto i=0; i < nf; ++i)
303  {
304  a(i, (i+1)%nf) = 0.5;
305  a(i,i) = 0.5;
306  }
307 
308  setInCache(A_,f,a);
309  return a;
310  }
311 
312 
316  Vector vectorArea(const Face f) const
317  {
318  Real3dPoint af(0.0,0.0,0.0);
319  const auto vertices = mySurfaceMesh->incidentVertices(f);
320  auto it = vertices.cbegin();
321  auto itnext = vertices.cbegin();
322  ++itnext;
323  while (it != vertices.cend())
324  {
325  auto xi = myEmbedder(f,*it);
326  auto xip = myEmbedder(f,*itnext);
327  af += xi.crossProduct(xip);
328  ++it;
329  ++itnext;
330  if (itnext == vertices.cend())
331  itnext = vertices.cbegin();
332  }
333  Eigen::Vector3d output = {af[0],af[1],af[2]};
334  return 0.5*output;
335  }
336 
340  double faceArea(const Face f) const
341  {
342  return vectorArea(f).norm();
343  }
344 
348  Vector faceNormal(const Face f) const
349  {
350  Vector v = vectorArea(f);
351  v.normalize();
352  return v;
353  }
354 
359  {
360  Vector v = faceNormal(f);
361  return {v(0),v(1),v(2)};
362  }
363 
367  DenseMatrix coGradient(const Face f) const
368  {
369  if (checkCache(COGRAD_,f))
370  return myGlobalCache[COGRAD_][f];
371  DenseMatrix op = E(f).transpose() * A(f);
372  setInCache(COGRAD_, f, op);
373  return op;
374  }
375 
378  DenseMatrix bracket(const Vector &n) const
379  {
380  DenseMatrix brack(3,3);
381  brack << 0.0 , -n(2), n(1),
382  n(2), 0.0 , -n(0),
383  -n(1) , n(0),0.0 ;
384  return brack;
385  }
386 
390  DenseMatrix gradient(const Face f) const
391  {
392  if (checkCache(GRAD_,f))
393  return myGlobalCache[GRAD_][f];
394 
395  DenseMatrix op = -1.0/faceArea(f) * bracket( faceNormal(f) ) * coGradient(f);
396 
397  setInCache(GRAD_,f,op);
398  return op;
399  }
400 
404  DenseMatrix flat(const Face f) const
405  {
406  if (checkCache(FLAT_,f))
407  return myGlobalCache[FLAT_][f];
408  DenseMatrix n = faceNormal(f);
409  DenseMatrix op = E(f)*( DenseMatrix::Identity(3,3) - n*n.transpose());
410  setInCache(FLAT_,f,op);
411  return op;
412  }
413 
417  DenseMatrix B(const Face f) const
418  {
419  if (checkCache(B_,f))
420  return myGlobalCache[B_][f];
421  DenseMatrix res = A(f) * X(f);
422  setInCache(B_,f,res);
423  return res;
424  }
425 
428  Vector centroid(const Face f) const
429  {
430  const auto nf = myFaceDegree[f];
431  return 1.0/(double)nf * X(f).transpose() * Vector::Ones(nf);
432  }
433 
437  {
438  const Vector c = centroid(f);
439  return {c(0),c(1),c(2)};
440  }
441 
445  DenseMatrix sharp(const Face f) const
446  {
447  if (checkCache(SHARP_,f))
448  return myGlobalCache[SHARP_][f];
449 
450  const auto nf = myFaceDegree[f];
451  DenseMatrix op = 1.0/faceArea(f) * bracket(faceNormal(f)) *
452  ( B(f).transpose() - centroid(f)* Vector::Ones(nf).transpose() );
453 
454  setInCache(SHARP_,f,op);
455  return op;
456  }
457 
461  DenseMatrix P(const Face f) const
462  {
463  if (checkCache(P_,f))
464  return myGlobalCache[P_][f];
465 
466  const auto nf = myFaceDegree[f];
467  DenseMatrix op = DenseMatrix::Identity(nf,nf) - flat(f)*sharp(f);
468 
469  setInCache(P_, f, op);
470  return op;
471  }
472 
477  DenseMatrix M(const Face f, const double lambda=1.0) const
478  {
479  if (checkCache(M_,f))
480  return myGlobalCache[M_][f];
481 
482  DenseMatrix Uf = sharp(f);
483  DenseMatrix Pf = P(f);
484  DenseMatrix op = faceArea(f) * Uf.transpose()*Uf + lambda * Pf.transpose()*Pf;
485 
486  setInCache(M_,f,op);
487  return op;
488  }
489 
505  DenseMatrix divergence(const Face f, const double lambda=1.0) const
506  {
507  if (checkCache(DIVERGENCE_,f))
508  return myGlobalCache[DIVERGENCE_][f];
509 
510  DenseMatrix op = -1.0 * D(f).transpose() * M(f);
511  setInCache(DIVERGENCE_,f,op);
512 
513  return op;
514  }
515 
519  DenseMatrix curl(const Face f) const
520  {
521  if (checkCache(CURL_,f))
522  return myGlobalCache[CURL_][f];
523 
524  DenseMatrix op = DenseMatrix::Identity(myFaceDegree[f],myFaceDegree[f]);
525 
526  setInCache(CURL_,f,op);
527  return op;
528  }
529 
530 
545  DenseMatrix laplaceBeltrami(const Face f, const double lambda=1.0) const
546  {
547  if (checkCache(L_,f))
548  return myGlobalCache[L_][f];
549 
550  DenseMatrix Df = D(f);
551  // Laplacian is a negative operator.
552  DenseMatrix op = -1.0 * Df.transpose() * M(f,lambda) * Df;
553 
554  setInCache(L_, f, op);
555  return op;
556  }
558 
559  // ----------------------- Vector calculus----------------------------------
560  //MARK: Vector Field Calculus
563 
564 public:
567  DenseMatrix Tv(const Vertex & v) const
568  {
569  Eigen::Vector3d nv = n_v(v);
570  ASSERT(std::abs(nv.norm() - 1.0) < 0.001);
571  const auto & N = getSurfaceMeshPtr()->neighborVertices(v);
572  auto neighbor = *N.begin();
573  Real3dPoint tangentVector = getSurfaceMeshPtr()->position(v) -
574  getSurfaceMeshPtr()->position(neighbor);
575  Eigen::Vector3d w = toVec3(tangentVector);
576  Eigen::Vector3d uu = project(w,nv).normalized();
577  Eigen::Vector3d vv = nv.cross(uu);
578 
579  DenseMatrix tanB(3,2);
580  tanB.col(0) = uu;
581  tanB.col(1) = vv;
582  return tanB;
583  }
584 
587  DenseMatrix Tf(const Face & f) const
588  {
589  Eigen::Vector3d nf = faceNormal(f);
590  ASSERT(std::abs(nf.norm() - 1.0) < 0.001);
591  const auto & N = getSurfaceMeshPtr()->incidentVertices(f);
592  auto v1 = *(N.begin());
593  auto v2 = *(N.begin() + 1);
594  Real3dPoint tangentVector =
596  Eigen::Vector3d w = toVec3(tangentVector);
597  Eigen::Vector3d uu = project(w,nf).normalized();
598  Eigen::Vector3d vv = nf.cross(uu);
599 
600  DenseMatrix tanB(3,2);
601  tanB.col(0) = uu;
602  tanB.col(1) = vv;
603  return tanB;
604  }
605 
611  Vector toExtrinsicVector(const Vertex v, const Vector & I) const
612  {
613  DenseMatrix T = Tv(v);
614  return T.col(0) * I(0) + T.col(1) * I(1);
615  }
616 
621  std::vector<Vector> toExtrinsicVectors(const std::vector<Vector> & I) const
622  {
623  std::vector<Vector> ext(mySurfaceMesh->nbVertices());
624  for (auto v = 0; v < mySurfaceMesh->nbVertices(); v++)
625  ext[v] = toExtrinsicVector(v,I[v]);
626  return ext;
627  }
628 
631  DenseMatrix Qvf(const Vertex & v, const Face & f) const
632  {
633  Eigen::Vector3d nf = faceNormal(f);
634  Eigen::Vector3d nv = n_v(v);
635  double c = nv.dot(nf);
636  ASSERT(std::abs( c + 1.0) > 0.0001);
637  //Special case for opposite nv and nf vectors.
638  if (std::abs( c + 1.0) < 0.00001)
639  return -Eigen::Matrix3d::Identity();
640 
641  auto vv = nv.cross(nf);
642  DenseMatrix skew = bracket(vv);
643  return Eigen::Matrix3d::Identity() + skew +
644  1.0 / (1.0 + c) * skew * skew;
645  }
646 
649  DenseMatrix Rvf(const Vertex & v, const Face & f) const
650  {
651  return Tf(f).transpose() * Qvf(v,f) * Tv(v);
652  }
653 
657  {
658  DenseMatrix N(myFaceDegree[f],3);
659  uint cpt = 0;
661  {
662  N.block(cpt,0,3,1) = n_v(v).transpose();
663  cpt++;
664  }
665  DenseMatrix GN = gradient(f) * N, Tf = T_f(f);
666 
667  return 0.5 * Tf.transpose() * (GN + GN.transpose()) * Tf;
668  }
669 
678  {
679  DenseMatrix uf_nabla(myFaceDegree[f], 2);
680  size_t cpt = 0;
681  for (auto v : mySurfaceMesh->incidentVertices(f))
682  {
683  uf_nabla.block(cpt,0,1,2) =
684  (Rvf(v,f) * uf.block(2 * cpt,0,2,1)).transpose();
685  ++cpt;
686  }
687  return uf_nabla;
688  }
689 
700  {
701  return Tf(f).transpose() * gradient(f) *
703  }
704 
716  {
717  return P(f) * D(f) * transportAndFormatVectorField(f,uf);
718  }
719 
726  {
727  return kroneckerWithI2(Tf(f).transpose() * gradient(f)) * blockConnection(f);
728  }
729 
736  {
737  return kroneckerWithI2(P(f) * D(f)) * blockConnection(f);
738  ;
739  }
740 
753  DenseMatrix connectionLaplacian(const Face & f, double lambda = 1.0) const
754  {
755  if (checkCache(CON_L_,f))
756  return myGlobalCache[CON_L_][f];
759  DenseMatrix L = -(faceArea(f) * G.transpose() * G + lambda * P.transpose() * P);
760  setInCache(CON_L_,f,L);
761  return L;
762  }
764 
765  // ----------------------- Global operators --------------------------------------
766  //MARK: Global Operators
769 
771  Form form0() const
772  {
773  return Form::Zero( nbVertices() );
774  }
777  {
778  SparseMatrix Id0( nbVertices(), nbVertices() );
779  Id0.setIdentity();
780  return Id0;
781  }
782 
798  SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
799  {
801  std::vector<Triplet> triplets;
802  for( auto f = 0; f < mySurfaceMesh->nbFaces(); ++f )
803  {
804  auto nf = myFaceDegree[f];
805  DenseMatrix Lap = this->laplaceBeltrami(f,lambda);
806  const auto vertices = mySurfaceMesh->incidentVertices(f);
807  for(auto i=0; i < nf; ++i)
808  for(auto j=0; j < nf; ++j)
809  {
810  auto v = Lap(i,j);
811  if (v!= 0.0)
812  triplets.emplace_back( Triplet( (SparseMatrix::StorageIndex)vertices[ i ], (SparseMatrix::StorageIndex)vertices[ j ],
813  Lap( i, j ) ) );
814  }
815  }
816  lapGlobal.setFromTriplets(triplets.begin(), triplets.end());
817  return lapGlobal;
818  }
819 
826  {
828  std::vector<Triplet> triplets;
829  for ( auto v = 0; v < mySurfaceMesh->nbVertices(); ++v )
830  {
831  auto faces = mySurfaceMesh->incidentFaces(v);
832  auto varea = 0.0;
833  for(auto f: faces)
834  varea += faceArea(f) /(double)myFaceDegree[f];
835  triplets.emplace_back(Triplet(v,v,varea));
836  }
837  M.setFromTriplets(triplets.begin(),triplets.end());
838  return M;
839  }
840 
846  {
848  for ( int k = 0; k < iM0.outerSize(); ++k )
849  for ( typename SparseMatrix::InnerIterator it( iM0, k ); it; ++it )
850  it.valueRef() = 1.0 / it.value();
851  return iM0;
852  }
853 
871  SparseMatrix globalConnectionLaplace(const double lambda = 1.0) const
872  {
873  auto nv = mySurfaceMesh->nbVertices();
874  SparseMatrix lapGlobal(2 * nv, 2 * nv);
875  SparseMatrix local(2 * nv, 2 * nv);
876  std::vector<Triplet> triplets;
877  for (auto f = 0; f < mySurfaceMesh->nbFaces(); f++)
878  {
879  auto nf = degree(f);
880  DenseMatrix Lap = connectionLaplacian(f,lambda);
881  const auto vertices = mySurfaceMesh->incidentVertices(f);
882  for (auto i = 0u; i < nf; ++i)
883  for (auto j = 0u; j < nf; ++j)
884  for (short k1 = 0; k1 < 2; k1++)
885  for (short k2 = 0; k2 < 2; k2++)
886  {
887  auto v = Lap(2 * i + k1, 2 * j + k2);
888  if (v != 0.0)
889  triplets.emplace_back(Triplet(2 * vertices[i] + k1,
890  2 * vertices[j] + k2, v));
891  }
892  }
893  lapGlobal.setFromTriplets(triplets.begin(), triplets.end());
894  return lapGlobal;
895  }
896 
904  {
905  auto nv = mySurfaceMesh->nbVertices();
906  SparseMatrix M(2 * nv, 2 * nv);
907  std::vector<Triplet> triplets;
908  for (auto v = 0; v < mySurfaceMesh->nbVertices(); ++v)
909  {
910  auto faces = mySurfaceMesh->incidentFaces(v);
911  auto varea = 0.0;
912  for (auto f : faces)
913  varea += faceArea(f) / (double)myFaceDegree[f];
914  triplets.emplace_back(Triplet(2 * v, 2 * v, varea));
915  triplets.emplace_back(Triplet(2 * v + 1, 2 * v + 1, varea));
916  }
917  M.setFromTriplets(triplets.begin(), triplets.end());
918  return M;
919  }
921 
922  // ----------------------- Cache mechanism --------------------------------------
925 
942  std::vector<DenseMatrix> getOperatorCacheMatrix(const std::function<DenseMatrix(Face)> &perFaceOperator) const
943  {
944  std::vector<DenseMatrix> cache;
945  for(auto f=0; f < mySurfaceMesh->nbFaces(); ++f)
946  cache.push_back(perFaceOperator(f));
947  return cache;
948  }
949 
966  std::vector<Vector> getOperatorCacheVector(const std::function<Vector(Face)> &perFaceVectorOperator) const
967  {
968  std::vector<Vector> cache;
969  for(auto f=0; f < mySurfaceMesh->nbFaces(); ++f)
970  cache.push_back(perFaceVectorOperator(f));
971  return cache;
972  }
973 
977  {
978  myGlobalCacheEnabled = true;
979  }
980 
984  {
985  myGlobalCacheEnabled = false;
986  myGlobalCache.clear();
987  }
988 
990 
991  // ----------------------- Common --------------------------------------
992 public:
995 
998  void init()
999  {
1000  updateFaceDegree();
1001  }
1002 
1006  size_t faceDegree(Face f) const
1007  {
1008  return myFaceDegree[f];
1009  }
1010 
1012  size_t nbVertices() const
1013  {
1014  return mySurfaceMesh->nbVertices();
1015  }
1016 
1018  size_t nbFaces() const
1019  {
1020  return mySurfaceMesh->nbFaces();
1021  }
1022 
1025  size_t degree(const Face f) const
1026  {
1027  return myFaceDegree[f];
1028  }
1029 
1032  {
1033  return mySurfaceMesh;
1034  }
1035 
1040  void selfDisplay ( std::ostream & out ) const
1041  {
1042  out << "[PolygonalCalculus]: ";
1044  out<< "internal cache enabled, ";
1045  else
1046  out<<"internal cache disabled, ";
1047  out <<"SurfaceMesh="<<*mySurfaceMesh;
1048  }
1049 
1054  bool isValid() const
1055  {
1056  return true;
1057  }
1058 
1060 
1061  // ------------------------- Protected Datas ------------------------------
1062  //MARK: Protected
1063 
1064 protected:
1067 
1070 
1073  {
1074  myFaceDegree.resize(mySurfaceMesh->nbFaces());
1075  for(auto f = 0; f < mySurfaceMesh->nbFaces(); ++f)
1076  {
1077  auto vertices = mySurfaceMesh->incidentVertices(f);
1078  auto nf = vertices.size();
1079  myFaceDegree[f] = nf;
1080  }
1081  }
1082 
1087  bool checkCache(OPERATOR key, const Face f) const
1088  {
1090  if (myGlobalCache[key].find(f) != myGlobalCache[key].end())
1091  return true;
1092  return false;
1093  }
1094 
1099  void setInCache(OPERATOR key, const Face f,
1100  const DenseMatrix &ope) const
1101  {
1103  myGlobalCache[key][f] = ope;
1104  }
1105 
1110  static Vector project(const Vector & u, const Vector & n)
1111  {
1112  return u - (u.dot(n) / n.squaredNorm()) * n;
1113  }
1114 
1119  static Vector toVector(const Eigen::Vector3d & x)
1120  {
1121  Vector X(3);
1122  for (int i = 0; i < 3; i++)
1123  X(i) = x(i);
1124  return X;
1125  }
1126 
1130  static Eigen::Vector3d toVec3(const Real3dPoint & x)
1131  {
1132  return Eigen::Vector3d(x(0),x(1),x(2));
1133  }
1134 
1139  static Real3dVector toReal3dVector(const Eigen::Vector3d & x)
1140  {
1141  return { x(0), x(1), x(2)};
1142  }
1143 
1144 
1151  {
1152  Vector n(3);
1153  n(0) = 0.;
1154  n(1) = 0.;
1155  n(2) = 0.;
1156  /* for (auto f : mySurfaceMesh->incidentFaces(v))
1157  n += vectorArea(f);
1158  return n.normalized();
1159  */
1160  auto faces = mySurfaceMesh->incidentFaces(v);
1161  for (auto f : faces)
1162  n += vectorArea(f);
1163 
1164  if (fabs(n.norm() - 0.0) < 0.00001)
1165  {
1166  //On non-manifold edges touching the boundary, n may be null.
1167  trace.warning()<<"[PolygonalCalculus] Trying to compute the normal vector at a boundary vertex incident to pnon-manifold edge, we return a random vector."<<std::endl;
1168  n << Vector::Random(3);
1169  }
1170  n = n.normalized();
1171  return n;
1172  }
1173 
1176  Eigen::Vector3d n_v(const Vertex & v) const
1177  {
1178  return toVec3(myVertexNormalEmbedder(v));
1179  }
1180 
1181  //Covariant operators routines
1184  {
1185  auto nf = degree(f);
1186  DenseMatrix RU_fO = DenseMatrix::Zero(nf * 2,nf * 2);
1187  size_t cpt = 0;
1188  for (auto v : getSurfaceMeshPtr()->incidentVertices(f))
1189  {
1190  auto Rv = Rvf(v,f);
1191  RU_fO.block<2,2>(2 * cpt,2 * cpt) = Rv;
1192  ++cpt;
1193  }
1194  return RU_fO;
1195  }
1196 
1199  {
1200  size_t h = M.rows();
1201  size_t w = M.cols();
1202  DenseMatrix MK = DenseMatrix::Zero(h * 2,w * 2);
1203  for (size_t j = 0; j < h; j++)
1204  for (size_t i = 0; i < w; i++)
1205  {
1206  MK(2 * j, 2 * i) = M(j, i);
1207  MK(2 * j + 1, 2 * i + 1) = M(j, i);
1208  }
1209  return MK;
1210  }
1211 
1212 
1213 
1215 
1216  // ------------------------- Internals ------------------------------------
1217  //MARK: Internals
1218 private:
1219 
1222 
1224  std::function<Real3dPoint(Face, Vertex)> myEmbedder;
1225 
1228 
1230  std::vector<size_t> myFaceDegree;
1231 
1234  mutable std::array<std::unordered_map<Face,DenseMatrix>, 15> myGlobalCache;
1235 
1236 }; // end of class PolygonalCalculus
1237 
1244 template <typename TP, typename TV>
1245 std::ostream&
1246 operator<< ( std::ostream & out, const PolygonalCalculus<TP,TV> & object )
1247 {
1248  object.selfDisplay( out );
1249  return out;
1250 }
1251 
1252 } // namespace DGtal
DGtal::SurfaceMesh::Vertex
Index Vertex
Definition: SurfaceMesh.h:108
DGtal::PolygonalCalculus::Face
MySurfaceMesh::Face Face
Face type.
Definition: PolygonalCalculus.h:87
DGtal::PolygonalCalculus::getOperatorCacheMatrix
std::vector< DenseMatrix > getOperatorCacheMatrix(const std::function< DenseMatrix(Face)> &perFaceOperator) const
Definition: PolygonalCalculus.h:942
DGtal::PolygonalCalculus::DIVERGENCE_
@ DIVERGENCE_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::updateFaceDegree
void updateFaceDegree()
Update the face degree cache.
Definition: PolygonalCalculus.h:1072
DGtal::PolygonalCalculus::vectorArea
Vector vectorArea(const Face f) const
Definition: PolygonalCalculus.h:316
DGtal::PolygonalCalculus::A_
@ A_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::Form
Vector Form
Global 0-form, 1-form, 2-form are Vector.
Definition: PolygonalCalculus.h:98
DGtal::SurfaceMesh::position
RealPoint & position(Vertex v)
Definition: SurfaceMesh.h:637
DGtal::ConstAlias
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:186
DGtal::PolygonalCalculus::toExtrinsicVectors
std::vector< Vector > toExtrinsicVectors(const std::vector< Vector > &I) const
Definition: PolygonalCalculus.h:621
DGtal::SurfaceMesh::RealVector
TRealVector RealVector
Definition: SurfaceMesh.h:94
DGtal::PolygonalCalculus::faceNormalAsDGtalVector
Real3dVector faceNormalAsDGtalVector(const Face f) const
Definition: PolygonalCalculus.h:358
DGtal::PolygonalCalculus::X_
@ X_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::globalLaplaceBeltrami
SparseMatrix globalLaplaceBeltrami(const double lambda=1.0) const
Definition: PolygonalCalculus.h:798
DGtal::PolygonalCalculus::selfDisplay
void selfDisplay(std::ostream &out) const
Definition: PolygonalCalculus.h:1040
DGtal::PolygonalCalculus::SparseMatrix
LinAlg::SparseMatrix SparseMatrix
Type of sparse matrix.
Definition: PolygonalCalculus.h:102
DGtal::PolygonalCalculus::L_
@ L_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::PolygonalCalculus
PolygonalCalculus(const ConstAlias< MySurfaceMesh > surf, bool globalInternalCacheEnabled=false)
Definition: PolygonalCalculus.h:116
DGtal::PolygonalCalculus::transportAndFormatVectorField
DenseMatrix transportAndFormatVectorField(const Face f, const Vector &uf)
Definition: PolygonalCalculus.h:677
DGtal::PolygonalCalculus::centroidAsDGtalPoint
Real3dPoint centroidAsDGtalPoint(const Face f) const
Definition: PolygonalCalculus.h:436
DGtal::trace
Trace trace
Definition: Common.h:154
DGtal::PolygonalCalculus::setInCache
void setInCache(OPERATOR key, const Face f, const DenseMatrix &ope) const
Definition: PolygonalCalculus.h:1099
DGtal::PolygonalCalculus::PolygonalCalculus
PolygonalCalculus(const ConstAlias< MySurfaceMesh > surf, const std::function< Real3dPoint(Face, Vertex)> &pos_embedder, const std::function< Vector(Vertex)> &normal_embedder, bool globalInternalCacheEnabled=false)
Definition: PolygonalCalculus.h:166
DGtal::Dimension
DGtal::uint32_t Dimension
Definition: Common.h:137
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::PolygonalCalculus::divergence
DenseMatrix divergence(const Face f, const double lambda=1.0) const
Definition: PolygonalCalculus.h:505
DGtal::PolygonalCalculus::faceArea
double faceArea(const Face f) const
Definition: PolygonalCalculus.h:340
DGtal::PolygonalCalculus::dimension
static const Dimension dimension
Concept checking.
Definition: PolygonalCalculus.h:76
DGtal::PolygonalCalculus::OPERATOR
OPERATOR
Enum for operators in the internal cache strategy.
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::myGlobalCache
std::array< std::unordered_map< Face, DenseMatrix >, 15 > myGlobalCache
Definition: PolygonalCalculus.h:1234
DGtal::PolygonalCalculus::globalInverseLumpedMassMatrix
SparseMatrix globalInverseLumpedMassMatrix() const
Definition: PolygonalCalculus.h:845
DGtal::PolygonalCalculus::nbFaces
size_t nbFaces() const
Definition: PolygonalCalculus.h:1018
DGtal::PolygonalCalculus::getSurfaceMeshPtr
const MySurfaceMesh * getSurfaceMeshPtr() const
Definition: PolygonalCalculus.h:1031
DGtal::PolygonalCalculus::P_
@ P_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::setEmbedder
void setEmbedder(const std::function< Real3dPoint(Face, Vertex)> &externalFunctor)
Definition: PolygonalCalculus.h:222
DGtal::PolygonalCalculus
Implements differential operators on polygonal surfaces from .
Definition: PolygonalCalculus.h:70
DGtal::PolygonalCalculus::doubledGlobalLumpedMassMatrix
SparseMatrix doubledGlobalLumpedMassMatrix() const
Definition: PolygonalCalculus.h:903
DGtal::EigenLinearAlgebraBackend::DenseVector
Eigen::VectorXd DenseVector
Definition: EigenSupport.h:98
DGtal::PolygonalCalculus::SHARP_
@ SHARP_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::covariantProjection
DenseMatrix covariantProjection(const Face f, const Vector &uf)
Definition: PolygonalCalculus.h:715
DGtal::EigenLinearAlgebraBackend::DenseMatrix
Eigen::MatrixXd DenseMatrix
Definition: EigenSupport.h:99
DGtal::PolygonalCalculus::toVector
static Vector toVector(const Eigen::Vector3d &x)
toVector convert Real3dPoint to Eigen::VectorXd
Definition: PolygonalCalculus.h:1119
DGtal::PolygonalCalculus::coGradient
DenseMatrix coGradient(const Face f) const
Definition: PolygonalCalculus.h:367
DGtal::EigenLinearAlgebraBackend::Triplet
Eigen::Triplet< double, SparseMatrix::StorageIndex > Triplet
Definition: EigenSupport.h:103
DGtal::PolygonalCalculus::globalConnectionLaplace
SparseMatrix globalConnectionLaplace(const double lambda=1.0) const
Definition: PolygonalCalculus.h:871
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
DGtal::operator<<
std::ostream & operator<<(std::ostream &out, const ATu0v1< TKSpace, TLinearAlgebra > &object)
DGtal::PolygonalCalculus::DenseMatrix
LinAlg::DenseMatrix DenseMatrix
Type of dense matrix.
Definition: PolygonalCalculus.h:100
DGtal::SurfaceMesh::incidentFaces
const Faces & incidentFaces(Vertex v) const
Definition: SurfaceMesh.h:313
DGtal::PolygonalCalculus::Triplet
LinAlg::Triplet Triplet
Type of sparse matrix triplet.
Definition: PolygonalCalculus.h:104
DGtal::PolygonalCalculus::E
DenseMatrix E(const Face f) const
Definition: PolygonalCalculus.h:281
DGtal::PolygonalCalculus::P
DenseMatrix P(const Face f) const
Definition: PolygonalCalculus.h:461
DGtal::PolygonalCalculus::D
DenseMatrix D(const Face f) const
Definition: PolygonalCalculus.h:261
DGtal::PolygonalCalculus::M_
@ M_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::checkCache
bool checkCache(OPERATOR key, const Face f) const
Definition: PolygonalCalculus.h:1087
DGtal::PolygonalCalculus::myVertexNormalEmbedder
std::function< Real3dVector(Vertex)> myVertexNormalEmbedder
Embedding function (vertex)->R^3 for the vertex normal.
Definition: PolygonalCalculus.h:1227
DGtal::PolygonalCalculus::D_
@ D_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::B_
@ B_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::toReal3dVector
static Real3dVector toReal3dVector(const Eigen::Vector3d &x)
toReal3dVector converts Eigen::Vector3d to Real3dVector.
Definition: PolygonalCalculus.h:1139
DGtal::PolygonalCalculus::sharp
DenseMatrix sharp(const Face f) const
Definition: PolygonalCalculus.h:445
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::PolygonalCalculus::MySurfaceMesh
SurfaceMesh< TRealPoint, TRealVector > MySurfaceMesh
Type of SurfaceMesh.
Definition: PolygonalCalculus.h:83
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::PolygonalCalculus::FLAT_
@ FLAT_
Definition: PolygonalCalculus.h:1069
DGtal::EigenLinearAlgebraBackend
Aim: Provide linear algebra backend using Eigen dense and sparse matrix as well as dense vector....
Definition: EigenSupport.h:96
DGtal::PolygonalCalculus::LinAlg
EigenLinearAlgebraBackend LinAlg
Linear Algebra Backend from Eigen.
Definition: PolygonalCalculus.h:94
DGtal::PolygonalCalculus::computeVertexNormal
Vector computeVertexNormal(const Vertex &v) const
Definition: PolygonalCalculus.h:1150
DGtal::PolygonalCalculus::PolygonalCalculus
PolygonalCalculus(const ConstAlias< MySurfaceMesh > surf, const std::function< Real3dVector(Vertex)> &embedder, bool globalInternalCacheEnabled=false)
Definition: PolygonalCalculus.h:146
DGtal::PolygonalCalculus::init
void init()
Definition: PolygonalCalculus.h:998
DGtal::PolygonalCalculus::Tf
DenseMatrix Tf(const Face &f) const
Definition: PolygonalCalculus.h:587
DGtal::PolygonalCalculus::isValid
bool isValid() const
Definition: PolygonalCalculus.h:1054
DGtal::PolygonalCalculus::Solver
LinAlg::SolverSimplicialLDLT Solver
Type of a sparse matrix solver.
Definition: PolygonalCalculus.h:107
DGtal::PolygonalCalculus::CURL_
@ CURL_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::project
static Vector project(const Vector &u, const Vector &n)
Definition: PolygonalCalculus.h:1110
DGtal::PolygonalCalculus::myEmbedder
std::function< Real3dPoint(Face, Vertex)> myEmbedder
Embedding function (face,vertex)->R^3 for the vertex position wrt. the face.
Definition: PolygonalCalculus.h:1224
DGtal::PolygonalCalculus::curl
DenseMatrix curl(const Face f) const
Definition: PolygonalCalculus.h:519
DGtal::PolygonalCalculus::toExtrinsicVector
Vector toExtrinsicVector(const Vertex v, const Vector &I) const
toExtrinsicVector
Definition: PolygonalCalculus.h:611
DGtal::PolygonalCalculus::disableInternalGlobalCache
void disableInternalGlobalCache()
Definition: PolygonalCalculus.h:983
DGtal::PolygonalCalculus::Real3dPoint
MySurfaceMesh::RealPoint Real3dPoint
Position type.
Definition: PolygonalCalculus.h:89
DGtal::PolygonalCalculus::BOOST_STATIC_ASSERT
BOOST_STATIC_ASSERT((dimension==3))
DGtal::PolygonalCalculus::blockConnection
DenseMatrix blockConnection(const Face &f) const
Definition: PolygonalCalculus.h:1183
DGtal::PolygonalCalculus::Real3dVector
MySurfaceMesh::RealVector Real3dVector
Real vector type.
Definition: PolygonalCalculus.h:91
DGtal::SurfaceMesh::Face
Index Face
Definition: SurfaceMesh.h:106
DGtal::PolygonalCalculus::mySurfaceMesh
const MySurfaceMesh * mySurfaceMesh
Underlying SurfaceMesh.
Definition: PolygonalCalculus.h:1221
DGtal::PolygonalCalculus::nbVertices
size_t nbVertices() const
Definition: PolygonalCalculus.h:1012
DGtal::EigenLinearAlgebraBackend::SolverSimplicialLDLT
Eigen::SimplicialLDLT< SparseMatrix > SolverSimplicialLDLT
Definition: EigenSupport.h:107
DGtal::PolygonalCalculus::M
DenseMatrix M(const Face f, const double lambda=1.0) const
Definition: PolygonalCalculus.h:477
DGtal::PolygonalCalculus::X
DenseMatrix X(const Face f) const
Definition: PolygonalCalculus.h:236
DGtal::PolygonalCalculus::toVec3
static Eigen::Vector3d toVec3(const Real3dPoint &x)
toVec3 convert Real3dPoint to Eigen::Vector3d
Definition: PolygonalCalculus.h:1130
DGtal::PolygonalCalculus::Rvf
DenseMatrix Rvf(const Vertex &v, const Face &f) const
Definition: PolygonalCalculus.h:649
DGtal::PolygonalCalculus::Vertex
MySurfaceMesh::Vertex Vertex
Vertex type.
Definition: PolygonalCalculus.h:85
DGtal::PolygonalCalculus::Qvf
DenseMatrix Qvf(const Vertex &v, const Face &f) const
Definition: PolygonalCalculus.h:631
DGtal::PolygonalCalculus::PolygonalCalculus
PolygonalCalculus()=delete
DGtal::PolygonalCalculus::B
DenseMatrix B(const Face f) const
Definition: PolygonalCalculus.h:417
DGtal::PolygonalCalculus::n_v
Eigen::Vector3d n_v(const Vertex &v) const
Definition: PolygonalCalculus.h:1176
DGtal::SurfaceMesh::neighborVertices
const Vertices & neighborVertices(Vertex v) const
Definition: SurfaceMesh.h:323
DGtal::PolygonalCalculus::Vector
LinAlg::DenseVector Vector
Type of Vector.
Definition: PolygonalCalculus.h:96
DGtal::PolygonalCalculus::identity0
SparseMatrix identity0() const
Definition: PolygonalCalculus.h:776
DGtal::PolygonalCalculus::myFaceDegree
std::vector< size_t > myFaceDegree
Cache containing the face degree.
Definition: PolygonalCalculus.h:1230
DGtal::SurfaceMesh::RealPoint
TRealPoint RealPoint
Definition: SurfaceMesh.h:93
DGtal::PolygonalCalculus::myGlobalCacheEnabled
bool myGlobalCacheEnabled
Global cache.
Definition: PolygonalCalculus.h:1233
DGtal::PolygonalCalculus::COGRAD_
@ COGRAD_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::faceNormal
Vector faceNormal(const Face f) const
Definition: PolygonalCalculus.h:348
DGtal::PolygonalCalculus::GRAD_
@ GRAD_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::Self
PolygonalCalculus< TRealPoint, TRealVector > Self
Self type.
Definition: PolygonalCalculus.h:80
DGtal::EigenLinearAlgebraBackend::SparseMatrix
Eigen::SparseMatrix< DenseVector::Scalar, Eigen::ColMajor, DenseVector::Index > SparseMatrix
Definition: EigenSupport.h:102
DGtal::PolygonalCalculus::flat
DenseMatrix flat(const Face f) const
Definition: PolygonalCalculus.h:404
DGtal::PolygonalCalculus::~PolygonalCalculus
~PolygonalCalculus()=default
DGtal::PolygonalCalculus::operator=
PolygonalCalculus & operator=(const PolygonalCalculus &other)=delete
DGtal::Trace::warning
std::ostream & warning()
DGtal::PolygonalCalculus::CON_L_
@ CON_L_
Definition: PolygonalCalculus.h:1069
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::PolygonalCalculus::kroneckerWithI2
DenseMatrix kroneckerWithI2(const DenseMatrix &M) const
Definition: PolygonalCalculus.h:1198
DGtal::SurfaceMesh::incidentVertices
const Vertices & incidentVertices(Face f) const
Definition: SurfaceMesh.h:307
DGtal::PolygonalCalculus::PolygonalCalculus
PolygonalCalculus(const ConstAlias< MySurfaceMesh > surf, const std::function< Real3dPoint(Face, Vertex)> &embedder, bool globalInternalCacheEnabled=false)
Definition: PolygonalCalculus.h:131
DGtal::PolygonalCalculus::bracket
DenseMatrix bracket(const Vector &n) const
Definition: PolygonalCalculus.h:378
DGtal::PolygonalCalculus::covariantGradient_f
DenseMatrix covariantGradient_f(const Face &f) const
Definition: PolygonalCalculus.h:725
DGtal::PolygonalCalculus::enableInternalGlobalCache
void enableInternalGlobalCache()
Definition: PolygonalCalculus.h:976
DGtal::PolygonalCalculus::shapeOperator
DenseMatrix shapeOperator(const Face f) const
Definition: PolygonalCalculus.h:656
DGtal::SurfaceMesh::nbFaces
Size nbFaces() const
Definition: SurfaceMesh.h:288
DGtal::PolygonalCalculus::E_
@ E_
Definition: PolygonalCalculus.h:1069
DGtal::PolygonalCalculus::covariantProjection_f
DenseMatrix covariantProjection_f(const Face &f) const
Definition: PolygonalCalculus.h:735
DGtal::PolygonalCalculus::Tv
DenseMatrix Tv(const Vertex &v) const
Definition: PolygonalCalculus.h:567
DGtal::PolygonalCalculus::A
DenseMatrix A(const Face f) const
Definition: PolygonalCalculus.h:295
DGtal::PolygonalCalculus::degree
size_t degree(const Face f) const
Definition: PolygonalCalculus.h:1025