VCG Library
Loading...
Searching...
No Matches
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
35namespace vcg {
36namespace tri {
37
39
41
43
47template <class MeshType>
48class UpdateCurvature
49{
50
51public:
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
66private:
67 // aux data struct used by PrincipalDirections
68 struct AdjVertex {
69 VertexType * vert;
70 float doubleArea;
71 bool isBorder;
72 };
73
74
75public:
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
400static 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
536static 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:476
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: color4.h:30
Generated by doxygen 1.9.6 on Mon Dec 1 2025 17:19:44