$darkmode
VCG Library
curvature.h
1 /****************************************************************************
2 * VCGLib o o *
3 * Visual and Computer Graphics Library o o *
4 * _ O _ *
5 * Copyright(C) 2004-2016 \/)\/ *
6 * Visual Computing Lab /\/| *
7 * ISTI - Italian National Research Council | *
8 * \ *
9 * All rights reserved. *
10 * *
11 * This program is free software; you can redistribute it and/or modify *
12 * it under the terms of the GNU General Public License as published by *
13 * the Free Software Foundation; either version 2 of the License, or *
14 * (at your option) any later version. *
15 * *
16 * This program is distributed in the hope that it will be useful, *
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
19 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
20 * for more details. *
21 * *
22 ****************************************************************************/
23 
24 #ifndef VCGLIB_UPDATE_CURVATURE_
25 #define VCGLIB_UPDATE_CURVATURE_
26 
27 #include <vcg/space/index/grid_static_ptr.h>
28 #include <vcg/simplex/face/jumping_pos.h>
29 #include <vcg/complex/algorithms/update/normal.h>
30 #include <vcg/complex/algorithms/point_sampling.h>
31 #include <vcg/complex/algorithms/intersection.h>
32 #include <vcg/complex/algorithms/inertia.h>
33 #include <Eigen/Core>
34 
35 namespace vcg {
36 namespace tri {
37 
39 
41 
43 
47 template <class MeshType>
48 class UpdateCurvature
49 {
50 
51 public:
52  typedef typename MeshType::FaceType FaceType;
53  typedef typename MeshType::FacePointer FacePointer;
54  typedef typename MeshType::FaceIterator FaceIterator;
55  typedef typename MeshType::VertexIterator VertexIterator;
56  typedef typename MeshType::VertContainer VertContainer;
57  typedef typename MeshType::VertexType VertexType;
58  typedef typename MeshType::VertexPointer VertexPointer;
59  typedef vcg::face::VFIterator<FaceType> VFIteratorType;
60  typedef typename MeshType::CoordType CoordType;
61  typedef typename CoordType::ScalarType ScalarType;
62  typedef typename MeshType::VertexType::CurScalarType CurScalarType;
63  typedef typename MeshType::VertexType::CurVecType CurVecType;
64 
65 
66 private:
67  // aux data struct used by PrincipalDirections
68  struct AdjVertex {
69  VertexType * vert;
70  float doubleArea;
71  bool isBorder;
72  };
73 
74 
75 public:
77 
78 /*
79  Compute principal direction and magniuto of curvature as describe in the paper:
80  @InProceedings{bb33922,
81  author = "G. Taubin",
82  title = "Estimating the Tensor of Curvature of a Surface from a
83  Polyhedral Approximation",
84  booktitle = "International Conference on Computer Vision",
85  year = "1995",
86  pages = "902--907",
87  URL = "http://dx.doi.org/10.1109/ICCV.1995.466840",
88  bibsource = "http://www.visionbib.com/bibliography/describe440.html#TT32253",
89  */
90  static void PrincipalDirections(MeshType &m)
91  {
92  tri::RequireVFAdjacency(m);
93  vcg::tri::UpdateNormal<MeshType>::PerVertexAngleWeighted(m);
94  vcg::tri::UpdateNormal<MeshType>::NormalizePerVertex(m);
95 
96  for (VertexIterator vi =m.vert.begin(); vi !=m.vert.end(); ++vi) {
97  if ( ! (*vi).IsD() && (*vi).VFp() != NULL) {
98 
99  VertexType * central_vertex = &(*vi);
100 
101  std::vector<float> weights;
102  std::vector<AdjVertex> vertices;
103 
104  vcg::face::JumpingPos<FaceType> pos((*vi).VFp(), central_vertex);
105 
106  // firstV is the first vertex of the 1ring neighboorhood
107  VertexType* firstV = pos.VFlip();
108  VertexType* tempV;
109  float totalDoubleAreaSize = 0.0f;
110 
111  // compute the area of each triangle around the central vertex as well as their total area
112  do
113  {
114  // this bring the pos to the next triangle counterclock-wise
115  pos.FlipF();
116  pos.FlipE();
117 
118  // tempV takes the next vertex in the 1ring neighborhood
119  tempV = pos.VFlip();
120  assert(tempV!=central_vertex);
121  AdjVertex v;
122 
123  v.isBorder = pos.IsBorder();
124  v.vert = tempV;
125  v.doubleArea = vcg::DoubleArea(*pos.F());
126  totalDoubleAreaSize += v.doubleArea;
127 
128  vertices.push_back(v);
129  }
130  while(tempV != firstV);
131 
132  // compute the weights for the formula computing matrix M
133  for (size_t i = 0; i < vertices.size(); ++i) {
134  if (vertices[i].isBorder) {
135  weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize);
136  } else {
137  weights.push_back(0.5f * (vertices[i].doubleArea + vertices[(i-1)%vertices.size()].doubleArea) / totalDoubleAreaSize);
138  }
139  assert(weights.back() < 1.0f);
140  }
141 
142  // compute I-NN^t to be used for computing the T_i's
143  Matrix33<ScalarType> Tp;
144  for (int i = 0; i < 3; ++i)
145  Tp[i][i] = 1.0f - powf(central_vertex->cN()[i],2);
146  Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->cN()[1]);
147  Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->cN()[1] * central_vertex->cN()[2]);
148  Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->cN()[0] * central_vertex->cN()[2]);
149 
150  // for all neighbors vi compute the directional curvatures k_i and the T_i
151  // compute M by summing all w_i k_i T_i T_i^t
152  Matrix33<ScalarType> tempMatrix;
153  Matrix33<ScalarType> M;
154  M.SetZero();
155  for (size_t i = 0; i < vertices.size(); ++i) {
156  CoordType edge = (central_vertex->cP() - vertices[i].vert->cP());
157  float curvature = (2.0f * (central_vertex->cN().dot(edge)) ) / edge.SquaredNorm();
158  CoordType T = (Tp*edge).normalized();
159  tempMatrix.ExternalProduct(T,T);
160  M += tempMatrix * weights[i] * curvature ;
161  }
162 
163  // compute vector W for the Householder matrix
164  CoordType W;
165  CoordType e1(1.0f,0.0f,0.0f);
166  if ((e1 - central_vertex->cN()).SquaredNorm() > (e1 + central_vertex->cN()).SquaredNorm())
167  W = e1 - central_vertex->cN();
168  else
169  W = e1 + central_vertex->cN();
170  W.Normalize();
171 
172  // compute the Householder matrix I - 2WW^t
173  Matrix33<ScalarType> Q;
174  Q.SetIdentity();
175  tempMatrix.ExternalProduct(W,W);
176  Q -= tempMatrix * 2.0f;
177 
178  // compute matrix Q^t M Q
179  Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
180 
181 // CoordType bl = Q.GetColumn(0);
182  CoordType T1 = Q.GetColumn(1);
183  CoordType T2 = Q.GetColumn(2);
184 
185  // find sin and cos for the Givens rotation
186  float s,c;
187  // Gabriel Taubin hint and Valentino Fiorin impementation
188  float alpha = QtMQ[1][1]-QtMQ[2][2];
189  float beta = QtMQ[2][1];
190 
191  float h[2];
192  float delta = sqrtf(4.0f*powf(alpha, 2) +16.0f*powf(beta, 2));
193  h[0] = (2.0f*alpha + delta) / (2.0f*beta);
194  h[1] = (2.0f*alpha - delta) / (2.0f*beta);
195 
196  float t[2];
197  float best_c, best_s;
198  float min_error = std::numeric_limits<ScalarType>::infinity();
199  for (int i=0; i<2; i++)
200  {
201  delta = sqrtf(powf(h[i], 2) + 4.0f);
202  t[0] = (h[i]+delta) / 2.0f;
203  t[1] = (h[i]-delta) / 2.0f;
204 
205  for (int j=0; j<2; j++)
206  {
207  float squared_t = powf(t[j], 2);
208  float denominator = 1.0f + squared_t;
209  s = (2.0f*t[j]) / denominator;
210  c = (1-squared_t) / denominator;
211 
212  float approximation = c*s*alpha + (powf(c, 2) - powf(s, 2))*beta;
213  float angle_similarity = fabs(acosf(c)/asinf(s));
214  float error = fabs(1.0f-angle_similarity)+fabs(approximation);
215  if (error<min_error)
216  {
217  min_error = error;
218  best_c = c;
219  best_s = s;
220  }
221  }
222  }
223  c = best_c;
224  s = best_s;
225 
226  Eigen::Matrix2f minor2x2;
227  Eigen::Matrix2f S;
228 
229 
230  // diagonalize M
231  minor2x2(0,0) = QtMQ[1][1];
232  minor2x2(0,1) = QtMQ[1][2];
233  minor2x2(1,0) = QtMQ[2][1];
234  minor2x2(1,1) = QtMQ[2][2];
235 
236  S(0,0) = S(1,1) = c;
237  S(0,1) = s;
238  S(1,0) = -1.0f * s;
239 
240  Eigen::Matrix2f StMS = S.transpose() * minor2x2 * S;
241 
242  // compute curvatures and curvature directions
243  float Principal_Curvature1 = (3.0f * StMS(0,0)) - StMS(1,1);
244  float Principal_Curvature2 = (3.0f * StMS(1,1)) - StMS(0,0);
245 
246  CoordType Principal_Direction1 = T1 * c - T2 * s;
247  CoordType Principal_Direction2 = T1 * s + T2 * c;
248 
249  (*vi).PD1().Import(Principal_Direction1);
250  (*vi).PD2().Import(Principal_Direction2);
251  (*vi).K1() = Principal_Curvature1;
252  (*vi).K2() = Principal_Curvature2;
253  }
254  }
255  }
256 
257 
258 
259 
260  class AreaData
261  {
262  public:
263  float A;
264  };
265 
273  typedef vcg::GridStaticPtr <FaceType, ScalarType > MeshGridType;
274  typedef vcg::GridStaticPtr <VertexType, ScalarType > PointsGridType;
275 
276  static void PrincipalDirectionsPCA(MeshType &m, ScalarType r, bool pointVSfaceInt = true,vcg::CallBackPos * cb = NULL)
277  {
278  std::vector<VertexType*> closests;
279  std::vector<ScalarType> distances;
280  std::vector<CoordType> points;
281  VertexIterator vi;
282  ScalarType area;
283  MeshType tmpM;
284  typename std::vector<CoordType>::iterator ii;
285  vcg::tri::TrivialSampler<MeshType> vs;
286  tri::UpdateNormal<MeshType>::PerVertexAngleWeighted(m);
287  tri::UpdateNormal<MeshType>::NormalizePerVertex(m);
288 
289  MeshGridType mGrid;
290  PointsGridType pGrid;
291 
292  // Fill the grid used
293  if(pointVSfaceInt)
294  {
295  area = Stat<MeshType>::ComputeMeshArea(m);
296  vcg::tri::SurfaceSampling<MeshType,vcg::tri::TrivialSampler<MeshType> >::Montecarlo(m,vs,1000 * area / (2*M_PI*r*r ));
297  vi = vcg::tri::Allocator<MeshType>::AddVertices(tmpM,m.vert.size());
298  for(size_t y = 0; y < m.vert.size(); ++y,++vi) (*vi).P() = m.vert[y].P();
299  pGrid.Set(tmpM.vert.begin(),tmpM.vert.end());
300  }
301  else
302  {
303  mGrid.Set(m.face.begin(),m.face.end());
304  }
305 
306  int jj = 0;
307  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
308  {
309  vcg::Matrix33<ScalarType> A, eigenvectors;
310  vcg::Point3<ScalarType> bp, eigenvalues;
311 // int nrot;
312 
313  // sample the neighborhood
314  if(pointVSfaceInt)
315  {
316  vcg::tri::GetInSphereVertex<
317  MeshType,
318  PointsGridType,std::vector<VertexType*>,
319  std::vector<ScalarType>,
320  std::vector<CoordType> >(tmpM,pGrid, (*vi).cP(),r ,closests,distances,points);
321 
322  A.Covariance(points,bp);
323  A*=area*area/1000;
324  }
325  else{
326  IntersectionBallMesh<MeshType,ScalarType>( m ,vcg::Sphere3<ScalarType>((*vi).cP(),r),tmpM );
327  vcg::Point3<ScalarType> _bary;
328  vcg::tri::Inertia<MeshType>::Covariance(tmpM,_bary,A);
329  }
330 
331 // Eigen::Matrix3f AA; AA=A;
332 // Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eig(AA);
333 // Eigen::Vector3f c_val = eig.eigenvalues();
334 // Eigen::Matrix3f c_vec = eig.eigenvectors();
335 
336 // Jacobi(A, eigenvalues , eigenvectors, nrot);
337 
338  Eigen::Matrix3d AA;
339  A.ToEigenMatrix(AA);
340  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(AA);
341  Eigen::Vector3d c_val = eig.eigenvalues();
342  Eigen::Matrix3d c_vec = eig.eigenvectors(); // eigenvector are stored as columns.
343  eigenvectors.FromEigenMatrix(c_vec);
344  eigenvalues.FromEigenVector(c_val);
345 // EV.transposeInPlace();
346 // ev.FromEigenVector(c_val);
347 
348  // get the estimate of curvatures from eigenvalues and eigenvectors
349  // find the 2 most tangent eigenvectors (by finding the one closest to the normal)
350  int best = 0; ScalarType bestv = fabs( (*vi).cN().dot(eigenvectors.GetColumn(0).normalized()) );
351  for(int i = 1 ; i < 3; ++i){
352  ScalarType prod = fabs((*vi).cN().dot(eigenvectors.GetColumn(i).normalized()));
353  if( prod > bestv){bestv = prod; best = i;}
354  }
355 
356  (*vi).PD1().Import(eigenvectors.GetColumn( (best+1)%3).normalized());
357  (*vi).PD2().Import(eigenvectors.GetColumn( (best+2)%3).normalized());
358 
359  // project them to the plane identified by the normal
360  vcg::Matrix33<CurScalarType> rot;
361  CurVecType NN = CurVecType::Construct((*vi).N());
362  CurScalarType angle;
363  angle = acos((*vi).PD1().dot(NN));
364  rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD1()^NN);
365  (*vi).PD1() = rot*(*vi).PD1();
366  angle = acos((*vi).PD2().dot(NN));
367  rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD2()^NN);
368  (*vi).PD2() = rot*(*vi).PD2();
369 
370 
371  // copmutes the curvature values
372  const ScalarType r5 = r*r*r*r*r;
373  const ScalarType r6 = r*r5;
374  (*vi).K1() = (2.0/5.0) * (4.0*M_PI*r5 + 15*eigenvalues[(best+2)%3]-45.0*eigenvalues[(best+1)%3])/(M_PI*r6);
375  (*vi).K2() = (2.0/5.0) * (4.0*M_PI*r5 + 15*eigenvalues[(best+1)%3]-45.0*eigenvalues[(best+2)%3])/(M_PI*r6);
376  if((*vi).K1() < (*vi).K2()) { std::swap((*vi).K1(),(*vi).K2());
377  std::swap((*vi).PD1(),(*vi).PD2());
378  if (cb)
379  {
380  (*cb)(int(100.0f * (float)jj / (float)m.vn),"Vertices Analysis");
381  ++jj;
382  } }
383  }
384 
385 
386  }
387 
389 
400 static void MeanAndGaussian(MeshType & m)
401 {
402  tri::RequireFFAdjacency(m);
403 
404  float area0, area1, area2, angle0, angle1, angle2;
405  FaceIterator fi;
406  VertexIterator vi;
407  typename MeshType::CoordType e01v ,e12v ,e20v;
408 
409  SimpleTempData<VertContainer, AreaData> TDAreaPtr(m.vert);
410  SimpleTempData<VertContainer, typename MeshType::CoordType> TDContr(m.vert);
411 
412  vcg::tri::UpdateNormal<MeshType>::PerVertexNormalized(m);
413  auto KH = vcg::tri::Allocator<MeshType>:: template GetPerVertexAttribute<ScalarType> (m, std::string("KH"));
414  auto KG = vcg::tri::Allocator<MeshType>:: template GetPerVertexAttribute<ScalarType> (m, std::string("KG"));
415 
416  //Compute AreaMix in H (vale anche per K)
417  for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD())
418  {
419  (TDAreaPtr)[*vi].A = 0.0;
420  (TDContr)[*vi] =typename MeshType::CoordType(0.0,0.0,0.0);
421  KH[*vi] = 0.0;
422  KG[*vi] = (ScalarType)(2.0 * M_PI);
423  }
424 
425  for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD())
426  {
427  // angles
428  angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
429  angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
430  angle2 = M_PI-(angle0+angle1);
431 
432  if((angle0 < M_PI/2) && (angle1 < M_PI/2) && (angle2 < M_PI/2)) // triangolo non ottuso
433  {
434  float e01 = SquaredDistance( (*fi).V(1)->cP() , (*fi).V(0)->cP() );
435  float e12 = SquaredDistance( (*fi).V(2)->cP() , (*fi).V(1)->cP() );
436  float e20 = SquaredDistance( (*fi).V(0)->cP() , (*fi).V(2)->cP() );
437 
438  area0 = ( e20*(1.0/tan(angle1)) + e01*(1.0/tan(angle2)) ) / 8.0;
439  area1 = ( e01*(1.0/tan(angle2)) + e12*(1.0/tan(angle0)) ) / 8.0;
440  area2 = ( e12*(1.0/tan(angle0)) + e20*(1.0/tan(angle1)) ) / 8.0;
441 
442  (TDAreaPtr)[(*fi).V(0)].A += area0;
443  (TDAreaPtr)[(*fi).V(1)].A += area1;
444  (TDAreaPtr)[(*fi).V(2)].A += area2;
445 
446  }
447  else // obtuse
448  {
449  if(angle0 >= M_PI/2)
450  {
451  (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
452  (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
453  (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
454  }
455  else if(angle1 >= M_PI/2)
456  {
457  (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
458  (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
459  (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
460  }
461  else
462  {
463  (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
464  (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
465  (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
466  }
467  }
468  }
469 
470 
471  for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD() )
472  {
473  angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
474  angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
475  angle2 = M_PI-(angle0+angle1);
476 
477  // Skip degenerate triangles.
478  if(angle0==0 || angle1==0 || angle2==0) continue;
479 
480  e01v = ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ;
481  e12v = ( (*fi).V(2)->cP() - (*fi).V(1)->cP() ) ;
482  e20v = ( (*fi).V(0)->cP() - (*fi).V(2)->cP() ) ;
483 
484  TDContr[(*fi).V(0)] += ( e20v * (1.0/tan(angle1)) - e01v * (1.0/tan(angle2)) ) / 4.0;
485  TDContr[(*fi).V(1)] += ( e01v * (1.0/tan(angle2)) - e12v * (1.0/tan(angle0)) ) / 4.0;
486  TDContr[(*fi).V(2)] += ( e12v * (1.0/tan(angle0)) - e20v * (1.0/tan(angle1)) ) / 4.0;
487 
488  KG[(*fi).V(0)] -= angle0;
489  KG[(*fi).V(1)] -= angle1;
490  KG[(*fi).V(2)] -= angle2;
491 
492 
493  for(int i=0;i<3;i++)
494  {
495  if(vcg::face::IsBorder((*fi), i))
496  {
497  CoordType e1,e2;
498  vcg::face::Pos<FaceType> hp(&*fi, i, (*fi).V(i));
499  vcg::face::Pos<FaceType> hp1=hp;
500 
501  hp1.FlipV();
502  e1=hp1.v->cP() - hp.v->cP();
503  hp1.FlipV();
504  hp1.NextB();
505  e2=hp1.v->cP() - hp.v->cP();
506  KG[(*fi).V(i)] -= math::Abs(Angle(e1,e2));
507  }
508  }
509  }
510 
511  for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD() /*&& !(*vi).IsB()*/)
512  {
513  if((TDAreaPtr)[*vi].A<=std::numeric_limits<ScalarType>::epsilon())
514  {
515  KH[(*vi)] = 0;
516  KG[(*vi)] = 0;
517  }
518  else
519  {
520  KH[(*vi)] = (((TDContr)[*vi].dot((*vi).cN())>0)?1.0:-1.0)*((TDContr)[*vi] / (TDAreaPtr) [*vi].A).Norm();
521  KG[(*vi)] /= (TDAreaPtr)[*vi].A;
522  }
523  }
524 }
525 
526 
528 
536 static void PerVertexAbsoluteMeanAndGaussian(MeshType & m)
537 {
538  tri::RequireVFAdjacency(m);
539  tri::RequireCompactness(m);
540  const bool areaNormalize = true;
541  const bool barycentricArea=false;
542  auto KH = vcg::tri::Allocator<MeshType>:: template GetPerVertexAttribute<ScalarType> (m, std::string("KH"));
543  auto KG = vcg::tri::Allocator<MeshType>:: template GetPerVertexAttribute<ScalarType> (m, std::string("KG"));
544  int faceCnt=0;
545  for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
546  {
547  VertexPointer v=&*vi;
548  VFIteratorType vfi(v);
549  ScalarType A = 0;
550 
551  KH[v] = 0;
552  ScalarType AngleDefect = (ScalarType)(2.0 * M_PI);;
553 
554  while (!vfi.End()) {
555  faceCnt++;
556  FacePointer f = vfi.F();
557  CoordType nf = TriangleNormal(*f);
558  int i = vfi.I();
559  VertexPointer v0 = f->V0(i), v1 = f->V1(i), v2 = f->V2(i);
560  assert (v==v0);
561 
562  ScalarType ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P() ));
563  ScalarType ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P() ));
564  ScalarType ang2 = M_PI - ang0 - ang1;
565 
566  ScalarType s01 = SquaredDistance(v1->P(), v0->P());
567  ScalarType s02 = SquaredDistance(v2->P(), v0->P());
568 
569  // voronoi cell of current vertex
570  if(barycentricArea)
571  A+=vcg::DoubleArea(*f)/6.0;
572  else
573  {
574  if (ang0 >= M_PI/2)
575  A += (0.5f * DoubleArea(*f) - (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 );
576  else if (ang1 >= M_PI/2)
577  A += (s01 * tan(ang0)) / 8.0;
578  else if (ang2 >= M_PI/2)
579  A += (s02 * tan(ang0)) / 8.0;
580  else // non obctuse triangle
581  A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0;
582  }
583  // gaussian curvature update
584  AngleDefect -= ang0;
585 
586  // mean curvature update
587  // Note that the standard abs mean curvature approximation would require
588  // to sum all the edges*diehedralAngle. Here with just VF adjacency
589  // we make a rough approximation that 1/2 of the edge len plus something
590  // that is half of the diedral angle
591  ang1 = math::Abs(Angle(nf, v1->N()+v0->N()));
592  ang2 = math::Abs(Angle(nf, v2->N()+v0->N()));
593  KH[v] += math::Sqrt(s01)*ang1 + math::Sqrt(s02)*ang2 ;
594  ++vfi;
595  }
596 
597  KH[v] /= 4.0;
598 
599  if(areaNormalize) {
600  if(A <= std::numeric_limits<float>::epsilon()) {
601  KH[v] = 0;
602  KG[v] = 0;
603  }
604  else {
605  KH[v] /= A;
606  KG[v] = AngleDefect / A;
607  }
608  }
609  }
610 }
611 
612 
613 
614 /*
615  Compute principal curvature directions and value with normal cycle:
616  @inproceedings{CohMor03,
617  author = {Cohen-Steiner, David and Morvan, Jean-Marie },
618  booktitle = {SCG '03: Proceedings of the nineteenth annual symposium on Computational geometry},
619  title - {Restricted delaunay triangulations and normal cycle}
620  year = {2003}
621 }
622  */
623 
624  static void PrincipalDirectionsNormalCycle(MeshType & m){
625  tri::RequireVFAdjacency(m);
626  tri::RequireFFAdjacency(m);
627 
628  typename MeshType::VertexIterator vi;
629 
630  for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
631  if(!((*vi).IsD())){
632  vcg::Matrix33<ScalarType> m33;m33.SetZero();
633  face::JumpingPos<typename MeshType::FaceType> p((*vi).VFp(),&(*vi));
634  p.FlipE();
635  typename MeshType::VertexType * firstv = p.VFlip();
636  assert(p.F()->V(p.VInd())==&(*vi));
637 
638  do{
639  if( p.F() != p.FFlip()){
640  Point3<ScalarType> normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P();
641  ScalarType edge_length = normalized_edge.Norm();
642  normalized_edge/=edge_length;
643  Point3<ScalarType> n1 = NormalizedTriangleNormal(*(p.F()));
644  Point3<ScalarType> n2 = NormalizedTriangleNormal(*(p.FFlip()));
645  ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge);
646  n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0));
647  ScalarType beta = math::Asin(n1n2);
648  m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0];
649  m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0];
650  m33[1][1] += beta*edge_length*normalized_edge[1]*normalized_edge[1];
651  m33[0][2] += beta*edge_length*normalized_edge[2]*normalized_edge[0];
652  m33[1][2] += beta*edge_length*normalized_edge[2]*normalized_edge[1];
653  m33[2][2] += beta*edge_length*normalized_edge[2]*normalized_edge[2];
654  }
655  p.NextFE();
656  }while(firstv != p.VFlip());
657 
658  if(m33.Determinant()==0.0){ // degenerate case
659  (*vi).K1() = (*vi).K2() = 0.0; continue;}
660 
661  m33[1][0] = m33[0][1];
662  m33[2][0] = m33[0][2];
663  m33[2][1] = m33[1][2];
664 
665  Eigen::Matrix3d it;
666  m33.ToEigenMatrix(it);
667  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(it);
668  Eigen::Vector3d c_val = eig.eigenvalues();
669  Eigen::Matrix3d c_vec = eig.eigenvectors();
670 
671  Point3<ScalarType> lambda;
672  Matrix33<ScalarType> vect;
673  vect.FromEigenMatrix(c_vec);
674  lambda.FromEigenVector(c_val);
675 
676  ScalarType bestNormal = 0;
677  int bestNormalIndex = -1;
678  for(int i = 0; i < 3; ++i)
679  {
680  float agreeWithNormal = fabs((*vi).N().Normalize().dot(vect.GetColumn(i)));
681  if( agreeWithNormal > bestNormal )
682  {
683  bestNormal= agreeWithNormal;
684  bestNormalIndex = i;
685  }
686  }
687  int maxI = (bestNormalIndex+2)%3;
688  int minI = (bestNormalIndex+1)%3;
689  if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI);
690 
691  (*vi).PD1().Import(vect.GetColumn(maxI));
692  (*vi).PD2().Import(vect.GetColumn(minI));
693  (*vi).K1() = lambda[2];
694  (*vi).K2() = lambda[1];
695  }
696  }
697 
698  static void PerVertexBasicRadialCrossField(MeshType &m, float anisotropyRatio = 1.0 )
699  {
700  tri::RequirePerVertexCurvatureDir(m);
701  CoordType c=m.bbox.Center();
702  float maxRad = m.bbox.Diag()/2.0f;
703 
704  for(size_t i=0;i<m.vert.size();++i) {
705  CoordType dd = m.vert[i].P()-c;
706  dd.Normalize();
707  m.vert[i].PD1().Import(dd^m.vert[i].N());
708  m.vert[i].PD1().Normalize();
709  m.vert[i].PD2().Import(m.vert[i].N()^CoordType::Construct(m.vert[i].PD1()));
710  m.vert[i].PD2().Normalize();
711  // Now the anisotropy
712  // the idea is that the ratio between the two direction is at most <anisotropyRatio>
713  // eg |PD1|/|PD2| < ratio
714  // and |PD1|^2 + |PD2|^2 == 1
715 
716  float q =Distance(m.vert[i].P(),c) / maxRad; // it is in the 0..1 range
717  const float minRatio = 1.0f/anisotropyRatio;
718  const float maxRatio = anisotropyRatio;
719  const float curRatio = minRatio + (maxRatio-minRatio)*q;
720  float pd1Len = sqrt(1.0/(1+curRatio*curRatio));
721  float pd2Len = curRatio * pd1Len;
722 // assert(fabs(curRatio - pd2Len/pd1Len)<0.0000001);
723 // assert(fabs(pd1Len*pd1Len + pd2Len*pd2Len - 1.0f)<0.0001);
724  m.vert[i].PD1() *= pd1Len;
725  m.vert[i].PD2() *= pd2Len;
726  }
727  }
728 
729 };
730 
731 } // end namespace tri
732 } // end namespace vcg
733 #endif
vcg::Point3
Definition: point3.h:43
vcg::tri::Allocator
Class to safely add and delete elements in a mesh.
Definition: allocate.h:97
vcg::tri::Allocator::AddVertices
static VertexIterator AddVertices(MeshType &m, size_t n, PointerUpdater< VertexPointer > &pu)
Add n vertices to the mesh. Function to add n vertices to the mesh. The elements are added always to ...
Definition: allocate.h:189
vcg::tri::Inertia::Covariance
static void Covariance(const MeshType &m, vcg::Point3< ScalarType > &bary, vcg::Matrix33< ScalarType > &C)
Definition: inertia.h:318
vcg::tri::SurfaceSampling
Main Class of the Sampling framework.
Definition: point_sampling.h:467
vcg::tri::TrivialSampler
A basic sampler class that show the required interface used by the SurfaceSampling class.
Definition: point_sampling.h:72
vcg::tri::UpdateCurvature::AreaData
Definition: curvature.h:261
vcg::tri::UpdateCurvature
Management, updating and computation of per-vertex and per-face normals.
Definition: curvature.h:49
vcg::tri::UpdateCurvature::MeshGridType
vcg::GridStaticPtr< FaceType, ScalarType > MeshGridType
Definition: curvature.h:273
vcg::tri::UpdateCurvature::PerVertexAbsoluteMeanAndGaussian
static void PerVertexAbsoluteMeanAndGaussian(MeshType &m)
Update the mean and the gaussian curvature of a vertex.
Definition: curvature.h:536
vcg::tri::UpdateCurvature::PrincipalDirections
static void PrincipalDirections(MeshType &m)
Compute principal direction and magnitudo of curvature.
Definition: curvature.h:90
vcg::tri::UpdateCurvature::MeanAndGaussian
static void MeanAndGaussian(MeshType &m)
Computes the discrete mean gaussian curvature.
Definition: curvature.h:400
vcg::tri::UpdateNormal::NormalizePerVertex
static void NormalizePerVertex(ComputeMeshType &m)
Normalize the length of the vertex normals.
Definition: normal.h:239
vcg::tri::UpdateNormal::PerVertexAngleWeighted
static void PerVertexAngleWeighted(ComputeMeshType &m)
Calculates the vertex normal as an angle weighted average. It does not need or exploit current face n...
Definition: normal.h:124
vcg::tri::UpdateNormal::PerVertexNormalized
static void PerVertexNormalized(ComputeMeshType &m)
Equivalent to PerVertex() and NormalizePerVertex()
Definition: normal.h:269
vcg::face::IsBorder
bool IsBorder(FaceType const &f, const int j)
Definition: topology.h:55
vcg
Definition: namespaces.dox:6
Generated by doxygen 1.9.1 on Wed Mar 13 2024 13:53:31