manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

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.
This commit is contained in:
Paolo Cignoni 2017-02-21 17:46:46 +01:00
parent 75bb6dff89
commit 82ddb476a4
1 changed files with 344 additions and 256 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

600
vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h
View File

@ -38,8 +38,6 @@ namespace vcg{
namespace tri{ namespace tri{
/** /**
This class describe Quadric based collapse operation. This class describe Quadric based collapse operation.
@ -87,48 +85,32 @@ namespace tri{
class TriEdgeCollapseQuadricParameter : public BaseParameterClass class TriEdgeCollapseQuadricParameter : public BaseParameterClass
{ {
public: public:
double BoundaryWeight; double BoundaryQuadricWeight = 0.5;
double CosineThr; bool FastPreserveBoundary = false;
bool FastPreserveBoundary; bool AreaCheck = false;
bool NormalCheck; bool HardQualityCheck = false;
double NormalThrRad; double HardQualityThr = 0.1;
bool OptimalPlacement; bool HardNormalCheck = false;
bool PreserveTopology; bool NormalCheck = false;
bool PreserveBoundary; double NormalThrRad = M_PI/2.0;
double QuadricEpsilon; double CosineThr = 0 ; // ~ cos(pi/2)
bool QualityCheck; bool OptimalPlacement =true;
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. bool SVDPlacement = false;
double QualityThr; // Collapsed that generate faces with quality LOWER than this value are penalized. So bool PreserveTopology =false;
bool QualityWeight; bool PreserveBoundary = false;
double QualityWeightFactor; double QuadricEpsilon = 1e-15;
double ScaleFactor; bool QualityCheck =true;
bool ScaleIndependent; 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 UseArea; 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.
bool UseVertexWeight; 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() TriEdgeCollapseQuadricParameter() {}
{
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();}
}; };
@ -150,8 +132,9 @@ public:
typedef TriEdgeCollapseQuadricParameter QParameter; typedef TriEdgeCollapseQuadricParameter QParameter;
typedef HelperType QH; 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<typename TriMeshType::VertexPointer> & WV(){ static std::vector<typename TriMeshType::VertexPointer> & WV(){
static std::vector<typename TriMeshType::VertexPointer> _WV; return _WV; static std::vector<typename TriMeshType::VertexPointer> _WV; return _WV;
} }
@ -166,141 +149,146 @@ public:
} }
inline bool IsFeasible(BaseParameterClass *_pp){ inline bool IsFeasible(BaseParameterClass *_pp){
QParameter *pp=(QParameter *)_pp; QParameter *pp=(QParameter *)_pp;
if(!pp->PreserveTopology) return true; if(!pp->PreserveTopology) return true;
bool res = ( EdgeCollapser<TriMeshType, VertexPair>::LinkConditions(this->pos) ); bool res = ( EdgeCollapser<TriMeshType, VertexPair>::LinkConditions(this->pos) );
if(!res) ++( TEC::FailStat::LinkConditionEdge() ); if(!res) ++( TEC::FailStat::LinkConditionEdge() );
return res; 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<TriMeshType,VertexPair>::Do(m, this->pos, ComputePosition(_pp)); // v0 is deleted and v1 take the new position
} }
void ComputePosition(BaseParameterClass *_pp)
// 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<typename TriMeshType::VertexPointer>::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; QParameter *pp=(QParameter *)_pp;
CoordType newPos = (this->pos.V(0)->P()+this->pos.V(1)->P())/2.0;
typename TriMeshType::VertexIterator vi; if(pp->OptimalPlacement==false)
typename TriMeshType::FaceIterator pf; newPos=this->pos.V(1)->P();
else
pp->CosineThr=cos(pp->NormalThrRad);
vcg::tri::UpdateTopology<TriMeshType>::VertexFace(m);
vcg::tri::UpdateFlags<TriMeshType>::FaceBorderFromVF(m);
if(pp->FastPreserveBoundary)
{ {
for(pf=m.face.begin();pf!=m.face.end();++pf) if((QH::Qd(this->pos.V(0)).Apply(newPos) + QH::Qd(this->pos.V(1)).Apply(newPos)) > 2.0*pp->QuadricEpsilon)
if( !(*pf).IsD() && (*pf).IsW() ) {
for(int j=0;j<3;++j) QuadricType q=QH::Qd(this->pos.V(0));
if((*pf).IsB(j)) q+=QH::Qd(this->pos.V(1));
{
(*pf).V(j)->ClearW(); Point3<QuadricType::ScalarType> x;
(*pf).V1(j)->ClearW(); if(pp->SVDPlacement)
} q.MinimumClosestToPoint(x,Point3d::Construct(newPos));
} else
q.Minimum(x);
if(pp->PreserveBoundary) 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<TriMeshType,VertexPair>::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(); typename TriMeshType::VertexIterator vi;
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)
for(vi=m.vert.begin();vi!=m.vert.end();++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<typename TriMeshType::VertexPointer>::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<TriMeshType>::VertexFace(m);
vcg::tri::UpdateFlags<TriMeshType>::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<FaceType> x; (*pf).V(j)->ClearW();
for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x){ (*pf).V1(j)->ClearW();
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()<x.V1()) && x.V1()->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()<x.V2()) && x.V2()->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(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<FaceType> 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()<x.V1()) && x.V1()->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()<x.V2()) && x.V2()->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) { // if the collapse is A-symmetric (e.g. u->v != v->u)
for(vi=m.vert.begin();vi!=m.vert.end();++vi) for(auto vi=m.vert.begin();vi!=m.vert.end();++vi)
if(!(*vi).IsD() && (*vi).IsRW()) if(!(*vi).IsD() && (*vi).IsRW())
{ {
vcg::face::VFIterator<FaceType> x; vcg::face::VFIterator<FaceType> x;
UnMarkAll(m); UnMarkAll(m);
for( x.F() = (*vi).VFp(), x.I() = (*vi).VFi(); x.F()!=0; ++ x) 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()))){
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)));
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()))){
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)));
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)?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 IsSymmetric(BaseParameterClass *_pp) {return ((QParameter *)_pp)->OptimalPlacement;}
static bool IsVertexStable(BaseParameterClass *_pp) {return !((QParameter *)_pp)->OptimalPlacement;} static bool IsVertexStable(BaseParameterClass *_pp) {return !((QParameter *)_pp)->OptimalPlacement;}
@ -315,60 +303,83 @@ public:
ScalarType ComputePriority(BaseParameterClass *_pp) ScalarType ComputePriority(BaseParameterClass *_pp)
{ {
QParameter *pp=(QParameter *)_pp; QParameter *pp=(QParameter *)_pp;
std::vector<CoordType> onVec; // vector with incident faces original normals
VertexType * v[2]; VertexType * v[2];
v[0] = this->pos.V(0); v[0] = this->pos.V(0);
v[1] = this->pos.V(1); v[1] = this->pos.V(1);
if(pp->NormalCheck){ // Compute maximal normal variation std::vector<CoordType> origNormalVec; // vector with incident faces original normals
// store the old normals for non-collapsed face in v0 if(pp->NormalCheck){ // Collect Original Normals
for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 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( x.V1()!=v[1] && x.V2()!=v[1] ) // skip faces with v1
onVec.push_back(TriangleNormal(*x.F()).Normalize()); origNormalVec.push_back(NormalizedTriangleNormal(*x.F()));
// store the old normals for non-collapsed face in v1
for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1 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( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0
onVec.push_back(TriangleNormal(*x.F()).Normalize()); 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<double>::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) //// Move the two vertexes into new position (storing the old ones)
CoordType OldPos0=v[0]->P(); CoordType OldPos0=v[0]->P();
CoordType OldPos1=v[1]->P(); CoordType OldPos1=v[1]->P();
CoordType newPos = ComputePosition(_pp); ComputePosition(_pp);
v[0]->P() = v[1]->P() = newPos; // 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<double>::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 ScalarType newQual = std::numeric_limits<ScalarType>::max(); //
int i=0; if(pp->QualityCheck){
double MinCos = std::numeric_limits<double>::max(); // minimo coseno di variazione di una normale della faccia for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0
// (e.g. max angle) Mincos varia da 1 (normali coincidenti) a if( x.V1()!=v[1] && x.V2()!=v[1] )
// -1 (normali opposte); newQual=std::min(newQual,QualityFace(*x.F()));
double MinQual = std::numeric_limits<double>::max(); for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1
for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0 if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0
if( x.V1()!=v[1] && x.V2()!=v[1] ) // skiping faces with v1 newQual=std::min(newQual,QualityFace(*x.F()));
{ }
if(pp->NormalCheck){
CoordType nn=NormalizedTriangleNormal(*x.F()); ScalarType newArea=0;
double ndiff=nn.dot(onVec[i++]); if(pp->AreaCheck){ // Collect Area
MinCos=std::min(MinCos,ndiff); for(VFIterator x(v[0]); !x.End(); ++x ) // for all faces in v0
} newArea += DoubleArea(*x.F());
if(pp->QualityCheck){ for(VFIterator x(v[1]); !x.End(); ++x ) // for all faces in v1
double qt= QualityFace(*x.F()); if( x.V1()!=v[0] && x.V2()!=v[0] ) // skip faces with v0
MinQual=std::min(MinQual,qt); 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
{
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);
}
}
QuadricType qq=QH::Qd(v[0]); QuadricType qq=QH::Qd(v[0]);
qq+=QH::Qd(v[1]); qq+=QH::Qd(v[1]);
@ -377,28 +388,40 @@ public:
assert(!math::IsNAN(QuadErr)); assert(!math::IsNAN(QuadErr));
// All collapses involving triangles with quality larger than <QualityThr> have no penalty; // All collapses involving triangles with quality larger than <QualityThr> have no penalty;
if(MinQual>pp->QualityThr) MinQual=pp->QualityThr; if(newQual>pp->QualityThr) newQual=pp->QualityThr;
if(pp->NormalCheck){ if(pp->NormalCheck){
// All collapses where the normal vary less than <NormalThr> (e.g. more than CosineThr) // All collapses where the normal vary less than <NormalThr> (e.g. more than CosineThr)
// have no penalty // have no particular penalty
if(MinCos>pp->CosineThr) MinCos=pp->CosineThr; 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); QuadErr= std::max(QuadErr,pp->QuadricEpsilon);
if(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; if( pp->UseVertexWeight ) QuadErr *= (QH::W(v[1])+QH::W(v[0]))/2;
ScalarType error; ScalarType error;
if(!pp->QualityCheck && !pp->NormalCheck) error = (ScalarType)(QuadErr); 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 / 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<ScalarType>::max();
if(pp->HardQualityCheck &&
(newQual < pp->HardQualityThr && newQual < origQual*0.9) )
error = std::numeric_limits<ScalarType>::max();
if(pp->HardNormalCheck)
if(CheckForFlip()) error = std::numeric_limits<ScalarType>::max();
// Restore old position of v0 and v1 // Restore old position of v0 and v1
v[0]->P()=OldPos0; v[0]->P()=OldPos0;
@ -408,17 +431,95 @@ public:
return this->_priority; return this->_priority;
} }
//
//static double MaxError() {return 1e100;} bool CheckForFlippedFaceOverVertex(VertexType *vp, ScalarType angleThrRad = math::ToRad(150.))
// {
std::map<VertexType *, CoordType> 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<VertexType *, CoordType> 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) inline void AddCollapseToHeap(HeapType & h_ret, VertexType *v0, VertexType *v1, BaseParameterClass *_pp)
{ {
QParameter *pp=(QParameter *)_pp; QParameter *pp=(QParameter *)_pp;
ScalarType maxAdmitErr = std::numeric_limits<ScalarType>::max();
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(v0,v1), this->GlobalMark(),_pp))); 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)){ if(!IsSymmetric(pp)){
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(v1,v0), this->GlobalMark(),_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 ) { for(VFIterator vfi(v[1]); !vfi.End(); ++vfi ) {
vfi.V1()->ClearV(); vfi.V1()->ClearV();
vfi.V2()->ClearV(); vfi.V2()->ClearV();
vfi.V1()->IMark() = this->GlobalMark();
vfi.V2()->IMark() = this->GlobalMark();
} }
// Second Loop // Second Loop
@ -457,8 +560,7 @@ public:
{ {
QParameter *pp=(QParameter *)_pp; QParameter *pp=(QParameter *)_pp;
QH::Init(); QH::Init();
// m.ClearFlags(); for(VertexIterator pv=m.vert.begin();pv!=m.vert.end();++pv)
for(VertexIterator pv=m.vert.begin();pv!=m.vert.end();++pv) // Azzero le quadriche
if( ! (*pv).IsD() && (*pv).IsW()) if( ! (*pv).IsD() && (*pv).IsW())
QH::Qd(*pv).SetZero(); QH::Qd(*pv).SetZero();
@ -487,9 +589,9 @@ public:
QuadricType bq; QuadricType bq;
// Border quadric record the squared distance from the plane orthogonal to the face and passing // Border quadric record the squared distance from the plane orthogonal to the face and passing
// through the edge. // through the edge.
borderPlane.SetDirection(facePlane.Direction() ^ ( (*fi).V1(j)->cP() - (*fi).V(j)->cP() ).normalized()); 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 if( (*fi).IsB(j) ) borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->BoundaryQuadricWeight )); // amplify border planes
else borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->BoundaryWeight/100.0)); // and consider much less quadric for quality else borderPlane.SetDirection(borderPlane.Direction()* (ScalarType)(pp->QualityQuadricWeight )); // and consider much less quadric for quality
borderPlane.SetOffset(borderPlane.Direction().dot((*fi).V(j)->cP())); borderPlane.SetOffset(borderPlane.Direction().dot((*fi).V(j)->cP()));
bq.ByPlane(borderPlane); bq.ByPlane(borderPlane);
@ -518,36 +620,24 @@ public:
} }
} }
CoordType ComputeMinimal()
{
// }
//
// CoordType ComputeMinimalOld()
//
//
//
//static void InitMesh(MESH_TYPE &m){
// pp->CosineThr=cos(pp->NormalThr);
// InitQuadric(m);
// //m.Topology();
// //OldInitQuadric(m,UseArea);
// }
//
CoordType ComputeMinimal()
{ {
typename TriMeshType::VertexType * v[2]; VertexType* &v0 = this->pos.V(0);
v[0] = this->pos.V(0); VertexType* &v1 = this->pos.V(1);
v[1] = this->pos.V(1); QuadricType q=QH::Qd(v0);
QuadricType q=QH::Qd(v[0]); q+=QH::Qd(v1);
q+=QH::Qd(v[1]);
Point3<QuadricType::ScalarType> x; Point3<QuadricType::ScalarType> x;
bool rt=q.Minimum(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. if(!rt) { // if the computation of the minimum fails we choose between the two edge points and the middle one.
Point3<QuadricType::ScalarType> x0=Point3d::Construct(v[0]->P()); Point3<QuadricType::ScalarType> x0=Point3d::Construct(v0->P());
Point3<QuadricType::ScalarType> x1=Point3d::Construct(v[1]->P()); Point3<QuadricType::ScalarType> x1=Point3d::Construct(v1->P());
x.Import((v[0]->P()+v[1]->P())/2); x.Import((v1->P()+v1->P())/2);
double qvx=q.Apply(x); double qvx=q.Apply(x);
double qv0=q.Apply(x0); double qv0=q.Apply(x0);
double qv1=q.Apply(x1); double qv1=q.Apply(x1);
@ -557,10 +647,8 @@ public:
return CoordType::Construct(x); return CoordType::Construct(x);
} }
//
//
}; };
} // namespace tri
} // namespace vcg } // namespace tri
} // namespace vcg
#endif #endif