From 82ddb476a488fcee0e2aebf5aaf7ba50df592f7b Mon Sep 17 00:00:00 2001 From: Paolo Cignoni Date: Tue, 21 Feb 2017 17:46:46 +0100 Subject: [PATCH] Heavy refactoring. Closing #12 Many changes, improved general robustness and added more options to customise the behaviour. Added control on quality quadric, Hard normal flipping check, SVDPlacement that find better optimal position and many other small optimizations. --- .../tri_edge_collapse_quadric.h | 600 ++++++++++-------- 1 file changed, 344 insertions(+), 256 deletions(-) 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