diff --git a/vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h b/vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h index eb6b86b1..35425c98 100644 --- a/vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h +++ b/vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h @@ -38,8 +38,6 @@ namespace vcg{ namespace tri{ - - /** This class describe Quadric based collapse operation. @@ -87,48 +85,32 @@ namespace tri{ class TriEdgeCollapseQuadricParameter : public BaseParameterClass { public: - double BoundaryWeight; - double CosineThr; - bool FastPreserveBoundary; - bool NormalCheck; - double NormalThrRad; - bool OptimalPlacement; - bool PreserveTopology; - bool PreserveBoundary; - double QuadricEpsilon; - bool QualityCheck; - bool QualityQuadric; // During the initialization manage all the edges as border edges adding a set of additional quadrics that are useful mostly for keeping face aspect ratio good. - double QualityThr; // Collapsed that generate faces with quality LOWER than this value are penalized. So - bool QualityWeight; - double QualityWeightFactor; - double ScaleFactor; - bool ScaleIndependent; - bool UseArea; - bool UseVertexWeight; + double BoundaryQuadricWeight = 0.5; + bool FastPreserveBoundary = false; + bool AreaCheck = false; + bool HardQualityCheck = false; + double HardQualityThr = 0.1; + bool HardNormalCheck = false; + bool NormalCheck = false; + double NormalThrRad = M_PI/2.0; + double CosineThr = 0 ; // ~ cos(pi/2) + bool OptimalPlacement =true; + bool SVDPlacement = false; + bool PreserveTopology =false; + bool PreserveBoundary = false; + double QuadricEpsilon = 1e-15; + bool QualityCheck =true; + double QualityThr =.3; // Collapsed that generate faces with quality LOWER than this value are penalized. So higher the value -> better the quality of the accepted triangles + bool QualityQuadric =false; // During the initialization manage all the edges as border edges adding a set of additional quadrics that are useful mostly for keeping face aspect ratio good. + double QualityQuadricWeight = 0.001f; // During the initialization manage all the edges as border edges adding a set of additional quadrics that are useful mostly for keeping face aspect ratio good. + bool QualityWeight=false; + double QualityWeightFactor=100.0; + double ScaleFactor=1.0; + bool ScaleIndependent=true; + bool UseArea =true; + bool UseVertexWeight=false; - void SetDefaultParams() - { - BoundaryWeight=.5; - CosineThr=cos(M_PI/2); - FastPreserveBoundary=false; - NormalCheck=false; - NormalThrRad=M_PI/2; - OptimalPlacement=true; - PreserveBoundary = false; - PreserveTopology = false; - QuadricEpsilon =1e-15; - QualityCheck=true; - QualityQuadric=false; - QualityThr=.3; // higher the value -> better the quality of the accepted triangles - QualityWeight=false; - QualityWeightFactor=100.0; - ScaleFactor=1.0; - ScaleIndependent=true; - UseArea=true; - UseVertexWeight=false; - } - - TriEdgeCollapseQuadricParameter() {this->SetDefaultParams();} + TriEdgeCollapseQuadricParameter() {} }; @@ -150,8 +132,9 @@ public: typedef TriEdgeCollapseQuadricParameter QParameter; typedef HelperType QH; + CoordType optimalPos; // Local storage of the once computed optimal position of the collapse. - // puntatori ai vertici che sono stati messi non-w per preservare il boundary + // Pointer to the vector that store the Write flags. Used to preserve them if you ask to preserve for the boundaries. static std::vector & WV(){ static std::vector _WV; return _WV; } @@ -166,141 +149,146 @@ public: } - inline bool IsFeasible(BaseParameterClass *_pp){ - QParameter *pp=(QParameter *)_pp; - if(!pp->PreserveTopology) return true; - - bool res = ( EdgeCollapser::LinkConditions(this->pos) ); - if(!res) ++( TEC::FailStat::LinkConditionEdge() ); - return res; - } - - CoordType ComputePosition(BaseParameterClass *_pp) - { - QParameter *pp=(QParameter *)_pp; - CoordType newPos = (this->pos.V(0)->P()+this->pos.V(1)->P())/2.0; - if(pp->OptimalPlacement) - { - if((QH::Qd(this->pos.V(0)).Apply(newPos) + QH::Qd(this->pos.V(1)).Apply(newPos)) > 200.0*pp->QuadricEpsilon) - newPos = ComputeMinimal(); - } - else newPos=this->pos.V(1)->P(); - return newPos; - } - - void Execute(TriMeshType &m, BaseParameterClass *_pp) - { - QH::Qd(this->pos.V(1))+=QH::Qd(this->pos.V(0)); - EdgeCollapser::Do(m, this->pos, ComputePosition(_pp)); // v0 is deleted and v1 take the new position + inline bool IsFeasible(BaseParameterClass *_pp){ + QParameter *pp=(QParameter *)_pp; + if(!pp->PreserveTopology) return true; + + bool res = ( EdgeCollapser::LinkConditions(this->pos) ); + if(!res) ++( TEC::FailStat::LinkConditionEdge() ); + return res; } - - - - // Final Clean up after the end of the simplification process - static void Finalize(TriMeshType &m, HeapType& /*h_ret*/, BaseParameterClass *_pp) - { - QParameter *pp=(QParameter *)_pp; - - // If we had the boundary preservation we should clean up the writable flags - if(pp->FastPreserveBoundary) - { - typename TriMeshType::VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD()) (*vi).SetW(); - } - if(pp->PreserveBoundary) - { - typename std::vector::iterator wvi; - for(wvi=WV().begin();wvi!=WV().end();++wvi) - if(!(*wvi)->IsD()) (*wvi)->SetW(); - } - } - - static void Init(TriMeshType &m, HeapType &h_ret, BaseParameterClass *_pp) + + void ComputePosition(BaseParameterClass *_pp) { QParameter *pp=(QParameter *)_pp; - - typename TriMeshType::VertexIterator vi; - typename TriMeshType::FaceIterator pf; - - pp->CosineThr=cos(pp->NormalThrRad); - - vcg::tri::UpdateTopology::VertexFace(m); - vcg::tri::UpdateFlags::FaceBorderFromVF(m); - - if(pp->FastPreserveBoundary) + CoordType newPos = (this->pos.V(0)->P()+this->pos.V(1)->P())/2.0; + if(pp->OptimalPlacement==false) + newPos=this->pos.V(1)->P(); + else { - for(pf=m.face.begin();pf!=m.face.end();++pf) - if( !(*pf).IsD() && (*pf).IsW() ) - for(int j=0;j<3;++j) - if((*pf).IsB(j)) - { - (*pf).V(j)->ClearW(); - (*pf).V1(j)->ClearW(); - } - } - - if(pp->PreserveBoundary) + if((QH::Qd(this->pos.V(0)).Apply(newPos) + QH::Qd(this->pos.V(1)).Apply(newPos)) > 2.0*pp->QuadricEpsilon) + { + QuadricType q=QH::Qd(this->pos.V(0)); + q+=QH::Qd(this->pos.V(1)); + + Point3 x; + if(pp->SVDPlacement) + q.MinimumClosestToPoint(x,Point3d::Construct(newPos)); + else + q.Minimum(x); + newPos = CoordType::Construct(x); + } + } + this->optimalPos = newPos; + } + + void Execute(TriMeshType &m, BaseParameterClass */*_pp*/) + { + CoordType newPos = this->optimalPos; + QH::Qd(this->pos.V(1))+=QH::Qd(this->pos.V(0)); // v0 is deleted and v1 take the new position + EdgeCollapser::Do(m, this->pos, newPos); + } + + // Final Clean up after the end of the simplification process + static void Finalize(TriMeshType &m, HeapType& /*h_ret*/, BaseParameterClass *_pp) + { + QParameter *pp=(QParameter *)_pp; + + // If we had the boundary preservation we should clean up the writable flags + if(pp->FastPreserveBoundary) { - WV().clear(); - for(pf=m.face.begin();pf!=m.face.end();++pf) - if( !(*pf).IsD() && (*pf).IsW() ) - for(int j=0;j<3;++j) - if((*pf).IsB(j)) - { - if((*pf).V(j)->IsW()) {(*pf).V(j)->ClearW(); WV().push_back((*pf).V(j));} - if((*pf).V1(j)->IsW()) {(*pf).V1(j)->ClearW();WV().push_back((*pf).V1(j));} - } - } - - InitQuadric(m,pp); - - // Initialize the heap with all the possible collapses - if(IsSymmetric(pp)) - { // if the collapse is symmetric (e.g. u->v == v->u) + typename TriMeshType::VertexIterator vi; for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && (*vi).IsRW()) + if(!(*vi).IsD()) (*vi).SetW(); + } + if(pp->PreserveBoundary) + { + typename std::vector::iterator wvi; + for(wvi=WV().begin();wvi!=WV().end();++wvi) + if(!(*wvi)->IsD()) (*wvi)->SetW(); + } + } + + static void Init(TriMeshType &m, HeapType &h_ret, BaseParameterClass *_pp) + { + QParameter *pp=(QParameter *)_pp; + pp->CosineThr=cos(pp->NormalThrRad); + h_ret.clear(); + vcg::tri::UpdateTopology::VertexFace(m); + vcg::tri::UpdateFlags::FaceBorderFromVF(m); + + if(pp->FastPreserveBoundary) + { + for(auto pf=m.face.begin();pf!=m.face.end();++pf) + if( !(*pf).IsD() && (*pf).IsW() ) + for(int j=0;j<3;++j) + if((*pf).IsB(j)) { - vcg::face::VFIterator x; - for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x){ - x.V1()->ClearV(); - x.V2()->ClearV(); - } - for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++x ) - { - assert(x.F()->V(x.I())==&(*vi)); - if((x.V0()IsRW() && !x.V1()->IsV()){ - x.V1()->SetV(); - h_ret.push_back(HeapElem(new MYTYPE(VertexPair(x.V0(),x.V1()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp ))); - } - if((x.V0()IsRW()&& !x.V2()->IsV()){ - x.V2()->SetV(); - h_ret.push_back(HeapElem(new MYTYPE(VertexPair(x.V0(),x.V2()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp ))); - } - } + (*pf).V(j)->ClearW(); + (*pf).V1(j)->ClearW(); } } - else + + if(pp->PreserveBoundary) + { + WV().clear(); + for(auto pf=m.face.begin();pf!=m.face.end();++pf) + if( !(*pf).IsD() && (*pf).IsW() ) + for(int j=0;j<3;++j) + if((*pf).IsB(j)) + { + if((*pf).V(j)->IsW()) {(*pf).V(j)->ClearW(); WV().push_back((*pf).V(j));} + if((*pf).V1(j)->IsW()) {(*pf).V1(j)->ClearW();WV().push_back((*pf).V1(j));} + } + } + + InitQuadric(m,pp); + + // Initialize the heap with all the possible collapses + if(IsSymmetric(pp)) + { // if the collapse is symmetric (e.g. u->v == v->u) + for(auto vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && (*vi).IsRW()) + { + vcg::face::VFIterator x; + for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x){ + x.V1()->ClearV(); + x.V2()->ClearV(); + } + for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++x ) + { + if((x.V0()IsRW() && !x.V1()->IsV()){ + x.V1()->SetV(); + h_ret.push_back(HeapElem(new MYTYPE(VertexPair(x.V0(),x.V1()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp ))); + } + if((x.V0()IsRW()&& !x.V2()->IsV()){ + x.V2()->SetV(); + h_ret.push_back(HeapElem(new MYTYPE(VertexPair(x.V0(),x.V2()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp ))); + } + } + } + } + else { // if the collapse is A-symmetric (e.g. u->v != v->u) - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && (*vi).IsRW()) - { - vcg::face::VFIterator x; - UnMarkAll(m); - for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x) - { - assert(x.F()->V(x.I())==&(*vi)); - if(x.V()->IsRW() && x.V1()->IsRW() && !IsMarked(m,x.F()->V1(x.I()))){ - h_ret.push_back( HeapElem( new MYTYPE( VertexPair (x.V(),x.V1()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp))); - } - if(x.V()->IsRW() && x.V2()->IsRW() && !IsMarked(m,x.F()->V2(x.I()))){ - h_ret.push_back( HeapElem( new MYTYPE( VertexPair (x.V(),x.V2()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp))); - } - } - } - } -} - static float HeapSimplexRatio(BaseParameterClass *_pp) {return IsSymmetric(_pp)?5.0f:9.0f;} + for(auto vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && (*vi).IsRW()) + { + vcg::face::VFIterator x; + UnMarkAll(m); + for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x) + { + if(x.V()->IsRW() && x.V1()->IsRW() && !IsMarked(m,x.F()->V1(x.I()))){ + h_ret.push_back( HeapElem( new MYTYPE( VertexPair (x.V(),x.V1()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp))); + } + if(x.V()->IsRW() && x.V2()->IsRW() && !IsMarked(m,x.F()->V2(x.I()))){ + h_ret.push_back( HeapElem( new MYTYPE( VertexPair (x.V(),x.V2()),TriEdgeCollapse< TriMeshType,VertexPair,MYTYPE>::GlobalMark(),_pp))); + } + } + } + } + } +// static float HeapSimplexRatio(BaseParameterClass *_pp) {return IsSymmetric(_pp)?5.0f:9.0f;} + static float HeapSimplexRatio(BaseParameterClass *_pp) {return IsSymmetric(_pp)?4.0f:8.0f;} static bool IsSymmetric(BaseParameterClass *_pp) {return ((QParameter *)_pp)->OptimalPlacement;} static bool IsVertexStable(BaseParameterClass *_pp) {return !((QParameter *)_pp)->OptimalPlacement;} @@ -315,60 +303,83 @@ public: ScalarType ComputePriority(BaseParameterClass *_pp) { QParameter *pp=(QParameter *)_pp; - std::vector onVec; // vector with incident faces original normals + VertexType * v[2]; v[0] = this->pos.V(0); v[1] = this->pos.V(1); - if(pp->NormalCheck){ // Compute maximal normal variation - // store the old normals for non-collapsed face in v0 + std::vector origNormalVec; // vector with incident faces original normals + if(pp->NormalCheck){ // Collect Original Normals for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 - if( x.V1()!=v[1] && x.V2()!=v[1] ) // skip faces with v1 - onVec.push_back(TriangleNormal(*x.F()).Normalize()); - // store the old normals for non-collapsed face in v1 + if( x.V1()!=v[1] && x.V2()!=v[1] ) // skip faces with v1 + origNormalVec.push_back(NormalizedTriangleNormal(*x.F())); for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 - if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 - onVec.push_back(TriangleNormal(*x.F()).Normalize()); + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + origNormalVec.push_back(NormalizedTriangleNormal(*x.F())); } + ScalarType origArea=0; + if(pp->AreaCheck){ // Collect Original Area + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + origArea += DoubleArea(*x.F()); + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + origArea += DoubleArea(*x.F()); + } + + ScalarType origQual= std::numeric_limits::max(); + if(pp->HardQualityCheck){ // Collect original quality + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + origQual=std::min(origQual, QualityFace(*x.F())); + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + origQual=std::min(origQual, QualityFace(*x.F())); + } + + //// Move the two vertexes into new position (storing the old ones) CoordType OldPos0=v[0]->P(); CoordType OldPos1=v[1]->P(); - CoordType newPos = ComputePosition(_pp); - v[0]->P() = v[1]->P() = newPos; + ComputePosition(_pp); + // Now Simulate the collapse + v[0]->P() = v[1]->P() = this->optimalPos; + + //// Rescan faces and compute the worst quality and normals that happens after collapse + ScalarType MinCos = std::numeric_limits::max(); // Cosine of the angle variation: -1 ~ very bad to 1~perfect + if(pp->NormalCheck){ + int i=0; + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + if( x.V1()!=v[1] && x.V2()!=v[1] ) // skipping faces with v1 + { + CoordType nn=NormalizedTriangleNormal(*x.F()); + MinCos=std::min(MinCos,nn.dot(origNormalVec[i++])); + } + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + { + CoordType nn=NormalizedTriangleNormal(*x.F()); + MinCos=std::min(MinCos,nn.dot(origNormalVec[i++])); + } + } - //// Rescan faces and compute quality and difference between normals - int i=0; - double MinCos = std::numeric_limits::max(); // minimo coseno di variazione di una normale della faccia - // (e.g. max angle) Mincos varia da 1 (normali coincidenti) a - // -1 (normali opposte); - double MinQual = std::numeric_limits::max(); - for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 - if( x.V1()!=v[1] && x.V2()!=v[1] ) // skiping faces with v1 - { - if(pp->NormalCheck){ - CoordType nn=NormalizedTriangleNormal(*x.F()); - double ndiff=nn.dot(onVec[i++]); - MinCos=std::min(MinCos,ndiff); - } - if(pp->QualityCheck){ - double qt= QualityFace(*x.F()); - MinQual=std::min(MinQual,qt); - } - } - for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 - if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 - { - if(pp->NormalCheck){ - CoordType nn=NormalizedTriangleNormal(*x.F()); - double ndiff=nn.dot(onVec[i++]); - MinCos=std::min(MinCos,ndiff); - } - if(pp->QualityCheck){ - double qt= QualityFace(*x.F()); - MinQual=std::min(MinQual,qt); - } - } + ScalarType newQual = std::numeric_limits::max(); // + if(pp->QualityCheck){ + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + if( x.V1()!=v[1] && x.V2()!=v[1] ) + newQual=std::min(newQual,QualityFace(*x.F())); + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + newQual=std::min(newQual,QualityFace(*x.F())); + } + + ScalarType newArea=0; + if(pp->AreaCheck){ // Collect Area + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + newArea += DoubleArea(*x.F()); + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + newArea += DoubleArea(*x.F()); + } QuadricType qq=QH::Qd(v[0]); qq+=QH::Qd(v[1]); @@ -377,28 +388,40 @@ public: assert(!math::IsNAN(QuadErr)); // All collapses involving triangles with quality larger than have no penalty; - if(MinQual>pp->QualityThr) MinQual=pp->QualityThr; + if(newQual>pp->QualityThr) newQual=pp->QualityThr; if(pp->NormalCheck){ // All collapses where the normal vary less than (e.g. more than CosineThr) - // have no penalty + // have no particular penalty if(MinCos>pp->CosineThr) MinCos=pp->CosineThr; - MinCos=fabs((MinCos+1)/2.0); // Now it is in the range 0..1 with 0 very dangerous! + MinCos=fabs((MinCos+1.0)/2.0); // Now it is in the range 0..1 with 0 very bad! + assert(MinCos >=0 && MinCos<1.1 ); } + QuadErr= std::max(QuadErr,pp->QuadricEpsilon); if(QuadErr <= pp->QuadricEpsilon) { - QuadErr = - 1/Distance(OldPos0,OldPos1); + QuadErr *= Distance(OldPos0,OldPos1); } if( pp->UseVertexWeight ) QuadErr *= (QH::W(v[1])+QH::W(v[0]))/2; ScalarType error; if(!pp->QualityCheck && !pp->NormalCheck) error = (ScalarType)(QuadErr); - if( pp->QualityCheck && !pp->NormalCheck) error = (ScalarType)(QuadErr / MinQual); + if( pp->QualityCheck && !pp->NormalCheck) error = (ScalarType)(QuadErr / newQual); if(!pp->QualityCheck && pp->NormalCheck) error = (ScalarType)(QuadErr / MinCos); - if( pp->QualityCheck && pp->NormalCheck) error = (ScalarType)(QuadErr / (MinQual*MinCos)); + if( pp->QualityCheck && pp->NormalCheck) error = (ScalarType)(QuadErr / (newQual*MinCos)); + + if(pp->AreaCheck && ((fabs(origArea-newArea)/(origArea+newArea))>0.01) ) + error = std::numeric_limits::max(); + + if(pp->HardQualityCheck && + (newQual < pp->HardQualityThr && newQual < origQual*0.9) ) + error = std::numeric_limits::max(); + + if(pp->HardNormalCheck) + if(CheckForFlip()) error = std::numeric_limits::max(); // Restore old position of v0 and v1 v[0]->P()=OldPos0; @@ -408,17 +431,95 @@ public: return this->_priority; } -// -//static double MaxError() {return 1e100;} -// + + bool CheckForFlippedFaceOverVertex(VertexType *vp, ScalarType angleThrRad = math::ToRad(150.)) + { + std::map edgeNormMap; + ScalarType maxAngle=0; + + for(VFIterator x(vp); !x.End(); ++x ) // for all faces in v1 + { + if(QualityFace(*x.F()) <0.01 ) return true; + for(int i=0;i<2;++i) + { + VertexType *vv= i==0?x.V1():x.V2(); + assert(vv!=vp); + auto ni = edgeNormMap.find(vv); + if(ni==edgeNormMap.end()) edgeNormMap[vv] = NormalizedTriangleNormal(*x.F()); + else maxAngle = std::max(maxAngle,AngleN(NormalizedTriangleNormal(*x.F()),ni->second)); + } + } + + return (maxAngle > angleThrRad); + } + + // This function return true if, after an edge collapse, + // among the surviving faces, there are two adjacent faces forming a + // diedral angle larger than the given threshold + // It assumes that the two vertexes of the collapsing edge + // have been already moved to the new position but the topolgy has not yet been changed (e.g. there are two zero-area faces) + + bool CheckForFlip(ScalarType angleThrRad = math::ToRad(150.)) + { + std::map edgeNormMap; + VertexType * v[2]; + v[0] = this->pos.V(0); + v[1] = this->pos.V(1); + ScalarType maxAngle=0; + assert (v[0]->P()==v[1]->P()); + + for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 + if( x.V1()!=v[1] && x.V2()!=v[1] ) // skip faces with v1 + { + if(QualityFace(*x.F()) <0.01 ) return true; + for(int i=0;i<2;++i) + { + VertexType *vv= (i==0)?x.V1():x.V2(); + assert(vv!=v[0]); + auto ni = edgeNormMap.find(vv); + if(ni==edgeNormMap.end()) edgeNormMap[vv] = NormalizedTriangleNormal(*x.F()); + else maxAngle = std::max(maxAngle,AngleN(NormalizedTriangleNormal(*x.F()),ni->second)); + } + } + for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 + if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0 + { + if(QualityFace(*x.F()) <0.01 ) return true; + for(int i=0;i<2;++i) + { + VertexType *vv= i==0?x.V1():x.V2(); + assert(vv!=v[1]); + auto ni = edgeNormMap.find(vv); + if(ni==edgeNormMap.end()) edgeNormMap[vv] = NormalizedTriangleNormal(*x.F()); + else maxAngle = std::max(maxAngle,AngleN(NormalizedTriangleNormal(*x.F()),ni->second)); + } + } + return (maxAngle > angleThrRad); + } + + + + inline void AddCollapseToHeap(HeapType & h_ret, VertexType *v0, VertexType *v1, BaseParameterClass *_pp) { QParameter *pp=(QParameter *)_pp; + ScalarType maxAdmitErr = std::numeric_limits::max(); h_ret.push_back(HeapElem(new MYTYPE(VertexPair(v0,v1), this->GlobalMark(),_pp))); - std::push_heap(h_ret.begin(),h_ret.end()); + if(h_ret.back().pri > maxAdmitErr) { + delete h_ret.back().locModPtr; + h_ret.pop_back(); + } + else + std::push_heap(h_ret.begin(),h_ret.end()); + if(!IsSymmetric(pp)){ h_ret.push_back(HeapElem(new MYTYPE(VertexPair(v1,v0), this->GlobalMark(),_pp))); - std::push_heap(h_ret.begin(),h_ret.end()); + if(h_ret.back().pri > maxAdmitErr) { + delete h_ret.back().locModPtr; + h_ret.pop_back(); + } + else + std::push_heap(h_ret.begin(),h_ret.end()); } } @@ -434,6 +535,8 @@ public: for(VFIterator vfi(v[1]); !vfi.End(); ++vfi ) { vfi.V1()->ClearV(); vfi.V2()->ClearV(); + vfi.V1()->IMark() = this->GlobalMark(); + vfi.V2()->IMark() = this->GlobalMark(); } // Second Loop @@ -457,8 +560,7 @@ public: { QParameter *pp=(QParameter *)_pp; QH::Init(); - // m.ClearFlags(); - for(VertexIterator pv=m.vert.begin();pv!=m.vert.end();++pv) // Azzero le quadriche + for(VertexIterator pv=m.vert.begin();pv!=m.vert.end();++pv) if( ! (*pv).IsD() && (*pv).IsW()) QH::Qd(*pv).SetZero(); @@ -487,9 +589,9 @@ public: QuadricType bq; // Border quadric record the squared distance from the plane orthogonal to the face and passing // through the edge. - borderPlane.SetDirection(facePlane.Direction() ^ ( (*fi).V1(j)->cP() - (*fi).V(j)->cP() ).normalized()); - if( (*fi).IsB(j) ) borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->BoundaryWeight )); // amplify border planes - else borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->BoundaryWeight/100.0)); // and consider much less quadric for quality + borderPlane.SetDirection(facePlane.Direction() ^ (( (*fi).V1(j)->cP() - (*fi).V(j)->cP() ).normalized())); + if( (*fi).IsB(j) ) borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->BoundaryQuadricWeight )); // amplify border planes + else borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->QualityQuadricWeight )); // and consider much less quadric for quality borderPlane.SetOffset(borderPlane.Direction().dot((*fi).V(j)->cP())); bq.ByPlane(borderPlane); @@ -518,36 +620,24 @@ public: } } - - -// -// -// -// -// -// -//static void InitMesh(MESH_TYPE &m){ -// pp->CosineThr=cos(pp->NormalThr); -// InitQuadric(m); -// //m.Topology(); -// //OldInitQuadric(m,UseArea); -// } -// - CoordType ComputeMinimal() + CoordType ComputeMinimal() + { + } + +CoordType ComputeMinimalOld() { - typename TriMeshType::VertexType * v[2]; - v[0] = this->pos.V(0); - v[1] = this->pos.V(1); - QuadricType q=QH::Qd(v[0]); - q+=QH::Qd(v[1]); + VertexType* &v0 = this->pos.V(0); + VertexType* &v1 = this->pos.V(1); + QuadricType q=QH::Qd(v0); + q+=QH::Qd(v1); Point3 x; bool rt=q.Minimum(x); if(!rt) { // if the computation of the minimum fails we choose between the two edge points and the middle one. - Point3 x0=Point3d::Construct(v[0]->P()); - Point3 x1=Point3d::Construct(v[1]->P()); - x.Import((v[0]->P()+v[1]->P())/2); + Point3 x0=Point3d::Construct(v0->P()); + Point3 x1=Point3d::Construct(v1->P()); + x.Import((v1->P()+v1->P())/2); double qvx=q.Apply(x); double qv0=q.Apply(x0); double qv1=q.Apply(x1); @@ -557,10 +647,8 @@ public: return CoordType::Construct(x); } -// -// - }; - } // namespace tri - } // namespace vcg + +} // namespace tri +} // namespace vcg #endif