vcglib/vcg/complex/local_optimization/tri_edge_collapse_quadric.h

565 lines
19 KiB
C
Raw Normal View History

2004-09-15 12:25:20 +02:00
/****************************************************************************
* VCGLib o o *
* Visual and Computer Graphics Library o o *
* _ O _ *
* Copyright(C) 2004 \/)\/ *
* Visual Computing Lab /\/| *
* ISTI - Italian National Research Council | *
* \ *
* All rights reserved. *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License (http://www.gnu.org/licenses/gpl.txt) *
* for more details. *
* *
****************************************************************************/
/****************************************************************************
History
$Log: not supported by cvs2svn $
Revision 1.7 2005/04/14 11:35:07 ponchio
*** empty log message ***
2005-04-14 13:35:09 +02:00
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)
2004-09-15 12:25:20 +02:00
****************************************************************************/
#ifndef __VCG_TRIMESHCOLLAPSE_QUADRIC__
#define __VCG_TRIMESHCOLLAPSE_QUADRIC__
#include<vcg/math/quadric.h>
#include<vcg/simplex/face/pos.h>
#include<vcg/complex/trimesh/update/flag.h>
#include<vcg/complex/trimesh/update/bounding.h>
#include<vcg/complex/local_optimization/tri_edge_collapse.h>
#include<vcg/complex/local_optimization.h>
namespace vcg{
namespace tri{
2004-09-15 12:25:20 +02:00
/**
This class describe Quadric based collapse operation.
Requirements:
Vertex
must have incremental mark
must have:
field QuadricType Q;
member
ScalarType W() const;
A per-vertex Weight that can be used in simplification
lower weight means that error is lowered,
standard: return W==1.0
void Merge(MESH_TYPE::vertex_type const & v);
Merges the attributes of the current vertex with the ones of v
(e.g. its weight with the one of the given vertex, the color ect).
Standard: void function;
Faces devono avere Shared Adjacency
durante la init serve FF per le quadriche di bordo
durante la semplificazione si usa VF
*/
template<class TriMeshType,class MYTYPE>
class TriEdgeCollapseQuadric: public TriEdgeCollapse< TriMeshType,MYTYPE>
{
public:
typedef typename vcg::tri::TriEdgeCollapse< TriMeshType, MYTYPE > TEC;
typedef typename TEC::EdgeType EdgeType;
2004-09-15 12:25:20 +02:00
typedef typename TriEdgeCollapse<TriMeshType, MYTYPE>::HeapType HeapType;
typedef typename TriEdgeCollapse<TriMeshType, MYTYPE>::HeapElem HeapElem;
typedef typename TriMeshType::CoordType CoordType;
typedef typename TriMeshType::ScalarType ScalarType;
typedef math::Quadric< double > QuadricType;
2004-09-15 12:25:20 +02:00
typedef typename TriMeshType::FaceType FaceType;
2005-04-14 13:35:09 +02:00
typedef typename TriMeshType::VertexType VertexType;
2004-09-15 12:25:20 +02:00
class QParameter
{
public:
double QualityThr; // all
double BoundaryWeight;
double NormalThr;
double CosineThr;
double QuadricEpsilon;
double ScaleFactor;
bool UseArea;
bool UseVertexWeight;
bool NormalCheck;
bool QualityCheck;
bool OptimalPlacement;
bool MemoryLess;
bool ComplexCheck;
bool ScaleIndependent;
//***********************
bool PreserveTopology;
bool PreserveBoundary;
bool MarkComplex;
bool FastPreserveBoundary;
bool SafeHeapUpdate;
};
static QParameter & Params(){static QParameter p; return p;}
2004-09-15 12:25:20 +02:00
enum Hint {
HNHasFFTopology = 0x0001, // La mesh arriva con la topologia ff gia'fatta
HNHasVFTopology = 0x0002, // La mesh arriva con la topologia bf gia'fatta
HNHasBorderFlag = 0x0004 // La mesh arriva con i flag di bordo gia' settati
};
static int & Hnt(){static int hnt; return hnt;} // the current hints
static void SetHint(Hint hn) { Hnt() |= hn; }
static void ClearHint(Hint hn) { Hnt()&=(~hn);}
static bool IsSetHint(Hint hn) { return (Hnt()&hn)!=0; }
// 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(const EdgeType &p, int i)
//:TEC(p,i){}
{
2005-04-14 13:35:09 +02:00
this->localMark = i;
this->pos=p;
this->_priority = ComputePriority();
}
2004-09-15 12:25:20 +02:00
inline bool IsFeasible(){
2005-04-14 13:35:09 +02:00
bool res = ( !Params().PreserveTopology || LinkConditions(this->pos) );
if(!res) ++( TriEdgeCollapse< TriMeshType,MYTYPE>::FailStat::LinkConditionEdge() );
return res;
}
2004-09-15 12:25:20 +02:00
void Execute(TriMeshType &m)
{ CoordType newPos;
if(Params().OptimalPlacement) newPos= ComputeMinimal();
2005-04-14 13:35:09 +02:00
else newPos=this->pos.V(1)->P();
this->pos.V(1)->q+=this->pos.V(0)->q;
int FaceDel=DoCollapse(this->pos, newPos); // v0 is deleted and v1 take the new position
2004-09-15 12:25:20 +02:00
m.fn-=FaceDel;
--m.vn;
}
static void Init(TriMeshType &m,HeapType&h_ret){
typename TriMeshType::VertexIterator vi;
typename TriMeshType::FaceIterator pf;
EdgeType av0,av1,av01;
Params().CosineThr=cos(Params().NormalThr);
2004-09-15 12:25:20 +02:00
if(!IsSetHint(HNHasVFTopology) ) vcg::tri::UpdateTopology<TriMeshType>::VertexFace(m);
if(Params().MarkComplex) {
vcg::tri::UpdateTopology<TriMeshType>::FaceFace(m);
vcg::tri::UpdateFlags<TriMeshType>::FaceBorderFromFF(m);
vcg::tri::UpdateTopology<TriMeshType>::VertexFace(m);
} // e' un po' piu' lenta ma marca i vertici complex
else
if(!IsSetHint(HNHasBorderFlag) )
vcg::tri::UpdateFlags<TriMeshType>::FaceBorderFromVF(m);
if(Params().FastPreserveBoundary)
{
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(Params().PreserveBoundary)
{
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);
// Initialize the heap with all the possible collapses
if(IsSymmetric())
{ // if the collapse is symmetric (e.g. u->v == v->u)
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
2004-09-15 12:25:20 +02:00
if((*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 )
{
2004-09-15 12:25:20 +02:00
assert(x.F()->V(x.I())==&(*vi));
if((x.V0()<x.V1()) && x.V1()->IsRW() && !x.V1()->IsV()){
x.V1()->SetV();
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType(x.V0(),x.V1()),TriEdgeCollapse< TriMeshType,MYTYPE>::GlobalMark() )));
2004-09-15 12:25:20 +02:00
}
if((x.V0()<x.V2()) && x.V2()->IsRW()&& !x.V2()->IsV()){
x.V2()->SetV();
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType(x.V0(),x.V2()),TriEdgeCollapse< TriMeshType,MYTYPE>::GlobalMark() )));
2004-09-15 12:25:20 +02:00
}
}
}
}
else
{ // if the collapse is A-symmetric (e.g. u->v != v->u)
2004-09-15 12:25:20 +02:00
for(vi=m.vert.begin();vi!=m.vert.end();++vi)
{
vcg::face::VFIterator<FaceType> x;
m.UnMarkAll();
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() && !m.IsMarked(x.F()->V1(x.I()))){
2005-04-14 13:35:09 +02:00
h_ret.push_back( HeapElem( new MYTYPE( EdgeType (x.V(),x.V1()),TriEdgeCollapse< TriMeshType,MYTYPE>::GlobalMark())));
}
if(x.V()->IsRW() && x.V2()->IsRW() && !m.IsMarked(x.F()->V2(x.I()))){
2005-04-14 13:35:09 +02:00
h_ret.push_back( HeapElem( new MYTYPE( EdgeType (x.V(),x.V2()),TriEdgeCollapse< TriMeshType,MYTYPE>::GlobalMark())));
}
}
}
}
2004-09-15 12:25:20 +02:00
make_heap(h_ret.begin(),h_ret.end());
}
static bool IsSymmetric() {return Params().OptimalPlacement;}
static bool IsVertexStable() {return !Params().OptimalPlacement;}
static void SetDefaultParams(){
Params().UseArea=true;
Params().UseVertexWeight=false;
Params().NormalCheck=false;
Params().NormalThr=M_PI/2;
2004-09-15 12:25:20 +02:00
Params().QualityCheck=true;
Params().QualityThr=.1;
Params().BoundaryWeight=.5;
Params().OptimalPlacement=true;
Params().ScaleIndependent=true;
Params().ComplexCheck=false;
Params().QuadricEpsilon =1e-15;
2004-09-15 12:25:20 +02:00
Params().ScaleFactor=1.0;
Params().PreserveTopology = false;
2004-09-15 12:25:20 +02:00
}
///*
// 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)
//*/
ScalarType ComputePriority() {
2004-09-15 12:25:20 +02:00
ScalarType error;
typename vcg::face::VFIterator<FaceType> x;
std::vector<CoordType> on; // original normals
typename TriMeshType::VertexType * v[2];
2005-04-14 13:35:09 +02:00
v[0] = this->pos.V(0);
v[1] = this->pos.V(1);
2004-09-15 12:25:20 +02:00
if(Params().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(x.F()->NormalizedNormal());
// 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(x.F()->NormalizedNormal());
}
//// Move the two vertexe into new position (storing the old ones)
CoordType OldPos0=v[0]->P();
CoordType OldPos1=v[1]->P();
if(Params().OptimalPlacement) { v[0]->P() = ComputeMinimal(); v[1]->P()=v[0]->P();}
else v[0]->P() = v[1]->P();
2004-09-15 12:25:20 +02:00
//// Rescan faces and compute quality and difference between normals
int i;
double ndiff,MinCos = 1e100; // 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
{
if(Params().NormalCheck){
nn=x.F()->NormalizedNormal();
ndiff=nn*on[i++];
if(ndiff<MinCos) MinCos=ndiff;
}
if(Params().QualityCheck){
qt= x.F()->QualityFace();
if(qt<MinQual) 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
{
if(Params().NormalCheck){
nn=x.F()->NormalizedNormal();
ndiff=nn*on[i++];
if(ndiff<MinCos) MinCos=ndiff;
}
if(Params().QualityCheck){
qt= x.F()->QualityFace();
if(qt<MinQual) MinQual=qt;
}
}
QuadricType qq=v[0]->q;
qq+=v[1]->q;
double QuadErr = Params().ScaleFactor*qq.Apply(v[1]->P());
// All collapses involving triangles with quality larger than <QualityThr> has no penalty;
if(MinQual>Params().QualityThr) MinQual=Params().QualityThr;
if(Params().NormalCheck){
// All collapses where the normal vary less than <NormalThr> (e.g. more than CosineThr)
// have no penalty
if(MinCos>Params().CosineThr) MinCos=Params().CosineThr;
MinCos=(MinCos+1)/2.0; // Now it is in the range 0..1 with 0 very dangerous!
}
if(QuadErr<Params().QuadricEpsilon) QuadErr=Params().QuadricEpsilon;
if( Params().UseVertexWeight ) QuadErr *= (v[1]->W()+v[0]->W())/2;
if(!Params().QualityCheck && !Params().NormalCheck) error = QuadErr;
if( Params().QualityCheck && !Params().NormalCheck) error = QuadErr / MinQual;
if(!Params().QualityCheck && Params().NormalCheck) error = QuadErr / MinCos;
if( Params().QualityCheck && Params().NormalCheck) error = QuadErr / (MinQual*MinCos);
//Rrestore old position of v0 and v1
v[0]->P()=OldPos0;
v[1]->P()=OldPos1;
this->_priority = error;
2005-04-14 13:35:09 +02:00
return this->_priority;
2004-09-15 12:25:20 +02:00
}
//
//static double MaxError() {return 1e100;}
//
inline void UpdateHeap(HeapType & h_ret)
{
2005-04-14 13:35:09 +02:00
this->GlobalMark()++;
VertexType *v[2];
2005-04-14 13:35:09 +02:00
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()){
vfi.V1()->ClearV();
vfi.V2()->ClearV();
++vfi;
}
// Second Loop
vfi = face::VFIterator<FaceType>(v[1]);
while (!vfi.End())
{
assert(!vfi.F()->IsD());
for (int j=0;j<3;j++)
{
if( !(vfi.V1()->IsV()) && vfi.V1()->IsRW())
{
vfi.V1()->SetV();
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V0(),vfi.V1()), this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric()){
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V1(),vfi.V0()), this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
}
}
if( !(vfi.V2()->IsV()) && vfi.V2()->IsRW())
{
vfi.V2()->SetV();
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V0(),vfi.V2()),this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric()){
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V2(),vfi.V0()), this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
}
}
if(Params().SafeHeapUpdate && vfi.V1()->IsRW() && vfi.V2()->IsRW() )
{
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V1(),vfi.V2()),this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
if(!IsSymmetric()){
2005-04-14 13:35:09 +02:00
h_ret.push_back(HeapElem(new MYTYPE(EdgeType (vfi.V2(),vfi.V1()), this->GlobalMark())));
std::push_heap(h_ret.begin(),h_ret.end());
}
}
}
++vfi;
}
}
2004-09-15 12:25:20 +02:00
static void InitQuadric(TriMeshType &m)
{
typename TriMeshType::FaceIterator pf;
typename TriMeshType::VertexIterator pv;
int j;
// m.ClearFlags();
for(pv=m.vert.begin();pv!=m.vert.end();++pv) // Azzero le quadriche
if( ! (*pv).IsD() && (*pv).IsW())
(*pv).q.Zero();
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())
{
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
if(!Params().UseArea)
p.Normalize();
2004-09-15 12:25:20 +02:00
p.SetOffset( p.Direction() * (*pf).V(0)->cP());
// Calcolo quadrica delle facce
q.ByPlane(p);
for(j=0;j<3;++j)
if( (*pf).V(j)->IsW() ) (*pf).V(j)->q += q; // Sommo la quadrica ai vertici
for(j=0;j<3;++j)
if( (*pf).IsB(j)) // 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() ).Normalize());
pb.SetDirection(pb.Direction()* Params().BoundaryWeight); // amplify border planes
pb.SetOffset(pb.Direction() * (*pf).V(j)->cP());
q.ByPlane(pb);
if( (*pf).V (j)->IsW() ) (*pf).V (j)->q += q; // Sommo le quadriche
if( (*pf).V1(j)->IsW() ) (*pf).V1(j)->q += q;
}
}
if(Params().ScaleIndependent)
{
vcg::tri::UpdateBounding<TriMeshType>::Box(m);
//Make all quadric independent from mesh size
Params().ScaleFactor = 1e8*pow(1.0/m.bbox.Diag(),6); // scaling factor
//Params().ScaleFactor *=Params().ScaleFactor ;
//Params().ScaleFactor *=Params().ScaleFactor ;
//printf("Scale factor =%f\n",Params().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(Params().ComplexCheck)
{
// secondo loop per diminuire quadriche complex (se non c'erano i complex si poteva fare in un giro solo)
//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(j=0;j<3;++j)
// if((*pf).IsCF(j)) // Complex!
// {
// if( (*pf).V (j)->IsW() ) (*pf).V (j)->q *= 0.01; // Scalo le quadriche
// if( (*pf).V1(j)->IsW() ) (*pf).V1(j)->q *= 0.01;
// }
// }
}
}
//
//
//
//
//
//
//static void InitMesh(MESH_TYPE &m){
// Params().CosineThr=cos(Params().NormalThr);
// InitQuadric(m);
// //m.Topology();
// //OldInitQuadric(m,UseArea);
// }
//
CoordType ComputeMinimal()
{
typename TriMeshType::VertexType * v[2];
2005-04-14 13:35:09 +02:00
v[0] = this->pos.V(0);
v[1] = this->pos.V(1);
2004-09-15 12:25:20 +02:00
QuadricType q=v[0]->q;
q+=v[1]->q;
CoordType x;
bool rt=q.Minimum(x);
if(!rt) {
x=(v[0]->P()+v[1]->P())/2;
double qvx=q.Apply(x);
double qv0=q.Apply(v[0]->P());
double qv1=q.Apply(v[1]->P());
if(qv0<qvx) x=v[0]->P();
if(qv1<qvx && qv1<qv0) x=v[1]->P();
}
return x;
}
//
//
};
} // namespace tri
} // namespace vcg
#endif