From 8a265c9d321f9027af7776af5fbd524c091d3d9a Mon Sep 17 00:00:00 2001 From: ganovelli Date: Mon, 17 Mar 2008 11:29:59 +0000 Subject: [PATCH] taubin and desbrun estimates added (-> see vcg/simplex/vertexplus/component.h [component_ocf.h|component_occ.h ] --- vcg/complex/trimesh/update/curvature.h | 432 ++++++++++++++++++------- 1 file changed, 315 insertions(+), 117 deletions(-) diff --git a/vcg/complex/trimesh/update/curvature.h b/vcg/complex/trimesh/update/curvature.h index d201fb1e..017539bd 100644 --- a/vcg/complex/trimesh/update/curvature.h +++ b/vcg/complex/trimesh/update/curvature.h @@ -23,6 +23,11 @@ /**************************************************************************** History $Log: not supported by cvs2svn $ +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 @@ -36,8 +41,11 @@ the vertex #define VCGLIB_UPDATE_CURVATURE_ #include +#include #include #include +#include +#include namespace vcg { namespace tri { @@ -47,140 +55,330 @@ namespace tri { /// Management, updating and computation of per-vertex and per-face normals. /// This class is used to compute or update the normals that can be stored in the vertex or face component of a mesh. -template +template class UpdateCurvature { public: -typedef ComputeMeshType MeshType; -typedef typename MeshType::VertexType VertexType; -typedef typename VertexType::NormalType NormalType; -typedef typename VertexType::ScalarType ScalarType; -typedef typename MeshType::VertexPointer VertexPointer; -typedef typename MeshType::VertexIterator VertexIterator; -typedef typename MeshType::FaceType FaceType; -typedef typename MeshType::FacePointer FacePointer; -typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::VertContainer VertContainer; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::CoordType CoordType; + typedef typename CoordType::ScalarType ScalarType; + + +private: + typedef struct AdjVertex { + VertexType * vert; + float doubleArea; + bool isBorder; + }; +public: + /* + Compute principal direction and magniuto of curvature as describe in the paper: + @InProceedings{bb33922, + author = "G. Taubin", + title = "Estimating the Tensor of Curvature of a Surface from a + Polyhedral Approximation", + booktitle = "International Conference on Computer Vision", + year = "1995", + pages = "902--907", + URL = "http://dx.doi.org/10.1109/ICCV.1995.466840", + bibsource = "http://www.visionbib.com/bibliography/describe440.html#TT32253", + } + */ + static void PrincipalDirections(MeshType &m) { + + assert(m.HasVFTopology()); + + vcg::tri::UpdateNormals::PerVertexNormalized(m); + + VertexIterator vi; + for (vi =m.vert.begin(); vi !=m.vert.end(); ++vi) { + if ( ! (*vi).IsD() && (*vi).VFp() != NULL) { + + VertexType * central_vertex = &(*vi); + + std::vector weights; + std::vector vertices; + + vcg::face::JumpingPos pos((*vi).VFp(), central_vertex); + + VertexType* firstV = pos.VFlip(); + VertexType* tempV; + float totalDoubleAreaSize = 0.0f; + + if (((firstV->P()-central_vertex->P())^(pos.VFlip()->P()-central_vertex->P()))*central_vertex->N()<=0.0f) + { + pos.Set(central_vertex->VFp(), central_vertex); + pos.FlipE(); + firstV = pos.VFlip(); + } + else pos.Set(central_vertex->VFp(), central_vertex); + + do + { + pos.NextE(); + tempV = pos.VFlip(); + + AdjVertex v; + + v.isBorder = pos.IsBorder(); + v.vert = tempV; + v.doubleArea = ((pos.F()->V(1)->P() - pos.F()->V(0)->P()) ^ (pos.F()->V(2)->P()- pos.F()->V(0)->P())).Norm();; + totalDoubleAreaSize += v.doubleArea; + + vertices.push_back(v); + } + while(tempV != firstV); + + for (int 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); + } + + Matrix33f Tp; + for (int i = 0; i < 3; ++i) + Tp[i][i] = 1.0f - powf(central_vertex->N()[i],2); + Tp[0][1] = Tp[1][0] = -1.0f * (central_vertex->N()[0] * central_vertex->N()[1]); + Tp[1][2] = Tp[2][1] = -1.0f * (central_vertex->N()[1] * central_vertex->N()[2]); + Tp[0][2] = Tp[2][0] = -1.0f * (central_vertex->N()[0] * central_vertex->N()[2]); + + Matrix33f tempMatrix; + Matrix33f M; + M.SetZero(); + for (int i = 0; i < vertices.size(); ++i) { + Point3f edge = (central_vertex->P() - vertices[i].vert->P()); + float curvature = (2.0f * (central_vertex->N() * edge) ) / edge.SquaredNorm(); + Point3f T = (Tp*edge).Normalize(); + tempMatrix.ExternalProduct(T,T); + M += tempMatrix * weights[i] * curvature ; + } + + Point3f W; + Point3f e1(1.0f,0.0f,0.0f); + if ((e1 - central_vertex->N()).SquaredNorm() > (e1 + central_vertex->N()).SquaredNorm()) + W = e1 - central_vertex->N(); + else + W = e1 + central_vertex->N(); + W.Normalize(); + + Matrix33f Q; + Q.SetIdentity(); + tempMatrix.ExternalProduct(W,W); + Q -= tempMatrix * 2.0f; + + Matrix33f Qt(Q); + Qt.Transpose(); + + Matrix33f QtMQ = (Qt * M * Q); + + Point3f bl = Q.GetColumn(0); + Point3f T1 = Q.GetColumn(1); + Point3f T2 = Q.GetColumn(2); + + float s,c; + // Gabriel Taubin hint and Valentino Fiorin impementation + float qt21 = QtMQ[2][1]; + float qt12 = QtMQ[1][2]; + 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 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[1], 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; + + 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 TDAreaPtr(m.vert); TDAreaPtr.Start(); - //reset the values to 0 - for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD()) - (*vi).Q() = 0.0; + //Calcola AreaMix in H (vale anche per K) + for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi) if(!(*vi).IsD()) + { + (TDAreaPtr)[*vi].A = 0; + (*vi).H() = 0; + (*vi).K() = (float)(2.0 * M_PI); + } - //--- compute Areamix - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - - // angles - angle0 = math::Abs(Angle( (*fi).V(1)->P()-(*fi).V(0)->P(),(*fi).V(2)->P()-(*fi).V(0)->P() )); - angle1 = math::Abs(Angle( (*fi).V(0)->P()-(*fi).V(1)->P(),(*fi).V(2)->P()-(*fi).V(1)->P() )); - 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)->P() , (*fi).V(0)->P() ); - float e12 = SquaredDistance( (*fi).V(2)->P() , (*fi).V(1)->P() ); - float e20 = SquaredDistance( (*fi).V(0)->P() , (*fi).V(2)->P() ); - - // voronoi area v[0] - area0 = ( e01*(1/tan(angle2)) + e20*(1/tan(angle1)) ) /8; - // voronoi area v[1] - area1 = ( e01*(1/tan(angle2)) + e12*(1/tan(angle0)) ) /8; - // voronoi area v[2] - area2 = ( e20*(1/tan(angle1)) + e20*(1/tan(angle0)) ) /8; - - (*fi).V(0)->Q() += area0; - (*fi).V(1)->Q() += area1; - (*fi).V(2)->Q() += area2; - } - else // triangolo ottuso - { - (*fi).V(0)->Q() += (*fi).DoubleArea() / 6; - (*fi).V(1)->Q() += (*fi).DoubleArea() / 6; - (*fi).V(2)->Q() += (*fi).DoubleArea() / 6; - } - } + 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)->P() , (*fi).V(0)->P() ); + float e12 = SquaredDistance( (*fi).V(2)->P() , (*fi).V(1)->P() ); + float e20 = SquaredDistance( (*fi).V(0)->P() , (*fi).V(2)->P() ); + + 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; - i = 0; - for(vi=m.vert.begin();vi!=m.vert.end();++vi,++i) if(!(*vi).IsD()) - { - area[i] = (*vi).Q(); - (*vi).Q() = (float)(2.0 * M_PI); - } - - if(false) - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - for(int i=0;i<3;i++) - { - if(vcg::face::IsBorder((*fi), i)) - { - typename MeshType::CoordType e1,e2; - vcg::face::Pos hp(&*fi,i,(*fi).V(i)); - //MeshType::hedgepos_type hp(&*fi,i,(*fi).V(i)); - vcg::face::Pos hp1=hp; - //MeshType::hedgepos_type hp1=hp; + } + else // triangolo ottuso + { + (TDAreaPtr)[(*fi).V(0)].A += vcg::DoubleArea((*fi)) / 6.0; + (TDAreaPtr)[(*fi).V(1)].A += vcg::DoubleArea((*fi)) / 6.0; + (TDAreaPtr)[(*fi).V(2)].A += vcg::DoubleArea((*fi)) / 6.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); + + e01 = ( (*fi).V(1)->P() - (*fi).V(0)->P() ) * (*fi).V(0)->N(); + e12 = ( (*fi).V(2)->P() - (*fi).V(1)->P() ) * (*fi).V(1)->N(); + e20 = ( (*fi).V(0)->P() - (*fi).V(2)->P() ) * (*fi).V(2)->N(); + + area0 = ( e20 * (1.0/tan(angle1)) + e01 * (1.0/tan(angle2)) ) / 4.0; + area1 = ( e01 * (1.0/tan(angle2)) + e12 * (1.0/tan(angle0)) ) / 4.0; + area2 = ( e12 * (1.0/tan(angle0)) + e20 * (1.0/tan(angle1)) ) / 4.0; + + (*fi).V(0)->H() += area0; + (*fi).V(1)->H() += area1; + (*fi).V(2)->H() += area2; - hp1.FlipV(); - - e1= hp1.v->P()-hp.v->P(); - hp1.FlipV(); - hp1.NextB(); - e2= hp1.v->P()-hp.v->P(); - (*fi).V(i)->Q() -=math::Abs(Angle(e1,e2)); + (*fi).V(0)->K() -= angle0; + (*fi).V(1)->K() -= angle1; + (*fi).V(2)->K() -= 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->P() - hp.v->P(); + hp1.FlipV(); + hp1.NextB(); + e2=hp1.v->P() - hp.v->P(); + (*fi).V(i)->K() -= 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).H() = 0; + (*vi).K() = 0; + } + else + { + (*vi).H() /= (TDAreaPtr)[*vi].A; + (*vi).K() /= (TDAreaPtr)[*vi].A; + } } - } - } - - - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - float angle0 = math::Abs(Angle( - (*fi).V(1)->P()-(*fi).V(0)->P(),(*fi).V(2)->P()-(*fi).V(0)->P() )); - float angle1 = math::Abs(Angle( - (*fi).V(0)->P()-(*fi).V(1)->P(),(*fi).V(2)->P()-(*fi).V(1)->P() )); - float angle2 = M_PI-(angle0+angle1); - - (*fi).V(0)->Q() -= angle0; - (*fi).V(1)->Q() -= angle1; - (*fi).V(2)->Q() -= angle2; - } - i=0; - for(vi=m.vert.begin(); vi!=m.vert.end(); ++vi,++i) if(!(*vi).IsD()) - { - (*vi).Q() /= area[i]; - (*vi).Q()=math::Clamp((*vi).Q(),-0.050f,0.050f); - - /* if ( (*vi).Q() < 0 ) - (*vi).Q() = log( -(*vi).Q() ); - else if( (*vi).Q() > 0 ) - (*vi).Q() = log( (*vi).Q() );*/ - } - - //--- DeInit - - delete[] area; + TDAreaPtr.Stop(); -} + } + }; + + } } #endif