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

Significantly refactored. Renamed a bunch of variables, removed useless parameters. Corrected various minor bugs. Changed the semantics of quality weighting (was bugged!)

This commit is contained in:
Paolo Cignoni 2016-03-24 14:17:34 +00:00
parent a58040cf9c
commit 73b0576edf
1 changed files with 233 additions and 293 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -20,56 +20,6 @@
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.14 2007/03/22 11:07:16 cignoni
Solved an issue related to different casting double-float between gcc 3 and gcc 4
Revision 1.13 2007/02/25 09:20:10 cignoni
Added Rad to the NormalThr Option and removed a bug in multiple exectuion of non optimal simplification (missing an isD check)
Revision 1.12 2007/01/19 09:13:14 cignoni
Added Finalize() method to the interface, corrected minor bugs on border preserving and postsimplification cleanup
Avoided double make_heap (it is done only in the local_optimization init)
Revision 1.11 2006/10/15 07:31:21 cignoni
typenames and qualifiers for gcc compliance
Revision 1.10 2006/10/09 20:12:55 cignoni
Heavyly restructured for meshlab inclusion. Now the access to the quadric elements are mediated by a static helper class.
Revision 1.9 2006/10/07 17:20:25 cignoni
Updated to the new style face->Normal() becomes Normal(face)
Revision 1.8 2005/10/02 23:19:36 cignoni
Changed the sign of the priority of a collapse. Now it is its the error as it should (and not -error)
Revision 1.7 2005/04/14 11:35:07 ponchio
*** empty log message ***
Revision 1.6 2005/01/19 10:35:28 cignoni
Better management of symmetric/asymmetric edge collapses
Revision 1.5 2004/12/10 01:07:15 cignoni
Moved param classes inside; added support for optimal placement and symmetric; added update heap also here (not only in the base class)
Revision 1.4 2004/11/23 10:34:23 cignoni
passed parameters by reference in many funcs and gcc cleaning
Revision 1.3 2004/10/25 07:07:56 ganovelli
A vcg.::Pos was used to implement the collapse type. CHanged
to vcg::Edge
Revision 1.2 2004/09/29 17:08:16 ganovelli
corrected error in -error (see localoptimization)
****************************************************************************/
#ifndef __VCG_TRIMESHCOLLAPSE_QUADRIC__
#define __VCG_TRIMESHCOLLAPSE_QUADRIC__
@ -81,6 +31,7 @@ corrected error in -error (see localoptimization)
#include<vcg/complex/algorithms/update/bounding.h>
#include<vcg/complex/algorithms/local_optimization/tri_edge_collapse.h>
#include<vcg/complex/algorithms/local_optimization.h>
#include<vcg/complex/algorithms/stat.h>
namespace vcg{
@ -147,9 +98,9 @@ public:
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; // all
double QualityThr; // Collapsed that generate faces with quality LOWER than this value are penalized. So
bool QualityWeight;
bool SafeHeapUpdate;
double QualityWeightFactor;
double ScaleFactor;
bool ScaleIndependent;
bool UseArea;
@ -168,9 +119,9 @@ public:
QuadricEpsilon =1e-15;
QualityCheck=true;
QualityQuadric=false;
QualityThr=.1;
QualityThr=.3; // higher the value -> better the quality of the accepted triangles
QualityWeight=false;
SafeHeapUpdate =false;
QualityWeightFactor=100.0;
ScaleFactor=1.0;
ScaleIndependent=true;
UseArea=true;
@ -186,27 +137,24 @@ class TriEdgeCollapseQuadric: public TriEdgeCollapse< TriMeshType, VertexPair, M
{
public:
typedef typename vcg::tri::TriEdgeCollapse< TriMeshType, VertexPair, MYTYPE > TEC;
// typedef typename TEC::EdgeType EdgeType;
typedef typename TriEdgeCollapse<TriMeshType, VertexPair, MYTYPE>::HeapType HeapType;
typedef typename TriEdgeCollapse<TriMeshType, VertexPair, MYTYPE>::HeapElem HeapElem;
typedef typename TriMeshType::CoordType CoordType;
typedef typename TriMeshType::ScalarType ScalarType;
typedef math::Quadric< double > QuadricType;
typedef typename TriMeshType::FaceType FaceType;
typedef typename TriMeshType::VertexType VertexType;
typedef typename TriMeshType::VertexIterator VertexIterator;
typedef typename TriMeshType::FaceIterator FaceIterator;
typedef typename vcg::face::VFIterator<FaceType> VFIterator;
typedef math::Quadric< double > QuadricType;
typedef TriEdgeCollapseQuadricParameter QParameter;
typedef HelperType QH;
// static QParameter & Params(){
// static QParameter p;
// return p;
// }
// puntatori ai vertici che sono stati messi non-w per preservare il boundary
static std::vector<typename TriMeshType::VertexPointer> & WV(){
static std::vector<typename TriMeshType::VertexPointer> _WV; return _WV;
};
}
inline TriEdgeCollapseQuadric(){}
@ -227,15 +175,23 @@ public:
return res;
}
void Execute(TriMeshType &m, BaseParameterClass *_pp)
CoordType ComputePosition(BaseParameterClass *_pp)
{
QParameter *pp=(QParameter *)_pp;
CoordType newPos;
if(pp->OptimalPlacement) newPos= static_cast<MYTYPE*>(this)->ComputeMinimal();
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, newPos); // v0 is deleted and v1 take the new position
EdgeCollapser<TriMeshType,VertexPair>::Do(m, this->pos, ComputePosition(_pp)); // v0 is deleted and v1 take the new position
}
@ -349,100 +305,105 @@ public:
static bool IsVertexStable(BaseParameterClass *_pp) {return !((QParameter *)_pp)->OptimalPlacement;}
///*
// Funzione principale di valutazione dell'errore del collasso.
// In pratica simula il collasso vero e proprio.
//
// Da ottimizzare il ciclo sulle normali (deve sparire on e si deve usare per face normals)
//*/
/** Evaluate the priority (error) for an edge collapse
*
* It simulate the collapse and compute the quadric error
* generated by this collapse. This error is weighted with
* - aspect ratio of involved triangles
* - normal variation
*/
ScalarType ComputePriority(BaseParameterClass *_pp)
{
QParameter *pp=(QParameter *)_pp;
ScalarType error;
typename vcg::face::VFIterator<FaceType> x;
std::vector<CoordType> on; // original normals
typename TriMeshType::VertexType * v[2];
std::vector<CoordType> 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
for(x.F() = v[0]->VFp(), x.I() = v[0]->VFi(); x.F()!=0; ++x ) // for all faces in v0
if(x.F()->V(0)!=v[1] && x.F()->V(1)!=v[1] && x.F()->V(2)!=v[1] ) // skip faces with v1
on.push_back(TriangleNormal(*x.F()).Normalize());
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
for(x.F() = v[1]->VFp(), x.I() = v[1]->VFi(); x.F()!=0; ++x ) // for all faces in v1
if(x.F()->V(0)!=v[0] && x.F()->V(1)!=v[0] && x.F()->V(2)!=v[0] ) // skip faces with v0
on.push_back(TriangleNormal(*x.F()).Normalize());
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());
}
//// Move the two vertexe into new position (storing the old ones)
//// Move the two vertexes into new position (storing the old ones)
CoordType OldPos0=v[0]->P();
CoordType OldPos1=v[1]->P();
if(pp->OptimalPlacement) { v[0]->P() = ComputeMinimal(); v[1]->P()=v[0]->P();}
else v[0]->P() = v[1]->P();
CoordType newPos = ComputePosition(_pp);
v[0]->P() = v[1]->P() = newPos;
//// Rescan faces and compute quality and difference between normals
int i;
double ndiff,MinCos = 1e100; // minimo coseno di variazione di una normale della faccia
int i=0;
double MinCos = std::numeric_limits<double>::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 qt, MinQual = 1e100;
CoordType nn;
for(x.F() = v[0]->VFp(), x.I() = v[0]->VFi(),i=0; x.F()!=0; ++x ) // for all faces in v0
if(x.F()->V(0)!=v[1] && x.F()->V(1)!=v[1] && x.F()->V(2)!=v[1] ) // skip faces with v1
double MinQual = std::numeric_limits<double>::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){
nn=TriangleNormal(*x.F()).Normalize();
ndiff=nn.dot(on[i++]);
if(ndiff<MinCos) MinCos=ndiff;
CoordType nn=NormalizedTriangleNormal(*x.F());
double ndiff=nn.dot(onVec[i++]);
MinCos=std::min(MinCos,ndiff);
}
if(pp->QualityCheck){
qt= QualityFace(*x.F());
if(qt<MinQual) MinQual=qt;
double qt= QualityFace(*x.F());
MinQual=std::min(MinQual,qt);
}
}
for(x.F() = v[1]->VFp(), x.I() = v[1]->VFi(),i=0; x.F()!=0; ++x ) // for all faces in v1
if(x.F()->V(0)!=v[0] && x.F()->V(1)!=v[0] && x.F()->V(2)!=v[0] ) // skip faces with v0
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){
nn=TriangleNormal(*x.F()).Normalize();
ndiff=nn.dot(on[i++]);
if(ndiff<MinCos) MinCos=ndiff;
CoordType nn=NormalizedTriangleNormal(*x.F());
double ndiff=nn.dot(onVec[i++]);
MinCos=std::min(MinCos,ndiff);
}
if(pp->QualityCheck){
qt= QualityFace(*x.F());
if(qt<MinQual) MinQual=qt;
double qt= QualityFace(*x.F());
MinQual=std::min(MinQual,qt);
}
}
QuadricType qq=QH::Qd(v[0]);
qq+=QH::Qd(v[1]);
Point3d tpd=Point3d::Construct(v[1]->P());
double QuadErr = pp->ScaleFactor*qq.Apply(tpd);
// All collapses involving triangles with quality larger than <QualityThr> has no penalty;
double QuadErr = pp->ScaleFactor*qq.Apply(Point3d::Construct(v[1]->P()));
assert(!math::IsNAN(QuadErr));
// All collapses involving triangles with quality larger than <QualityThr> have no penalty;
if(MinQual>pp->QualityThr) MinQual=pp->QualityThr;
if(pp->NormalCheck){
// All collapses where the normal vary less than <NormalThr> (e.g. more than CosineThr)
// have no penalty
if(MinCos>pp->CosineThr) MinCos=pp->CosineThr;
MinCos=(MinCos+1)/2.0; // Now it is in the range 0..1 with 0 very dangerous!
MinCos=fabs((MinCos+1)/2.0); // Now it is in the range 0..1 with 0 very dangerous!
}
if(QuadErr<pp->QuadricEpsilon) QuadErr=pp->QuadricEpsilon;
QuadErr= std::max(QuadErr,pp->QuadricEpsilon);
if(QuadErr <= pp->QuadricEpsilon)
{
QuadErr = - 1/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 / MinCos);
if( pp->QualityCheck && pp->NormalCheck) error = (ScalarType)(QuadErr / (MinQual*MinCos));
//Rrestore old position of v0 and v1
// Restore old position of v0 and v1
v[0]->P()=OldPos0;
v[1]->P()=OldPos1;
this->_priority = error;
return this->_priority;
}
@ -450,119 +411,90 @@ public:
//
//static double MaxError() {return 1e100;}
//
inline void UpdateHeap(HeapType & h_ret,BaseParameterClass *_pp)
inline void AddCollapseToHeap(HeapType & h_ret, VertexType *v0, VertexType *v1, BaseParameterClass *_pp)
{
QParameter *pp=(QParameter *)_pp;
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(v0,v1), this->GlobalMark(),_pp)));
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());
}
}
inline void UpdateHeap(HeapType & h_ret, BaseParameterClass *_pp)
{
this->GlobalMark()++;
VertexType *v[2];
v[0]= this->pos.V(0);
v[1]= this->pos.V(1);
v[1]->IMark() = this->GlobalMark();
// First loop around the remaining vertex to unmark visited flags
vcg::face::VFIterator<FaceType> vfi(v[1]);
while (!vfi.End()){
// First loop around the surviving vertex to unmark the Visit flags
for(VFIterator vfi(v[1]); !vfi.End(); ++vfi ) {
vfi.V1()->ClearV();
vfi.V2()->ClearV();
++vfi;
}
// Second Loop
vfi = face::VFIterator<FaceType>(v[1]);
while (!vfi.End())
{
assert(!vfi.F()->IsD());
for(VFIterator vfi(v[1]); !vfi.End(); ++vfi ) {
if( !(vfi.V1()->IsV()) && vfi.V1()->IsRW())
{
vfi.V1()->SetV();
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(vfi.V0(),vfi.V1()), this->GlobalMark(),_pp)));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric(pp)){
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(vfi.V1(),vfi.V0()), this->GlobalMark(),_pp)));
std::push_heap(h_ret.begin(),h_ret.end());
}
AddCollapseToHeap(h_ret,vfi.V0(),vfi.V1(),_pp);
}
if( !(vfi.V2()->IsV()) && vfi.V2()->IsRW())
{
vfi.V2()->SetV();
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(vfi.V0(),vfi.V2()),this->GlobalMark(),_pp)));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric(pp)){
h_ret.push_back( HeapElem(new MYTYPE(VertexPair(vfi.V2(),vfi.V0()), this->GlobalMark(),_pp) ) );
std::push_heap(h_ret.begin(),h_ret.end());
AddCollapseToHeap(h_ret,vfi.V2(),vfi.V0(),_pp);
}
}
if(pp->SafeHeapUpdate && vfi.V1()->IsRW() && vfi.V2()->IsRW() )
{
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(vfi.V1(),vfi.V2()),this->GlobalMark(),_pp)));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric(pp)){
h_ret.push_back(HeapElem(new MYTYPE(VertexPair(vfi.V2(),vfi.V1()), this->GlobalMark(),_pp)));
std::push_heap(h_ret.begin(),h_ret.end());
}
}
++vfi;
}
if(vfi.V1()->IsRW() && vfi.V2()->IsRW() )
AddCollapseToHeap(h_ret,vfi.V1(),vfi.V2(),_pp);
} // end second loop around surviving vertex.
}
static void InitQuadric(TriMeshType &m,BaseParameterClass *_pp)
{
{
QParameter *pp=(QParameter *)_pp;
typename TriMeshType::FaceIterator pf;
typename TriMeshType::VertexIterator pv;
int j;
QH::Init();
// m.ClearFlags();
for(pv=m.vert.begin();pv!=m.vert.end();++pv) // Azzero le quadriche
for(VertexIterator pv=m.vert.begin();pv!=m.vert.end();++pv) // Azzero le quadriche
if( ! (*pv).IsD() && (*pv).IsW())
QH::Qd(*pv).SetZero();
for(pf=m.face.begin();pf!=m.face.end();++pf)
if( !(*pf).IsD() && (*pf).IsR() )
if((*pf).V(0)->IsR() &&(*pf).V(1)->IsR() &&(*pf).V(2)->IsR())
for(FaceIterator fi=m.face.begin();fi!=m.face.end();++fi)
if( !(*fi).IsD() && (*fi).IsR() )
if((*fi).V(0)->IsR() &&(*fi).V(1)->IsR() &&(*fi).V(2)->IsR())
{
QuadricType q;
Plane3<ScalarType,false> p;
// Calcolo piano
p.SetDirection( ( (*pf).V(1)->cP() - (*pf).V(0)->cP() ) ^ ( (*pf).V(2)->cP() - (*pf).V(0)->cP() ));
// Se normalizzo non dipende dall'area
Plane3<ScalarType,false> facePlane;
facePlane.SetDirection( ( (*fi).V(1)->cP() - (*fi).V(0)->cP() ) ^ ( (*fi).V(2)->cP() - (*fi).V(0)->cP() ));
if(!pp->UseArea)
p.Normalize();
facePlane.Normalize();
facePlane.SetOffset( facePlane.Direction().dot((*fi).V(0)->cP()));
p.SetOffset( p.Direction().dot((*pf).V(0)->cP()));
QuadricType q;
q.ByPlane(facePlane);
// Calcolo quadrica delle facce
q.ByPlane(p);
// The basic < add face quadric to each vertex > loop
for(int j=0;j<3;++j)
if( (*fi).V(j)->IsW() )
QH::Qd((*fi).V(j)) += q;
for(j=0;j<3;++j)
if( (*pf).V(j)->IsW() )
for(int j=0;j<3;++j)
if( (*fi).IsB(j) || pp->QualityQuadric )
{
if(pp->QualityWeight)
q*=(*pf).V(j)->Q();
QH::Qd((*pf).V(j)) += q; // Sommo la quadrica ai vertici
}
Plane3<ScalarType,false> borderPlane;
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.SetOffset(borderPlane.Direction().dot((*fi).V(j)->cP()));
bq.ByPlane(borderPlane);
for(j=0;j<3;++j)
if( (*pf).IsB(j) || pp->QualityQuadric ) // Bordo!
{
Plane3<ScalarType,false> pb; // Piano di bordo
// Calcolo la normale al piano di bordo e la sua distanza
// Nota che la lunghezza dell'edge DEVE essere Normalizzata
// poiche' la pesatura in funzione dell'area e'gia fatta in p.Direction()
// Senza la normalize il bordo e' pesato in funzione della grandezza della mesh (mesh grandi non decimano sul bordo)
pb.SetDirection(p.Direction() ^ ( (*pf).V1(j)->cP() - (*pf).V(j)->cP() ).normalized());
if( (*pf).IsB(j) ) pb.SetDirection(pb.Direction()* (ScalarType)pp->BoundaryWeight); // amplify border planes
else pb.SetDirection(pb.Direction()* (ScalarType)(pp->BoundaryWeight/100.0)); // and consider much less quadric for quality
pb.SetOffset(pb.Direction().dot((*pf).V(j)->cP()));
q.ByPlane(pb);
if( (*pf).V (j)->IsW() ) QH::Qd((*pf).V (j)) += q; // Sommo le quadriche
if( (*pf).V1(j)->IsW() ) QH::Qd((*pf).V1(j)) += q;
if( (*fi).V (j)->IsW() ) QH::Qd((*fi).V (j)) += bq;
if( (*fi).V1(j)->IsW() ) QH::Qd((*fi).V1(j)) += bq;
}
}
@ -571,12 +503,20 @@ public:
vcg::tri::UpdateBounding<TriMeshType>::Box(m);
//Make all quadric independent from mesh size
pp->ScaleFactor = 1e8*pow(1.0/m.bbox.Diag(),6); // scaling factor
//pp->ScaleFactor *=pp->ScaleFactor ;
//pp->ScaleFactor *=pp->ScaleFactor ;
//printf("Scale factor =%f\n",pp->ScaleFactor );
//printf("bb (%5.2f %5.2f %5.2f)-(%5.2f %5.2f %5.2f) Diag %f\n",m.bbox.min[0],m.bbox.min[1],m.bbox.min[2],m.bbox.max[0],m.bbox.max[1],m.bbox.max[2],m.bbox.Diag());
}
}
if(pp->QualityWeight) // we map quality range into a squared 01 and than this into the 1..QualityWeightFactor range
{
float minQ, maxQ;
tri::Stat<TriMeshType>::ComputePerVertexQualityMinMax(m,minQ,maxQ);
for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
if( ! (*vi).IsD() && (*vi).IsW())
{
const double quality01squared = pow((vi->Q()-minQ)/(maxQ-minQ),2.0);
QH::Qd(*vi) *= 1.0 + quality01squared * (pp->QualityWeightFactor-1.0);
}
}
}
@ -616,7 +556,7 @@ public:
}
return CoordType::Construct(x);
}
}
//
//