00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055 #ifndef VCGLIB_UPDATE_CURVATURE_
00056 #define VCGLIB_UPDATE_CURVATURE_
00057
00058 #include <vcg/space/index/grid_static_ptr.h>
00059 #include <vcg/math/base.h>
00060 #include <vcg/math/matrix.h>
00061 #include <vcg/simplex/face/topology.h>
00062 #include <vcg/simplex/face/pos.h>
00063 #include <vcg/simplex/face/jumping_pos.h>
00064 #include <vcg/container/simple_temporary_data.h>
00065 #include <vcg/complex/trimesh/update/normal.h>
00066 #include <vcg/complex/trimesh/point_sampling.h>
00067 #include <vcg/complex/trimesh/append.h>
00068 #include <vcg/complex/intersection.h>
00069 #include <vcg/complex/trimesh/inertia.h>
00070 #include <vcg/math/matrix33.h>
00071
00072
00073 namespace vcg {
00074 namespace tri {
00075
00077
00079
00081
00085 template <class MeshType>
00086 class UpdateCurvature
00087 {
00088
00089 public:
00090 typedef typename MeshType::FaceType FaceType;
00091 typedef typename MeshType::FacePointer FacePointer;
00092 typedef typename MeshType::FaceIterator FaceIterator;
00093 typedef typename MeshType::VertexIterator VertexIterator;
00094 typedef typename MeshType::VertContainer VertContainer;
00095 typedef typename MeshType::VertexType VertexType;
00096 typedef typename MeshType::VertexPointer VertexPointer;
00097 typedef vcg::face::VFIterator<FaceType> VFIteratorType;
00098 typedef typename MeshType::CoordType CoordType;
00099 typedef typename CoordType::ScalarType ScalarType;
00100
00101
00102 private:
00103 typedef struct AdjVertex {
00104 VertexType * vert;
00105 float doubleArea;
00106 bool isBorder;
00107 };
00108
00109
00110 public:
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 static void PrincipalDirections(MeshType &m) {
00126
00127 assert(m.HasVFTopology());
00128
00129 vcg::tri::UpdateNormals<MeshType>::PerVertexNormalized(m);
00130
00131 VertexIterator vi;
00132 for (vi =m.vert.begin(); vi !=m.vert.end(); ++vi) {
00133 if ( ! (*vi).IsD() && (*vi).VFp() != NULL) {
00134
00135 VertexType * central_vertex = &(*vi);
00136
00137 std::vector<float> weights;
00138 std::vector<AdjVertex> vertices;
00139
00140 vcg::face::JumpingPos<FaceType> pos((*vi).VFp(), central_vertex);
00141
00142
00143 VertexType* firstV = pos.VFlip();
00144 VertexType* tempV;
00145 float totalDoubleAreaSize = 0.0f;
00146
00147
00148 do
00149 {
00150
00151 pos.FlipF();
00152 pos.FlipE();
00153
00154
00155 tempV = pos.VFlip();
00156 assert(tempV!=central_vertex);
00157 AdjVertex v;
00158
00159 v.isBorder = pos.IsBorder();
00160 v.vert = tempV;
00161 v.doubleArea = vcg::DoubleArea(*pos.F());
00162 totalDoubleAreaSize += v.doubleArea;
00163
00164 vertices.push_back(v);
00165 }
00166 while(tempV != firstV);
00167
00168
00169 for (size_t i = 0; i < vertices.size(); ++i) {
00170 if (vertices[i].isBorder) {
00171 weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize);
00172 } else {
00173 weights.push_back(0.5f * (vertices[i].doubleArea + vertices[(i-1)%vertices.size()].doubleArea) / totalDoubleAreaSize);
00174 }
00175 assert(weights.back() < 1.0f);
00176 }
00177
00178
00179 Matrix33<ScalarType> Tp;
00180 for (int i = 0; i < 3; ++i)
00181 Tp[i][i] = 1.0f - powf(central_vertex->cN()[i],2);
00182 Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->cN()[1]);
00183 Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->cN()[1] * central_vertex->cN()[2]);
00184 Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->cN()[0] * central_vertex->cN()[2]);
00185
00186
00187
00188 Matrix33<ScalarType> tempMatrix;
00189 Matrix33<ScalarType> M;
00190 M.SetZero();
00191 for (size_t i = 0; i < vertices.size(); ++i) {
00192 CoordType edge = (central_vertex->cP() - vertices[i].vert->cP());
00193 float curvature = (2.0f * (central_vertex->cN().dot(edge)) ) / edge.SquaredNorm();
00194 CoordType T = (Tp*edge).normalized();
00195 tempMatrix.ExternalProduct(T,T);
00196 M += tempMatrix * weights[i] * curvature ;
00197 }
00198
00199
00200 CoordType W;
00201 CoordType e1(1.0f,0.0f,0.0f);
00202 if ((e1 - central_vertex->cN()).SquaredNorm() > (e1 + central_vertex->cN()).SquaredNorm())
00203 W = e1 - central_vertex->cN();
00204 else
00205 W = e1 + central_vertex->cN();
00206 W.Normalize();
00207
00208
00209 Matrix33<ScalarType> Q;
00210 Q.SetIdentity();
00211 tempMatrix.ExternalProduct(W,W);
00212 Q -= tempMatrix * 2.0f;
00213
00214
00215 Matrix33<ScalarType> QtMQ = (Q.transpose() * M * Q);
00216
00217 CoordType bl = Q.GetColumn(0);
00218 CoordType T1 = Q.GetColumn(1);
00219 CoordType T2 = Q.GetColumn(2);
00220
00221
00222 float s,c;
00223
00224 float alpha = QtMQ[1][1]-QtMQ[2][2];
00225 float beta = QtMQ[2][1];
00226
00227 float h[2];
00228 float delta = sqrtf(4.0f*powf(alpha, 2) +16.0f*powf(beta, 2));
00229 h[0] = (2.0f*alpha + delta) / (2.0f*beta);
00230 h[1] = (2.0f*alpha - delta) / (2.0f*beta);
00231
00232 float t[2];
00233 float best_c, best_s;
00234 float min_error = std::numeric_limits<ScalarType>::infinity();
00235 for (int i=0; i<2; i++)
00236 {
00237 delta = sqrtf(powf(h[i], 2) + 4.0f);
00238 t[0] = (h[i]+delta) / 2.0f;
00239 t[1] = (h[i]-delta) / 2.0f;
00240
00241 for (int j=0; j<2; j++)
00242 {
00243 float squared_t = powf(t[j], 2);
00244 float denominator = 1.0f + squared_t;
00245 s = (2.0f*t[j]) / denominator;
00246 c = (1-squared_t) / denominator;
00247
00248 float approximation = c*s*alpha + (powf(c, 2) - powf(s, 2))*beta;
00249 float angle_similarity = fabs(acosf(c)/asinf(s));
00250 float error = fabs(1.0f-angle_similarity)+fabs(approximation);
00251 if (error<min_error)
00252 {
00253 min_error = error;
00254 best_c = c;
00255 best_s = s;
00256 }
00257 }
00258 }
00259 c = best_c;
00260 s = best_s;
00261
00262 vcg::ndim::MatrixMNf minor2x2 (2,2);
00263 vcg::ndim::MatrixMNf S (2,2);
00264
00265
00266
00267 minor2x2[0][0] = QtMQ[1][1];
00268 minor2x2[0][1] = QtMQ[1][2];
00269 minor2x2[1][0] = QtMQ[2][1];
00270 minor2x2[1][1] = QtMQ[2][2];
00271
00272 S[0][0] = S[1][1] = c;
00273 S[0][1] = s;
00274 S[1][0] = -1.0f * s;
00275
00276 vcg::ndim::MatrixMNf StMS(S.transpose() * minor2x2 * S);
00277
00278
00279 float Principal_Curvature1 = (3.0f * StMS[0][0]) - StMS[1][1];
00280 float Principal_Curvature2 = (3.0f * StMS[1][1]) - StMS[0][0];
00281
00282 CoordType Principal_Direction1 = T1 * c - T2 * s;
00283 CoordType Principal_Direction2 = T1 * s + T2 * c;
00284
00285 (*vi).PD1() = Principal_Direction1;
00286 (*vi).PD2() = Principal_Direction2;
00287 (*vi).K1() = Principal_Curvature1;
00288 (*vi).K2() = Principal_Curvature2;
00289 }
00290 }
00291 }
00292
00293
00294
00295
00296 class AreaData
00297 {
00298 public:
00299 float A;
00300 };
00301
00309 typedef vcg::GridStaticPtr <FaceType, ScalarType > MeshGridType;
00310 typedef vcg::GridStaticPtr <VertexType, ScalarType > PointsGridType;
00311
00312 static void PrincipalDirectionsPCA(MeshType &m, ScalarType r, bool pointVSfaceInt = true) {
00313 std::vector<VertexType*> closests;
00314 std::vector<ScalarType> distances;
00315 std::vector<CoordType> points;
00316 VertexIterator vi;
00317 ScalarType area;
00318 MeshType tmpM;
00319 typename std::vector<CoordType>::iterator ii;
00320 vcg::tri::TrivialSampler<MeshType> vs;
00321
00322 MeshGridType mGrid;
00323 PointsGridType pGrid;
00324
00325
00326 if(pointVSfaceInt){
00327 area = Stat<MeshType>::ComputeMeshArea(m);
00328 vcg::tri::SurfaceSampling<MeshType,vcg::tri::TrivialSampler<MeshType> >::Montecarlo(m,vs,1000 * area / (2*M_PI*r*r ));
00329 vi = vcg::tri::Allocator<MeshType>::AddVertices(tmpM,m.vert.size());
00330 for(size_t y = 0; y < m.vert.size(); ++y,++vi) (*vi).P() = m.vert[y].P();
00331 pGrid.Set(tmpM.vert.begin(),tmpM.vert.end());
00332 } else{ mGrid.Set(m.face.begin(),m.face.end()); }
00333
00334 for(vi = m.vert.begin(); vi != m.vert.end(); ++vi){
00335 vcg::Matrix33<ScalarType> A,eigenvectors;
00336 vcg::Point3<ScalarType> bp,eigenvalues;
00337 int nrot;
00338
00339
00340 if(pointVSfaceInt)
00341 {
00342 vcg::tri::GetInSphereVertex<
00343 MeshType,
00344 PointsGridType,std::vector<VertexType*>,
00345 std::vector<ScalarType>,
00346 std::vector<CoordType> >(tmpM,pGrid, (*vi).cP(),r ,closests,distances,points);
00347
00348 A.Covariance(points,bp);
00349 A*=area*area/1000;
00350 }
00351 else{
00352 IntersectionBallMesh<MeshType,ScalarType>( m ,vcg::Sphere3<ScalarType>((*vi).cP(),r),tmpM );
00353 vcg::Point3<ScalarType> _bary;
00354 vcg::tri::Inertia<MeshType>::Covariance(tmpM,_bary,A);
00355 }
00356
00357 Jacobi(A, eigenvalues , eigenvectors, nrot);
00358
00359
00360
00361 int best = 0; ScalarType bestv = fabs( (*vi).cN().dot(eigenvectors.GetColumn(0).normalized()) );
00362 for(int i = 1 ; i < 3; ++i){
00363 ScalarType prod = fabs((*vi).cN().dot(eigenvectors.GetColumn(i).normalized()));
00364 if( prod > bestv){bestv = prod; best = i;}
00365 }
00366
00367 (*vi).PD1() = eigenvectors.GetColumn( (best+1)%3).normalized();
00368 (*vi).PD2() = eigenvectors.GetColumn( (best+2)%3).normalized();
00369
00370
00371 vcg::Matrix33<ScalarType> rot;
00372 ScalarType angle = acos((*vi).PD1().dot((*vi).N()));
00373 rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD1()^(*vi).N());
00374 (*vi).PD1() = rot*(*vi).PD1();
00375 angle = acos((*vi).PD2().dot((*vi).N()));
00376 rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD2()^(*vi).N());
00377 (*vi).PD2() = rot*(*vi).PD2();
00378
00379
00380
00381 const ScalarType r5 = r*r*r*r*r;
00382 const ScalarType r6 = r*r5;
00383 (*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);
00384 (*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);
00385 if((*vi).K1() < (*vi).K2()) { std::swap((*vi).K1(),(*vi).K2());
00386 std::swap((*vi).PD1(),(*vi).PD2());
00387 }
00388 }
00389
00390
00391 }
00393
00399 static void MeanAndGaussian(MeshType & m)
00400 {
00401 assert(HasFFAdjacency(m));
00402 float area0, area1, area2, angle0, angle1, angle2;
00403 FaceIterator fi;
00404 VertexIterator vi;
00405 typename MeshType::CoordType e01v ,e12v ,e20v;
00406
00407 SimpleTempData<VertContainer, AreaData> TDAreaPtr(m.vert);
00408 SimpleTempData<VertContainer, typename MeshType::CoordType> TDContr(m.vert);
00409
00410 vcg::tri::UpdateNormals<MeshType>::PerVertexNormalized(m);
00411
00412 for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD())
00413 {
00414 (TDAreaPtr)[*vi].A = 0.0;
00415 (TDContr)[*vi] =typename MeshType::CoordType(0.0,0.0,0.0);
00416 (*vi).Kh() = 0.0;
00417 (*vi).Kg() = (float)(2.0 * M_PI);
00418 }
00419
00420 for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD())
00421 {
00422
00423 angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
00424 angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
00425 angle2 = M_PI-(angle0+angle1);
00426
00427 if((angle0 < M_PI/2) && (angle1 < M_PI/2) && (angle2 < M_PI/2))
00428 {
00429 float e01 = SquaredDistance( (*fi).V(1)->cP() , (*fi).V(0)->cP() );
00430 float e12 = SquaredDistance( (*fi).V(2)->cP() , (*fi).V(1)->cP() );
00431 float e20 = SquaredDistance( (*fi).V(0)->cP() , (*fi).V(2)->cP() );
00432
00433 area0 = ( e20*(1.0/tan(angle1)) + e01*(1.0/tan(angle2)) ) / 8.0;
00434 area1 = ( e01*(1.0/tan(angle2)) + e12*(1.0/tan(angle0)) ) / 8.0;
00435 area2 = ( e12*(1.0/tan(angle0)) + e20*(1.0/tan(angle1)) ) / 8.0;
00436
00437 (TDAreaPtr)[(*fi).V(0)].A += area0;
00438 (TDAreaPtr)[(*fi).V(1)].A += area1;
00439 (TDAreaPtr)[(*fi).V(2)].A += area2;
00440
00441 }
00442 else
00443 {
00444 if(angle0 >= M_PI/2)
00445 {
00446 (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
00447 (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00448 (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00449 }
00450 else if(angle1 >= M_PI/2)
00451 {
00452 (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00453 (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
00454 (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00455 }
00456 else
00457 {
00458 (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00459 (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 8.0;
00460 (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea<typename MeshType::FaceType>((*fi)) / 4.0;
00461 }
00462 }
00463 }
00464
00465
00466 for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD() )
00467 {
00468 angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) ));
00469 angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) ));
00470 angle2 = M_PI-(angle0+angle1);
00471
00472
00473 if(angle0==0 || angle1==0 || angle1==0) continue;
00474
00475 e01v = ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ;
00476 e12v = ( (*fi).V(2)->cP() - (*fi).V(1)->cP() ) ;
00477 e20v = ( (*fi).V(0)->cP() - (*fi).V(2)->cP() ) ;
00478
00479 TDContr[(*fi).V(0)] += ( e20v * (1.0/tan(angle1)) - e01v * (1.0/tan(angle2)) ) / 4.0;
00480 TDContr[(*fi).V(1)] += ( e01v * (1.0/tan(angle2)) - e12v * (1.0/tan(angle0)) ) / 4.0;
00481 TDContr[(*fi).V(2)] += ( e12v * (1.0/tan(angle0)) - e20v * (1.0/tan(angle1)) ) / 4.0;
00482
00483 (*fi).V(0)->Kg() -= angle0;
00484 (*fi).V(1)->Kg() -= angle1;
00485 (*fi).V(2)->Kg() -= angle2;
00486
00487
00488 for(int i=0;i<3;i++)
00489 {
00490 if(vcg::face::IsBorder((*fi), i))
00491 {
00492 CoordType e1,e2;
00493 vcg::face::Pos<FaceType> hp(&*fi, i, (*fi).V(i));
00494 vcg::face::Pos<FaceType> hp1=hp;
00495
00496 hp1.FlipV();
00497 e1=hp1.v->cP() - hp.v->cP();
00498 hp1.FlipV();
00499 hp1.NextB();
00500 e2=hp1.v->cP() - hp.v->cP();
00501 (*fi).V(i)->Kg() -= math::Abs(Angle(e1,e2));
00502 }
00503 }
00504 }
00505
00506 for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD() )
00507 {
00508 if((TDAreaPtr)[*vi].A<=std::numeric_limits<ScalarType>::epsilon())
00509 {
00510 (*vi).Kh() = 0;
00511 (*vi).Kg() = 0;
00512 }
00513 else
00514 {
00515 (*vi).Kh() = (((TDContr)[*vi].dot((*vi).cN())>0)?1.0:-1.0)*((TDContr)[*vi] / (TDAreaPtr) [*vi].A).Norm();
00516 (*vi).Kg() /= (TDAreaPtr)[*vi].A;
00517 }
00518 }
00519 }
00520
00521
00523
00530 static float VertexCurvature(VertexPointer v, bool norm = true)
00531 {
00532
00533 assert(FaceType::HasVFAdjacency());
00534 assert(VertexType::HasVFAdjacency());
00535
00536 VFIteratorType vfi(v);
00537 float A = 0;
00538
00539 v->Kh() = 0;
00540 v->Kg() = 2 * M_PI;
00541
00542 while (!vfi.End()) {
00543 if (!vfi.F()->IsD()) {
00544 FacePointer f = vfi.F();
00545 int i = vfi.I();
00546 VertexPointer v0 = f->V0(i), v1 = f->V1(i), v2 = f->V2(i);
00547
00548 float ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P() ));
00549 float ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P() ));
00550 float ang2 = M_PI - ang0 - ang1;
00551
00552 float s01 = SquaredDistance(v1->P(), v0->P());
00553 float s02 = SquaredDistance(v2->P(), v0->P());
00554
00555
00556 if (ang0 >= M_PI/2)
00557 A += (0.5f * DoubleArea(*f) - (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 );
00558 else if (ang1 >= M_PI/2)
00559 A += (s01 * tan(ang0)) / 8.0;
00560 else if (ang2 >= M_PI/2)
00561 A += (s02 * tan(ang0)) / 8.0;
00562 else
00563 A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0;
00564
00565
00566 v->Kg() -= ang0;
00567
00568
00569 ang1 = math::Abs(Angle(f->N(), v1->N()));
00570 ang2 = math::Abs(Angle(f->N(), v2->N()));
00571 v->Kh() += ( (math::Sqrt(s01) / 2.0) * ang1 +
00572 (math::Sqrt(s02) / 2.0) * ang2 );
00573 }
00574
00575 ++vfi;
00576 }
00577
00578 v->Kh() /= 4.0f;
00579
00580 if(norm) {
00581 if(A <= std::numeric_limits<float>::epsilon()) {
00582 v->Kh() = 0;
00583 v->Kg() = 0;
00584 }
00585 else {
00586 v->Kh() /= A;
00587 v->Kg() /= A;
00588 }
00589 }
00590
00591 return A;
00592 }
00593
00594 static void VertexCurvature(MeshType & m){
00595
00596 for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
00597 VertexCurvature(&*vi,false);
00598 }
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612 static void PrincipalDirectionsNormalCycles(MeshType & m){
00613 assert(VertexType::HasVFAdjacency());
00614 assert(FaceType::HasFFAdjacency());
00615 assert(FaceType::HasFaceNormal());
00616
00617
00618 typename MeshType::VertexIterator vi;
00619
00620 for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
00621 if(!((*vi).IsD())){
00622 vcg::Matrix33<ScalarType> m33;m33.SetZero();
00623 face::JumpingPos<typename MeshType::FaceType> p((*vi).VFp(),&(*vi));
00624 p.FlipE();
00625 typename MeshType::VertexType * firstv = p.VFlip();
00626 assert(p.F()->V(p.VInd())==&(*vi));
00627
00628
00629 do{
00630 if( p.F() != p.FFlip()){
00631 Point3<ScalarType> normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P();
00632 ScalarType edge_length = normalized_edge.Norm();
00633 normalized_edge/=edge_length;
00634 Point3<ScalarType> n1 = p.F()->cN();n1.Normalize();
00635 Point3<ScalarType> n2 = p.FFlip()->cN();n2.Normalize();
00636 ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge);
00637 n1n2 = math::Max<ScalarType >(math::Min<ScalarType> ( 1.0,n1n2),-1.0);
00638 ScalarType beta = math::Asin(n1n2);
00639 m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0];
00640 m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0];
00641 m33[1][1] += beta*edge_length*normalized_edge[1]*normalized_edge[1];
00642 m33[0][2] += beta*edge_length*normalized_edge[2]*normalized_edge[0];
00643 m33[1][2] += beta*edge_length*normalized_edge[2]*normalized_edge[1];
00644 m33[2][2] += beta*edge_length*normalized_edge[2]*normalized_edge[2];
00645 }
00646 p.NextFE();
00647 }while(firstv != p.VFlip());
00648
00649 if(m33.Determinant()==0.0){
00650 (*vi).K1() = (*vi).K2() = 0.0; continue;}
00651
00652 m33[1][0] = m33[0][1];
00653 m33[2][0] = m33[0][2];
00654 m33[2][1] = m33[1][2];
00655
00656 Point3<ScalarType> lambda;
00657 Matrix33<ScalarType> vect;
00658 int n_rot;
00659 Jacobi(m33,lambda, vect,n_rot);
00660
00661 vect.transposeInPlace();
00662 ScalarType normal = std::numeric_limits<ScalarType>::min();
00663 int normI = 0;
00664 for(int i = 0; i < 3; ++i)
00665 if( fabs((*vi).N().Normalize().dot(vect.GetRow(i))) > normal )
00666 {
00667 normal= fabs((*vi).N().Normalize().dot(vect.GetRow(i)));
00668 normI = i;
00669 }
00670 int maxI = (normI+2)%3;
00671 int minI = (normI+1)%3;
00672 if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI);
00673
00674 (*vi).PD1() = *(Point3<ScalarType>*)(& vect[maxI][0]);
00675 (*vi).PD2() = *(Point3<ScalarType>*)(& vect[minI][0]);
00676 (*vi).K1() = lambda[maxI];
00677 (*vi).K2() = lambda[minI];
00678 }
00679 }
00680 };
00681
00682 }
00683 }
00684 #endif
