diff --git a/vcg/complex/algorithms/update/curvature.h b/vcg/complex/algorithms/update/curvature.h index b842011e..98c5faee 100644 --- a/vcg/complex/algorithms/update/curvature.h +++ b/vcg/complex/algorithms/update/curvature.h @@ -19,55 +19,20 @@ * GNU General Public License (http://www.gnu.org/licenses/gpl.txt) * * for more details. * * * -****************************************************************************/ -/**************************************************************************** - History -$Log: not supported by cvs2svn $ -Revision 1.8 2008/05/14 10:03:29 ganovelli -Point3f->Coordtype - -Revision 1.7 2008/04/23 16:37:15 onnis -VertexCurvature method added. - -Revision 1.6 2008/04/04 10:26:12 cignoni -Cleaned up names, now Kg() gives back Gaussian Curvature (k1*k2), while Kh() gives back Mean Curvature 1/2(k1+k2) - -Revision 1.5 2008/03/25 11:00:56 ganovelli -fixed bugs sign of principal direction and mean curvature value - -Revision 1.4 2008/03/17 11:29:59 ganovelli -taubin and desbrun estimates added (-> see vcg/simplex/vertex/component.h [component_ocf.h|component_occ.h ] - -Revision 1.3 2006/02/27 18:02:57 ponchio -Area -> doublearea/2 - -added some typename - -Revision 1.2 2005/10/25 09:17:41 spinelli -correct IsBorder - -Revision 1.1 2005/02/22 16:40:29 ganovelli -created. This version writes the gaussian curvature on the Q() member of -the vertex - ****************************************************************************/ #ifndef VCGLIB_UPDATE_CURVATURE_ #define VCGLIB_UPDATE_CURVATURE_ #include -#include -#include #include #include #include -#include #include #include -#include #include #include -#include +#include #include namespace vcg { @@ -100,11 +65,12 @@ public: private: - struct AdjVertex { - VertexType * vert; - float doubleArea; - bool isBorder; - }; + // aux data struct used by PrincipalDirections + struct AdjVertex { + VertexType * vert; + float doubleArea; + bool isBorder; + }; public: @@ -122,173 +88,174 @@ public: URL = "http://dx.doi.org/10.1109/ICCV.1995.466840", bibsource = "http://www.visionbib.com/bibliography/describe440.html#TT32253", */ - static void PrincipalDirections(MeshType &m) { + static void PrincipalDirections(MeshType &m) + { - assert(tri::HasPerFaceVFAdjacency(m) && tri::HasPerVertexVFAdjacency(m)); + assert(tri::HasPerFaceVFAdjacency(m) && tri::HasPerVertexVFAdjacency(m)); - vcg::tri::UpdateNormal::PerVertexNormalized(m); + vcg::tri::UpdateNormal::PerVertexAngleWeighted(m); + vcg::tri::UpdateNormal::NormalizePerVertex(m); - VertexIterator vi; - for (vi =m.vert.begin(); vi !=m.vert.end(); ++vi) { - if ( ! (*vi).IsD() && (*vi).VFp() != NULL) { + for (VertexIterator vi =m.vert.begin(); vi !=m.vert.end(); ++vi) { + if ( ! (*vi).IsD() && (*vi).VFp() != NULL) { - VertexType * central_vertex = &(*vi); + VertexType * central_vertex = &(*vi); - std::vector weights; - std::vector vertices; + std::vector weights; + std::vector vertices; - vcg::face::JumpingPos pos((*vi).VFp(), central_vertex); + vcg::face::JumpingPos pos((*vi).VFp(), central_vertex); - // firstV is the first vertex of the 1ring neighboorhood - VertexType* firstV = pos.VFlip(); - VertexType* tempV; - float totalDoubleAreaSize = 0.0f; + // firstV is the first vertex of the 1ring neighboorhood + VertexType* firstV = pos.VFlip(); + VertexType* tempV; + float totalDoubleAreaSize = 0.0f; - // compute the area of each triangle around the central vertex as well as their total area - do - { - // this bring the pos to the next triangle counterclock-wise - pos.FlipF(); - pos.FlipE(); + // compute the area of each triangle around the central vertex as well as their total area + do + { + // this bring the pos to the next triangle counterclock-wise + pos.FlipF(); + pos.FlipE(); - // tempV takes the next vertex in the 1ring neighborhood - tempV = pos.VFlip(); - assert(tempV!=central_vertex); - AdjVertex v; + // tempV takes the next vertex in the 1ring neighborhood + tempV = pos.VFlip(); + assert(tempV!=central_vertex); + AdjVertex v; - v.isBorder = pos.IsBorder(); - v.vert = tempV; - v.doubleArea = vcg::DoubleArea(*pos.F()); - totalDoubleAreaSize += v.doubleArea; + v.isBorder = pos.IsBorder(); + v.vert = tempV; + v.doubleArea = vcg::DoubleArea(*pos.F()); + totalDoubleAreaSize += v.doubleArea; - vertices.push_back(v); - } - while(tempV != firstV); + vertices.push_back(v); + } + while(tempV != firstV); - // compute the weights for the formula computing matrix M - for (size_t i = 0; i < vertices.size(); ++i) { - if (vertices[i].isBorder) { - weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize); - } else { - weights.push_back(0.5f * (vertices[i].doubleArea + vertices[(i-1)%vertices.size()].doubleArea) / totalDoubleAreaSize); - } - assert(weights.back() < 1.0f); - } + // compute the weights for the formula computing matrix M + for (size_t i = 0; i < vertices.size(); ++i) { + if (vertices[i].isBorder) { + weights.push_back(vertices[i].doubleArea / totalDoubleAreaSize); + } else { + weights.push_back(0.5f * (vertices[i].doubleArea + vertices[(i-1)%vertices.size()].doubleArea) / totalDoubleAreaSize); + } + assert(weights.back() < 1.0f); + } - // compute I-NN^t to be used for computing the T_i's - Matrix33 Tp; - for (int i = 0; i < 3; ++i) - Tp[i][i] = 1.0f - powf(central_vertex->cN()[i],2); - Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->cN()[1]); - Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->cN()[1] * central_vertex->cN()[2]); - Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->cN()[0] * central_vertex->cN()[2]); + // compute I-NN^t to be used for computing the T_i's + Matrix33 Tp; + for (int i = 0; i < 3; ++i) + Tp[i][i] = 1.0f - powf(central_vertex->cN()[i],2); + Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->cN()[1]); + Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->cN()[1] * central_vertex->cN()[2]); + Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->cN()[0] * central_vertex->cN()[2]); - // for all neighbors vi compute the directional curvatures k_i and the T_i - // compute M by summing all w_i k_i T_i T_i^t - Matrix33 tempMatrix; - Matrix33 M; - M.SetZero(); - for (size_t i = 0; i < vertices.size(); ++i) { - CoordType edge = (central_vertex->cP() - vertices[i].vert->cP()); - float curvature = (2.0f * (central_vertex->cN().dot(edge)) ) / edge.SquaredNorm(); - CoordType T = (Tp*edge).normalized(); - tempMatrix.ExternalProduct(T,T); - M += tempMatrix * weights[i] * curvature ; - } + // for all neighbors vi compute the directional curvatures k_i and the T_i + // compute M by summing all w_i k_i T_i T_i^t + Matrix33 tempMatrix; + Matrix33 M; + M.SetZero(); + for (size_t i = 0; i < vertices.size(); ++i) { + CoordType edge = (central_vertex->cP() - vertices[i].vert->cP()); + float curvature = (2.0f * (central_vertex->cN().dot(edge)) ) / edge.SquaredNorm(); + CoordType T = (Tp*edge).normalized(); + tempMatrix.ExternalProduct(T,T); + M += tempMatrix * weights[i] * curvature ; + } - // compute vector W for the Householder matrix - CoordType W; - CoordType e1(1.0f,0.0f,0.0f); - if ((e1 - central_vertex->cN()).SquaredNorm() > (e1 + central_vertex->cN()).SquaredNorm()) - W = e1 - central_vertex->cN(); - else - W = e1 + central_vertex->cN(); - W.Normalize(); + // compute vector W for the Householder matrix + CoordType W; + CoordType e1(1.0f,0.0f,0.0f); + if ((e1 - central_vertex->cN()).SquaredNorm() > (e1 + central_vertex->cN()).SquaredNorm()) + W = e1 - central_vertex->cN(); + else + W = e1 + central_vertex->cN(); + W.Normalize(); - // compute the Householder matrix I - 2WW^t - Matrix33 Q; - Q.SetIdentity(); - tempMatrix.ExternalProduct(W,W); - Q -= tempMatrix * 2.0f; + // compute the Householder matrix I - 2WW^t + Matrix33 Q; + Q.SetIdentity(); + tempMatrix.ExternalProduct(W,W); + Q -= tempMatrix * 2.0f; - // compute matrix Q^t M Q - Matrix33 QtMQ = (Q.transpose() * M * Q); + // compute matrix Q^t M Q + Matrix33 QtMQ = (Q.transpose() * M * Q); - CoordType bl = Q.GetColumn(0); - CoordType T1 = Q.GetColumn(1); - CoordType T2 = Q.GetColumn(2); + CoordType bl = Q.GetColumn(0); + CoordType T1 = Q.GetColumn(1); + CoordType T2 = Q.GetColumn(2); - // find sin and cos for the Givens rotation - float s,c; - // Gabriel Taubin hint and Valentino Fiorin impementation - float alpha = QtMQ[1][1]-QtMQ[2][2]; - float beta = QtMQ[2][1]; + // find sin and cos for the Givens rotation + float s,c; + // Gabriel Taubin hint and Valentino Fiorin impementation + float alpha = QtMQ[1][1]-QtMQ[2][2]; + float beta = QtMQ[2][1]; - float h[2]; - float delta = sqrtf(4.0f*powf(alpha, 2) +16.0f*powf(beta, 2)); - h[0] = (2.0f*alpha + delta) / (2.0f*beta); - h[1] = (2.0f*alpha - delta) / (2.0f*beta); + float h[2]; + float delta = sqrtf(4.0f*powf(alpha, 2) +16.0f*powf(beta, 2)); + h[0] = (2.0f*alpha + delta) / (2.0f*beta); + h[1] = (2.0f*alpha - delta) / (2.0f*beta); - float t[2]; - float best_c, best_s; - float min_error = std::numeric_limits::infinity(); - for (int i=0; i<2; i++) - { - delta = sqrtf(powf(h[i], 2) + 4.0f); - t[0] = (h[i]+delta) / 2.0f; - t[1] = (h[i]-delta) / 2.0f; + float t[2]; + float best_c, best_s; + float min_error = std::numeric_limits::infinity(); + for (int i=0; i<2; i++) + { + delta = sqrtf(powf(h[i], 2) + 4.0f); + t[0] = (h[i]+delta) / 2.0f; + t[1] = (h[i]-delta) / 2.0f; - for (int j=0; j<2; j++) - { - float squared_t = powf(t[j], 2); - float denominator = 1.0f + squared_t; - s = (2.0f*t[j]) / denominator; - c = (1-squared_t) / denominator; + for (int j=0; j<2; j++) + { + float squared_t = powf(t[j], 2); + float denominator = 1.0f + squared_t; + s = (2.0f*t[j]) / denominator; + c = (1-squared_t) / denominator; - float approximation = c*s*alpha + (powf(c, 2) - powf(s, 2))*beta; - float angle_similarity = fabs(acosf(c)/asinf(s)); - float error = fabs(1.0f-angle_similarity)+fabs(approximation); - if (error MeshGridType; - typedef vcg::GridStaticPtr PointsGridType; - - static void PrincipalDirectionsPCA(MeshType &m, ScalarType r, bool pointVSfaceInt = true,vcg::CallBackPos * cb = NULL) { - std::vector closests; - std::vector distances; - std::vector points; - VertexIterator vi; - ScalarType area; - MeshType tmpM; - typename std::vector::iterator ii; - vcg::tri::TrivialSampler vs; - - MeshGridType mGrid; - PointsGridType pGrid; - - // Fill the grid used - if(pointVSfaceInt){ - area = Stat::ComputeMeshArea(m); - vcg::tri::SurfaceSampling >::Montecarlo(m,vs,1000 * area / (2*M_PI*r*r )); - vi = vcg::tri::Allocator::AddVertices(tmpM,m.vert.size()); - for(size_t y = 0; y < m.vert.size(); ++y,++vi) (*vi).P() = m.vert[y].P(); - pGrid.Set(tmpM.vert.begin(),tmpM.vert.end()); - } else{ mGrid.Set(m.face.begin(),m.face.end()); } - int jj = 0; - for(vi = m.vert.begin(); vi != m.vert.end(); ++vi){ - vcg::Matrix33 A,eigenvectors; - vcg::Point3 bp,eigenvalues; - int nrot; - - // sample the neighborhood - if(pointVSfaceInt) - { - vcg::tri::GetInSphereVertex< - MeshType, - PointsGridType,std::vector, - std::vector, - std::vector >(tmpM,pGrid, (*vi).cP(),r ,closests,distances,points); - - A.Covariance(points,bp); - A*=area*area/1000; - } - else{ - IntersectionBallMesh( m ,vcg::Sphere3((*vi).cP(),r),tmpM ); - vcg::Point3 _bary; - vcg::tri::Inertia::Covariance(tmpM,_bary,A); - } - - Jacobi(A, eigenvalues , eigenvectors, nrot); - - // get the estimate of curvatures from eigenvalues and eigenvectors - // find the 2 most tangent eigenvectors (by finding the one closest to the normal) - int best = 0; ScalarType bestv = fabs( (*vi).cN().dot(eigenvectors.GetColumn(0).normalized()) ); - for(int i = 1 ; i < 3; ++i){ - ScalarType prod = fabs((*vi).cN().dot(eigenvectors.GetColumn(i).normalized())); - if( prod > bestv){bestv = prod; best = i;} - } - - (*vi).PD1() = eigenvectors.GetColumn( (best+1)%3).normalized(); - (*vi).PD2() = eigenvectors.GetColumn( (best+2)%3).normalized(); - - // project them to the plane identified by the normal - vcg::Matrix33 rot; - ScalarType angle = acos((*vi).PD1().dot((*vi).N())); - rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD1()^(*vi).N()); - (*vi).PD1() = rot*(*vi).PD1(); - angle = acos((*vi).PD2().dot((*vi).N())); - rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD2()^(*vi).N()); - (*vi).PD2() = rot*(*vi).PD2(); - - - // copmutes the curvature values - const ScalarType r5 = r*r*r*r*r; - const ScalarType r6 = r*r5; - (*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); - (*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); - if((*vi).K1() < (*vi).K2()) { std::swap((*vi).K1(),(*vi).K2()); - std::swap((*vi).PD1(),(*vi).PD2()); - if (cb) - { - (*cb)(int(100.0f * (float)jj / (float)m.vn),"Vertices Analysis"); - ++jj; - } } - } - - - } - /// \brief Computes the discrete gaussian curvature. - -/** For further details, please, refer to: \n - -- Discrete Differential-Geometry Operators for Triangulated 2-Manifolds Mark Meyer, - Mathieu Desbrun, Peter Schroder, Alan H. Barr VisMath '02, Berlin + /** Curvature meseaure as described in the paper: +Robust principal curvatures on Multiple Scales, Yong-Liang Yang, Yu-Kun Lai, Shi-Min Hu Helmut Pottmann +SGP 2004 +If pointVSfaceInt==true the covariance is computed by montecarlo sampling on the mesh (faster) +If pointVSfaceInt==false the covariance is computed by (analytic)integration over the surface (slower) */ - static void MeanAndGaussian(MeshType & m) + + typedef vcg::GridStaticPtr MeshGridType; + typedef vcg::GridStaticPtr PointsGridType; + + static void PrincipalDirectionsPCA(MeshType &m, ScalarType r, bool pointVSfaceInt = true,vcg::CallBackPos * cb = NULL) + { + std::vector closests; + std::vector distances; + std::vector points; + VertexIterator vi; + ScalarType area; + MeshType tmpM; + typename std::vector::iterator ii; + vcg::tri::TrivialSampler vs; + tri::UpdateNormal::PerVertexAngleWeighted(m); + tri::UpdateNormal::NormalizePerVertex(m); + + MeshGridType mGrid; + PointsGridType pGrid; + + // Fill the grid used + if(pointVSfaceInt) { - assert(HasFFAdjacency(m)); - float area0, area1, area2, angle0, angle1, angle2; - FaceIterator fi; - VertexIterator vi; - typename MeshType::CoordType e01v ,e12v ,e20v; - - SimpleTempData TDAreaPtr(m.vert); - SimpleTempData TDContr(m.vert); - - vcg::tri::UpdateNormal::PerVertexNormalized(m); - //Compute AreaMix in H (vale anche per K) - for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD()) - { - (TDAreaPtr)[*vi].A = 0.0; - (TDContr)[*vi] =typename MeshType::CoordType(0.0,0.0,0.0); - (*vi).Kh() = 0.0; - (*vi).Kg() = (float)(2.0 * M_PI); - } - - for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD()) - { - // angles - angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) )); - angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) )); - angle2 = M_PI-(angle0+angle1); - - if((angle0 < M_PI/2) && (angle1 < M_PI/2) && (angle2 < M_PI/2)) // triangolo non ottuso - { - float e01 = SquaredDistance( (*fi).V(1)->cP() , (*fi).V(0)->cP() ); - float e12 = SquaredDistance( (*fi).V(2)->cP() , (*fi).V(1)->cP() ); - float e20 = SquaredDistance( (*fi).V(0)->cP() , (*fi).V(2)->cP() ); - - area0 = ( e20*(1.0/tan(angle1)) + e01*(1.0/tan(angle2)) ) / 8.0; - area1 = ( e01*(1.0/tan(angle2)) + e12*(1.0/tan(angle0)) ) / 8.0; - area2 = ( e12*(1.0/tan(angle0)) + e20*(1.0/tan(angle1)) ) / 8.0; - - (TDAreaPtr)[(*fi).V(0)].A += area0; - (TDAreaPtr)[(*fi).V(1)].A += area1; - (TDAreaPtr)[(*fi).V(2)].A += area2; - - } - else // obtuse - { - if(angle0 >= M_PI/2) - { - (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 4.0; - (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 8.0; - (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 8.0; - } - else if(angle1 >= M_PI/2) - { - (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 8.0; - (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 4.0; - (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 8.0; - } - else - { - (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 8.0; - (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 8.0; - (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 4.0; - } - } - } - - - for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD() ) - { - angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) )); - angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) )); - angle2 = M_PI-(angle0+angle1); - - // Skip degenerate triangles. - if(angle0==0 || angle1==0 || angle1==0) continue; - - e01v = ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ; - e12v = ( (*fi).V(2)->cP() - (*fi).V(1)->cP() ) ; - e20v = ( (*fi).V(0)->cP() - (*fi).V(2)->cP() ) ; - - TDContr[(*fi).V(0)] += ( e20v * (1.0/tan(angle1)) - e01v * (1.0/tan(angle2)) ) / 4.0; - TDContr[(*fi).V(1)] += ( e01v * (1.0/tan(angle2)) - e12v * (1.0/tan(angle0)) ) / 4.0; - TDContr[(*fi).V(2)] += ( e12v * (1.0/tan(angle0)) - e20v * (1.0/tan(angle1)) ) / 4.0; - - (*fi).V(0)->Kg() -= angle0; - (*fi).V(1)->Kg() -= angle1; - (*fi).V(2)->Kg() -= angle2; - - - for(int i=0;i<3;i++) - { - if(vcg::face::IsBorder((*fi), i)) - { - CoordType e1,e2; - vcg::face::Pos hp(&*fi, i, (*fi).V(i)); - vcg::face::Pos hp1=hp; - - hp1.FlipV(); - e1=hp1.v->cP() - hp.v->cP(); - hp1.FlipV(); - hp1.NextB(); - e2=hp1.v->cP() - hp.v->cP(); - (*fi).V(i)->Kg() -= math::Abs(Angle(e1,e2)); - } - } - } - - for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD() /*&& !(*vi).IsB()*/) - { - if((TDAreaPtr)[*vi].A<=std::numeric_limits::epsilon()) - { - (*vi).Kh() = 0; - (*vi).Kg() = 0; - } - else - { - (*vi).Kh() = (((TDContr)[*vi].dot((*vi).cN())>0)?1.0:-1.0)*((TDContr)[*vi] / (TDAreaPtr) [*vi].A).Norm(); - (*vi).Kg() /= (TDAreaPtr)[*vi].A; - } - } + area = Stat::ComputeMeshArea(m); + vcg::tri::SurfaceSampling >::Montecarlo(m,vs,1000 * area / (2*M_PI*r*r )); + vi = vcg::tri::Allocator::AddVertices(tmpM,m.vert.size()); + for(size_t y = 0; y < m.vert.size(); ++y,++vi) (*vi).P() = m.vert[y].P(); + pGrid.Set(tmpM.vert.begin(),tmpM.vert.end()); } + else + { + mGrid.Set(m.face.begin(),m.face.end()); + } + + int jj = 0; + for(vi = m.vert.begin(); vi != m.vert.end(); ++vi) + { + vcg::Matrix33 A, eigenvectors; + vcg::Point3 bp, eigenvalues; +// int nrot; + + // sample the neighborhood + if(pointVSfaceInt) + { + vcg::tri::GetInSphereVertex< + MeshType, + PointsGridType,std::vector, + std::vector, + std::vector >(tmpM,pGrid, (*vi).cP(),r ,closests,distances,points); + + A.Covariance(points,bp); + A*=area*area/1000; + } + else{ + IntersectionBallMesh( m ,vcg::Sphere3((*vi).cP(),r),tmpM ); + vcg::Point3 _bary; + vcg::tri::Inertia::Covariance(tmpM,_bary,A); + } + +// Eigen::Matrix3f AA; AA=A; +// Eigen::SelfAdjointEigenSolver eig(AA); +// Eigen::Vector3f c_val = eig.eigenvalues(); +// Eigen::Matrix3f c_vec = eig.eigenvectors(); + +// Jacobi(A, eigenvalues , eigenvectors, nrot); + + Eigen::Matrix3d AA; + A.ToEigenMatrix(AA); + Eigen::SelfAdjointEigenSolver eig(AA); + Eigen::Vector3d c_val = eig.eigenvalues(); + Eigen::Matrix3d c_vec = eig.eigenvectors(); // eigenvector are stored as columns. + eigenvectors.FromEigenMatrix(c_vec); + eigenvalues.FromEigenVector(c_val); +// EV.transposeInPlace(); +// ev.FromEigenVector(c_val); + + // get the estimate of curvatures from eigenvalues and eigenvectors + // find the 2 most tangent eigenvectors (by finding the one closest to the normal) + int best = 0; ScalarType bestv = fabs( (*vi).cN().dot(eigenvectors.GetColumn(0).normalized()) ); + for(int i = 1 ; i < 3; ++i){ + ScalarType prod = fabs((*vi).cN().dot(eigenvectors.GetColumn(i).normalized())); + if( prod > bestv){bestv = prod; best = i;} + } + + (*vi).PD1() = eigenvectors.GetColumn( (best+1)%3).normalized(); + (*vi).PD2() = eigenvectors.GetColumn( (best+2)%3).normalized(); + + // project them to the plane identified by the normal + vcg::Matrix33 rot; + ScalarType angle = acos((*vi).PD1().dot((*vi).N())); + rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD1()^(*vi).N()); + (*vi).PD1() = rot*(*vi).PD1(); + angle = acos((*vi).PD2().dot((*vi).N())); + rot.SetRotateRad( - (M_PI*0.5 - angle),(*vi).PD2()^(*vi).N()); + (*vi).PD2() = rot*(*vi).PD2(); + + + // copmutes the curvature values + const ScalarType r5 = r*r*r*r*r; + const ScalarType r6 = r*r5; + (*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); + (*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); + if((*vi).K1() < (*vi).K2()) { std::swap((*vi).K1(),(*vi).K2()); + std::swap((*vi).PD1(),(*vi).PD2()); + if (cb) + { + (*cb)(int(100.0f * (float)jj / (float)m.vn),"Vertices Analysis"); + ++jj; + } } + } + + + } + +/// \brief Computes the discrete mean gaussian curvature. +/** +The algorithm used is the one Desbrun et al. that is based on a discrete analysis of the angles of the faces around a vertex. + +It requires FaceFace Adjacency; + +For further details, please, refer to: \n + +Discrete Differential-Geometry Operators for Triangulated 2-Manifolds
+Mark Meyer, Mathieu Desbrun, Peter Schroder, Alan H. Barr
+ VisMath '02, Berlin +*/ +static void MeanAndGaussian(MeshType & m) +{ + assert(HasFFAdjacency(m)); + float area0, area1, area2, angle0, angle1, angle2; + FaceIterator fi; + VertexIterator vi; + typename MeshType::CoordType e01v ,e12v ,e20v; + + SimpleTempData TDAreaPtr(m.vert); + SimpleTempData TDContr(m.vert); + + vcg::tri::UpdateNormal::PerVertexNormalized(m); + //Compute AreaMix in H (vale anche per K) + for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD()) + { + (TDAreaPtr)[*vi].A = 0.0; + (TDContr)[*vi] =typename MeshType::CoordType(0.0,0.0,0.0); + (*vi).Kh() = 0.0; + (*vi).Kg() = (float)(2.0 * M_PI); + } + + for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD()) + { + // angles + angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) )); + angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) )); + angle2 = M_PI-(angle0+angle1); + + if((angle0 < M_PI/2) && (angle1 < M_PI/2) && (angle2 < M_PI/2)) // triangolo non ottuso + { + float e01 = SquaredDistance( (*fi).V(1)->cP() , (*fi).V(0)->cP() ); + float e12 = SquaredDistance( (*fi).V(2)->cP() , (*fi).V(1)->cP() ); + float e20 = SquaredDistance( (*fi).V(0)->cP() , (*fi).V(2)->cP() ); + + area0 = ( e20*(1.0/tan(angle1)) + e01*(1.0/tan(angle2)) ) / 8.0; + area1 = ( e01*(1.0/tan(angle2)) + e12*(1.0/tan(angle0)) ) / 8.0; + area2 = ( e12*(1.0/tan(angle0)) + e20*(1.0/tan(angle1)) ) / 8.0; + + (TDAreaPtr)[(*fi).V(0)].A += area0; + (TDAreaPtr)[(*fi).V(1)].A += area1; + (TDAreaPtr)[(*fi).V(2)].A += area2; + + } + else // obtuse + { + if(angle0 >= M_PI/2) + { + (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 4.0; + (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 8.0; + (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 8.0; + } + else if(angle1 >= M_PI/2) + { + (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 8.0; + (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 4.0; + (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 8.0; + } + else + { + (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 8.0; + (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 8.0; + (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 4.0; + } + } + } + + + for(fi=m.face.begin();fi!=m.face.end();++fi) if( !(*fi).IsD() ) + { + angle0 = math::Abs(Angle( (*fi).P(1)-(*fi).P(0),(*fi).P(2)-(*fi).P(0) )); + angle1 = math::Abs(Angle( (*fi).P(0)-(*fi).P(1),(*fi).P(2)-(*fi).P(1) )); + angle2 = M_PI-(angle0+angle1); + + // Skip degenerate triangles. + if(angle0==0 || angle1==0 || angle1==0) continue; + + e01v = ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ; + e12v = ( (*fi).V(2)->cP() - (*fi).V(1)->cP() ) ; + e20v = ( (*fi).V(0)->cP() - (*fi).V(2)->cP() ) ; + + TDContr[(*fi).V(0)] += ( e20v * (1.0/tan(angle1)) - e01v * (1.0/tan(angle2)) ) / 4.0; + TDContr[(*fi).V(1)] += ( e01v * (1.0/tan(angle2)) - e12v * (1.0/tan(angle0)) ) / 4.0; + TDContr[(*fi).V(2)] += ( e12v * (1.0/tan(angle0)) - e20v * (1.0/tan(angle1)) ) / 4.0; + + (*fi).V(0)->Kg() -= angle0; + (*fi).V(1)->Kg() -= angle1; + (*fi).V(2)->Kg() -= angle2; + + + for(int i=0;i<3;i++) + { + if(vcg::face::IsBorder((*fi), i)) + { + CoordType e1,e2; + vcg::face::Pos hp(&*fi, i, (*fi).V(i)); + vcg::face::Pos hp1=hp; + + hp1.FlipV(); + e1=hp1.v->cP() - hp.v->cP(); + hp1.FlipV(); + hp1.NextB(); + e2=hp1.v->cP() - hp.v->cP(); + (*fi).V(i)->Kg() -= math::Abs(Angle(e1,e2)); + } + } + } + + for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD() /*&& !(*vi).IsB()*/) + { + if((TDAreaPtr)[*vi].A<=std::numeric_limits::epsilon()) + { + (*vi).Kh() = 0; + (*vi).Kg() = 0; + } + else + { + (*vi).Kh() = (((TDContr)[*vi].dot((*vi).cN())>0)?1.0:-1.0)*((TDContr)[*vi] / (TDAreaPtr) [*vi].A).Norm(); + (*vi).Kg() /= (TDAreaPtr)[*vi].A; + } + } +} /// \brief Update the mean and the gaussian curvature of a vertex. @@ -531,12 +529,8 @@ public: Based on the paper "Optimizing 3d triangulations using discrete curvature analysis" */ - static float VertexCurvature(VertexPointer v, bool norm = true) + static float ComputeSingleVertexCurvature(VertexPointer v, bool norm = true) { - // VFAdjacency required! - assert(FaceType::HasVFAdjacency()); - assert(VertexType::HasVFAdjacency()); - VFIteratorType vfi(v); float A = 0; @@ -595,10 +589,14 @@ public: return A; } - static void VertexCurvature(MeshType & m){ + static void PerVertex(MeshType & m) + { + // VFAdjacency required! + assert(FaceType::HasVFAdjacency()); + assert(VertexType::HasVFAdjacency()); for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) - VertexCurvature(&*vi,false); + ComputeSingleVertexCurvature(&*vi,false); } @@ -613,12 +611,11 @@ public: } */ - static void PrincipalDirectionsNormalCycles(MeshType & m){ + static void PrincipalDirectionsNormalCycle(MeshType & m){ assert(VertexType::HasVFAdjacency()); assert(FaceType::HasFFAdjacency()); assert(FaceType::HasFaceNormal()); - typename MeshType::VertexIterator vi; for(vi = m.vert.begin(); vi != m.vert.end(); ++vi) @@ -638,7 +635,7 @@ public: Point3 n1 = p.F()->cN();n1.Normalize(); Point3 n2 = p.FFlip()->cN();n2.Normalize(); ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge); - n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0)); + n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0)); ScalarType beta = math::Asin(n1n2); m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0]; m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0]; @@ -657,28 +654,36 @@ public: m33[2][0] = m33[0][2]; m33[2][1] = m33[1][2]; + Eigen::Matrix3d it; + m33.ToEigenMatrix(it); + Eigen::SelfAdjointEigenSolver eig(it); + Eigen::Vector3d c_val = eig.eigenvalues(); + Eigen::Matrix3d c_vec = eig.eigenvectors(); + Point3 lambda; Matrix33 vect; - int n_rot; - Jacobi(m33,lambda, vect,n_rot); + vect.FromEigenMatrix(c_vec); + lambda.FromEigenVector(c_val); - vect.transposeInPlace(); - ScalarType normal = std::numeric_limits::min(); - int normI = 0; + ScalarType bestNormal = 0; + int bestNormalIndex = -1; for(int i = 0; i < 3; ++i) - if( fabs((*vi).N().Normalize().dot(vect.GetRow(i))) > normal ) + { + float agreeWithNormal = fabs((*vi).N().Normalize().dot(vect.GetColumn(i))); + if( agreeWithNormal > bestNormal ) { - normal= fabs((*vi).N().Normalize().dot(vect.GetRow(i))); - normI = i; + bestNormal= agreeWithNormal; + bestNormalIndex = i; } - int maxI = (normI+2)%3; - int minI = (normI+1)%3; + } + int maxI = (bestNormalIndex+2)%3; + int minI = (bestNormalIndex+1)%3; if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI); (*vi).PD1() = *(Point3*)(& vect[maxI][0]); (*vi).PD2() = *(Point3*)(& vect[minI][0]); - (*vi).K1() = lambda[maxI]; - (*vi).K2() = lambda[minI]; + (*vi).K1() = lambda[2]; + (*vi).K2() = lambda[1]; } } }; diff --git a/vcg/complex/algorithms/update/curvature_fitting.h b/vcg/complex/algorithms/update/curvature_fitting.h index 4e8caf58..f3985f5a 100644 --- a/vcg/complex/algorithms/update/curvature_fitting.h +++ b/vcg/complex/algorithms/update/curvature_fitting.h @@ -28,27 +28,21 @@ #define VCGLIB_UPDATE_CURVATURE_FITTING #include -#include -#include #include #include #include #include #include #include -#include #include #include -#include +#include #include #include #include #include #include -// GG include -#include -#include namespace vcg { @@ -116,17 +110,17 @@ class Quadric return 2.0*c()*v + b()*u + e(); } - double duv(double u, double v) + double duv(double /*u*/, double /*v*/) { return b(); } - double duu(double u, double v) + double duu(double /*u*/, double /*v*/) { return 2.0*a(); } - double dvv(double u, double v) + double dvv(double /*u*/, double /*v*/) { return 2.0*c(); } @@ -231,8 +225,9 @@ class Quadric static void computeCurvature(MeshType & m) { + Allocator::CompactVertexVector(m); - assert(tri::HasPerVertexVFAdjacency(m) && tri::HasPerFaceVFAdjacency(m) ); + if(!HasFVAdjacency(m)) throw vcg::MissingComponentException("FVAdjacency"); vcg::tri::UpdateTopology::VertexFace(m); @@ -357,17 +352,17 @@ class Quadric return 2.0*c()*v + b()*u + e(); } - double duv(double u, double v) + double duv(double /*u*/, double /*v*/) { return b(); } - double duu(double u, double v) + double duu(double /*u*/, double /*v*/) { return 2.0*a(); } - double dvv(double u, double v) + double dvv(double /*u*/, double /*v*/) { return 2.0*c(); }