added variational shape approximation partition of the mesh both

for triangle meshes and for vertex meshes
This commit is contained in:
ganovelli 2011-06-14 17:10:57 +00:00
parent a815890b15
commit 09e0ccc62c
1 changed files with 921 additions and 0 deletions

View File

@ -0,0 +1,921 @@
#ifndef VCGLIB_REGION_GROWING_VERTEX_
#define VCGLIB_REGION_GROWING_VERTEX_
#include <list>
#include <vector>
#include <deque>
#include <vcg/complex/allocate.h>
#include <vcg/complex/append.h>
#include <vcg/space/fitting3.h>
#include <vcg/space/color4.h>
//#include <vcg/complex/algorithms/update/quality.h>
namespace vcg{
namespace tri{
/**
Element and associated error
*/
template <class ElemType>
struct ElemError{
ElemType * f;
float val;
ElemError(){};
ElemError(ElemType * _f,double _val):f(_f),val(_val){}
const bool operator <(const ElemError & o) const {return val < o.val;}
};
/**
Base class to define a region of a mesh
*/
template<class MeshType, class ElemType>
struct RegionBase{
RegionBase():isd(false),size(0){}
typedef typename MeshType MeshType;
typedef typename ElemType::CoordType CoordType;
typedef typename ElemType::CoordType::ScalarType ScalarType;
typedef typename std::list<RegionBase*>::iterator AdjIterator;
typedef ElemError<ElemType> ElemError;
typedef typename std::vector<ElemType*>::iterator ElemIte;
// elem belonging to the region
std::vector<ElemType*> elems;
// adjacent regions
std::list<RegionBase*> adj,nx_adj;
ElemError worst;
int id,size;
int cleansize;
ScalarType epsilon;
ScalarType approx_err;
ScalarType approx_var;
int changed;
bool isd;
void Clean(){
adj.sort();
adj.unique();
nx_adj.clear();
}
void Connect( RegionBase * tr){
adj.push_back(tr);
}
void UpdateError( ){
ElemIte vi;
float err=0,var=0,max_err=-1.0;
for(vi=elems.begin(); vi != elems.end(); ++vi)
{
float e =Evaluate(**vi);
err+=e;
if(e > max_err) {
worst.val = e;
worst.f = *vi;
max_err = e;
}
}
approx_err =err;
var = 0.0;
approx_var = var;
}
ScalarType Mergeable( RegionBase & tr){
/*
two regions are mergeable if:
- they have the same planes
- they have an extreme in common
or
- one is surrounded by the other which is much bigger
*/
//if(!Similar(this->p,tr.p))
// return false;
vcg::Plane3<ScalarType> plane;
ScalarType fit_error,new_error=0;
if( (this->adj.size() == 1) && (tr.elems.size()/float(this->elems.size()) > 5) ||
(tr.adj.size() == 1) && (this->elems.size()/float(tr.elems.size()) > 5)
)
return true;
std::vector<ElemType*> merged;
merged.insert(merged.end(),this->elems.begin(),this->elems.end());
merged.insert(merged.end(),tr.elems.begin(),tr.elems.end());
Fit(merged,plane,fit_error);
for(std::vector<ElemType*> ::iterator vi= merged.begin(); vi != merged.end(); ++vi)
new_error+= Evaluate(**vi);
new_error/=merged.size();
return new_error < this->approx_err * (elems.size()/float(merged.size()))+tr.approx_err * (tr.elems.size()/float(merged.size())) ;
}
/// pure virtual function to define how to fit the primitive onto this region
virtual void Fit(std::vector<ElemType*>& verts, vcg::Plane3f & plane, float & err) = 0;
/// pure virtual function to define how to fit the evaluate the actio of adding and element "e" to this region
virtual ScalarType Evaluate( ElemType & e) = 0;
/// pure virtual function to init the primitive with the element
virtual void Init(ElemType*) = 0;
};
/**
Base class to define a region grower, i.e. a floo algorithm to partion the mesh in regions
*/
template < class RegionType>
struct RegionGrower{
typedef typename RegionType RegionType;
typedef typename RegionType::MeshType MeshType;
typedef typename RegionType::ElemType ElemType;
typedef typename RegionType::ScalarType ScalarType;
typedef typename std::list<RegionType> ::iterator TriRegIterator;
typedef typename std::list<RegionBase<MeshType,ElemType>* >::iterator AdjIterator;
typedef typename RegionType::ElemError ElemError;
struct CandiElem{
ElemType * f;
float val;
RegionType * r;
CandiElem(){};
CandiElem(ElemType * _f,double _val,RegionType * _r):f(_f),val(_val),r(_r){}
const bool operator <(const CandiElem & o) const {return val < o.val;}
};
struct ErrorEval {
void Add(const float & v){
{
if(v < samples[i_min]) i_min = 0; else {++i_min;
if(v > samples[i_max]) i_max = 0;else ++i_max;
}
if(i_min == samples.size()) {i_min = 0; for(int i= 0; i < samples.size()-1; ++i) if(samples[i]<samples[i_min]) i_min = i; ++i_min; }
if(i_max == samples.size()) {i_max = 0; for(int i= 0; i < samples.size()-1; ++i) if(samples[i]>samples[i_max]) i_max = i; ++i_max;}
}
samples.pop_back();samples.push_front(v);
boxes.pop_back(); boxes.push_front(vcg::Point2f(samples[i_min],samples[i_max]));
++ns ;
}
float BoxOverlap(){
float maxsize = std::max( boxes.back()[1]-boxes.back()[0], (*boxes.begin())[1]-(*boxes.begin())[0]);
float overlap = std::max(0.f, std::min(boxes.back()[1],(*boxes.begin())[1])-std::max(boxes.back()[0],(*boxes.begin())[0]));
assert(overlap <= maxsize);
return (maxsize > 0.f)?overlap / maxsize:0.0;
}
float RelativeDecrease(){
if(ns<4) return std::numeric_limits<float>::max();
return fabs(samples[std::min<unsigned int>(samples.size()-1,ns)]-samples[0]);
}
void Init(int n ){
samples.clear();
boxes.clear();
for(int i = 0 ; i < n; ++i) samples.push_front(std::numeric_limits<float>::max());
for(int i = 0 ; i < n; ++i) boxes.push_front(vcg::Point2f(n,n-i));
i_max = i_min = 0;
ns = 0;
}
private:
int i_min,i_max; // index of min and max element in the queue
std::deque<float> samples;
std::deque<vcg::Point2f> boxes;
int ns;
};
RegionGrower():lastAdded(NULL),ave_error(-1),ave_var(-1),n_steps(0){}
std::list<RegionType> regions;
std::vector<RegionType*> toTunnel;
std::vector<RegionType*> workingset;
int n_faces,target_max_regions;
ElemType * lastAdded;
ScalarType ave_error // average error (over single regions' error)
,ave_var// average variance (over single regions' error)
,changed // faces that have changed from previous step ([0..1))
,err
,target_error; // taget error
ElemType * worst;
ScalarType worst_err ;
MeshType * m;
ErrorEval erroreval;
unsigned int n_steps;// number of steps done
std::vector<CandiElem> elemheap;
std::vector<ElemError > elemerr;
void Init(MeshType & mesh, int n_seeds,int max_regions, float max_err){
erroreval.Init(10);
m = &mesh;
target_error = mesh.bbox.Diag()*max_err;
target_max_regions = max_regions;
regions.clear();
}
/// add a region
void AddRegion(const RegionType & r){regions.push_back(r);}
/// remove a region
void DeleteRegion(const typename std::list<RegionType>::iterator & ri){std::remove(ri);}
void PushHeap(std::vector<ElemType*> & candi, RegionType & r){
typename std::vector<ElemType*>::iterator ci;
for(ci = candi.begin(); ci != candi.end(); ++ci)
{
elemheap.push_back(CandiElem( *ci,-r.Evaluate(*(*ci)), &r));
push_heap(elemheap.begin(),elemheap.end());
}
}
/// make two regions adjacent
void Connect(RegionType * r1,RegionType * r2){
assert(r1!=r2);
r1->Connect(r2);
r2->Connect(r1);
}
/// initialize a region
void CreateRegion(ElemType * vi){
AddRegion(RegionType());
RegionType & tr =regions.back();
AddElemToRegion(tr,vi);
tr.Refit();
tr.color = vcg::Color4b::Scatter(2000,(int)regions.size());
}
void AddElemToRegion( RegionType & r, ElemType * v){
r.elems.push_back(v);
B(v) = (RegionType*) &r;
if(B(v)!=B_old(v)) ++r.changed;
B_old(v) = B(v);
++r.size;
}
/// for each region take the candidates and fill in elemheap
void Refill( ){
elemheap.clear();
typename std::list<RegionType>::iterator ri;
std::vector<ElemType*> candi;
for(ri = regions.begin(); ri != regions.end(); ++ri) if(!(*ri).isd)
{
candi.clear();
Candidates((*ri),candi);
PushHeap(candi,*ri);
}
std::make_heap(elemheap.begin(),elemheap.end());
}
/// check if the iteration converged
bool IsRelaxed(){
ScalarType _ave_error=0;
ScalarType _ave_var= 0;
ScalarType _changed = 0.0;
int nr=0;
typename std::list<RegionType>::iterator ri;
worst=NULL;;
worst_err = -1;
for(ri = regions.begin(); ri != regions.end(); ++ri) if(!(*ri).isd){
++nr;
n_faces+=(*ri).elems.size();
(*ri).UpdateError();
_ave_error+=(*ri).approx_err;
_ave_var+=(*ri).approx_var;
_changed+=(*ri).changed;
(*ri).changed=0;
if((*ri).worst.val*(*ri).size > worst_err){
worst = (*ri).worst.f;
worst_err = (*ri).worst.val*(*ri).size;
}
}
_ave_error/=nr;
_ave_var/=nr;
_changed/=n_faces;
n_faces = 0;
erroreval.Add(_ave_error);
printf("Err: %f ov: %f Dec: %f \n",_ave_error,erroreval.BoxOverlap(),erroreval.RelativeDecrease());
return (erroreval.BoxOverlap() > 0.95) || (erroreval.RelativeDecrease() < 0.01*_ave_error);
}
/// merge two regions
void Merge(RegionType & r0,RegionBase<typename RegionType::MeshType,typename RegionType::ElemType> & r1){
assert(!r1.isd);
typename RegionType::ElemIte vi;
AdjIterator ai;
for(vi = r1.elems.begin();vi != r1.elems.end(); ++vi)
AddElemToRegion(r0,(*vi));
for(ai= r1.adj.begin(); ai != r1.adj.end();++ai)
if( !(*ai)->isd && (*ai)!=&r0)
r0.nx_adj.push_back(*ai);
r1.elems.clear();
r1.isd = true;
r0.Refit();
}
/// clean degeneracies: look for very small regions in the middle of a single regions (and similar cases)
bool IsToTunnel(RegionType * r){
return (r->elems.size()<2);
}
/// clean degeneracies: look for very small regions in the middle of a single regions (and similar cases)
void ComputeToTunnel(){
typename std::list<RegionType>::iterator ri;
for(ri = regions.begin(); ri != regions.end(); ++ri)
if(IsToTunnel(&(*ri)))
toTunnel.push_back(&(*ri));
}
void GrowStepOnce(){
if(elemheap.empty()) return;
CandiElem cf;
std::vector<ElemType*> toAdd;
std::pop_heap(elemheap.begin(),elemheap.end());
cf = elemheap.back();
//printf("err:%f\n",cf.val);
elemheap.pop_back();
if (B(cf.f)==NULL){
AddElemToRegion( *cf.r,cf.f); // ads to region
lastAdded = &*cf.f;
for(int i =0; i < cf.f->SizeNeigh();++i)
if((B(cf.f ->Neigh(i)) == NULL) )
toAdd.push_back(cf.f->Neigh(i));
else
if(B(cf.f ->Neigh(i)) != B(cf.f))
Connect((RegionType*)B(cf.f->Neigh(i)),(RegionType*)cf.r);
PushHeap(toAdd,*cf.r);
}
else
{
if ( B(cf.f) != (RegionType*)(cf.r) )
Connect((RegionType*)B(cf.f),(RegionType*)cf.r);
}
}
/// Execute and iteration
void GrowStep(){
CandiElem cf;
typename std::list<RegionType>::iterator ri;
n_steps++;
for(ri = regions.begin(); ri != regions.end(); ++ri) if(!(*ri).isd)
(*ri).Clean();
/*printf("Heap size %d\n", elemheap.size());*/
while(!elemheap.empty() ){
std::vector<ElemType*> toAdd;
std::pop_heap(elemheap.begin(),elemheap.end());
cf = elemheap.back();
//printf("err:%f\n",cf.val);
elemheap.pop_back();
if (B(cf.f)==NULL){
AddElemToRegion( *cf.r,cf.f); // ads to region
lastAdded = &*cf.f;
for(int i =0; i < cf.f->SizeNeigh();++i)
if((B(cf.f ->Neigh(i)) == NULL) )
toAdd.push_back(cf.f->Neigh(i));
else
if(B(cf.f ->Neigh(i)) != B(cf.f))
Connect((RegionType*)B(cf.f->Neigh(i)),(RegionType*)cf.r);
PushHeap(toAdd,*cf.r);
}
else
{
if ( B(cf.f) != (RegionType*)(cf.r) )
Connect((RegionType*)B(cf.f),(RegionType*)cf.r);
}
}
int h = (int) elemheap.size();
//printf("----> %d\n",h);
}
/// try to merge similare adjacent regions
unsigned int MergeStep(){
TriRegIterator tri;
unsigned int merged = 0;
typename RegionType::AdjIterator ai;
for(tri = regions.begin(); tri != regions.end(); ++tri) if(!(*tri).isd) (*tri).Clean();
for(tri = regions.begin(); tri != regions.end(); ++tri)if(!(*tri).isd)
for(ai = (*tri).adj.begin(); ai != (*tri).adj.end(); ++ai) if(!(*ai)->isd)
if((*tri).Mergeable(*(*ai))){
Merge((*tri),*(*ai));
merged++;
}
for(tri = regions.begin(); tri != regions.end(); ++tri){
for(ai = (*tri).nx_adj.begin(); ai != (*tri).nx_adj.end();++ai)
if(!(*ai)->isd)
(*tri).adj.push_back(*ai);
(*tri).adj.sort();
(*tri).adj.unique();
}
return merged;
}
/// re-init all for a new itearation
bool Restart(){
std::vector<ElemType*> candi;
TriRegIterator ri;
elemheap.clear();
printf("#regions: %d",this->regions.size());
ComputeToTunnel();
if(IsRelaxed()){
printf("R worst_err:%f, target_error %f \n",worst_err ,target_error);
if( (worst_err <= target_error) || (regions.size() >= target_max_regions))
return false;
else
{ printf("Add Region \n");
erroreval.Init(10);
ElemError wrs;
wrs.f = worst;
wrs.val = worst_err;
printf("worst triangle error %f\n",worst_err);
CreateRegion(wrs.f);// CreateRegion changes wr->r
// reset state variables
ave_error=-1;
ave_var=-1;
err=0.0;
changed=0;
}
printf("\n");
}
//
printf("tunnelling %d regions\n",toTunnel.size());
std::vector<RegionType*>::iterator it;
for(it = toTunnel.begin(); it != toTunnel.end(); ++it)
{
RegionType * r = *it;
typename RegionType::ElemType * e;
typename RegionType::ElemIte vi;
for(vi = r->elems.begin(); vi != r->elems.end(); ++vi) B(*vi) = NULL;
r->elems.clear();
do{e = RandomElem();} while( (B(e)!=NULL) && !IsToTunnel(B(e)));
AddElemToRegion(*r,e);
r->Refit();
//r->color = vcg::Color4b::Scatter(2000,(int)regions.size());
//r->Init(e);
Candidates(*r ,candi); // take the candidatees
PushHeap(candi,(*r )); // put the faces on to the heap
}
toTunnel.clear();
//
for(ri = regions.begin(); ri != regions.end(); )
if((*ri).isd)
ri = regions.erase(ri);
else
++ri;
for(ri = regions.begin(); ri != regions.end(); ++ri)
{
candi.clear();
(*ri).Refit(); // fit a plane to the region
Restart(*ri); // clear stuff in the region, move the seed to the best fitting to the plabne
Candidates(*ri,candi); // take the (three) faces candidatees
PushHeap(candi,(*ri)); // put the faces on to the heap
}
return true;
}
/// re-init a specific region
void Restart(RegionType &r){
if(!r.elems.empty()){
r.Refit();
r.size=0;
//float b_err = r.Evaluate(*(*r.elems.begin())),err;
ElemType* b_vert =(*r.elems.begin());
typename RegionType::ElemIte vi;
//for( vi = r.elems.begin(); vi != r.elems.end(); ++vi)
//{
// err = r.Evaluate(**vi);
// if(err < b_err)
// {
// b_err = err;
// b_vert = *vi;
// }
//}
b_vert = r.elems[r.elems.size()/2];
for( vi = r.elems.begin(); vi != r.elems.end(); ++vi) B(*vi) = NULL;
r.elems.clear();
r.adj.clear();
AddElemToRegion(r,b_vert);
r.Init(b_vert);
}
}
// after init, run all the process
void MakeCharts(){
this->Refill();
while(this->Restart()){
//do{
this->GrowStep();
//} while(Restart());
}
//while(this->MergeStep());
}
/// pick a random element
virtual typename RegionType::ElemType * RandomElem() = 0;
/// given an element returns the region to which it is associated (or NULL )
virtual RegionType * &B(ElemType *) = 0;
/// given an element returns the region to which it was associated (or NULL )
virtual RegionType * &B_old(ElemType *) = 0;
/// return in c the canidates to be added to the region r
virtual void Candidates(RegionType & r, std::vector< typename RegionType::ElemType*> & c) = 0;
};
template < class RegionType>
struct RegionGrowerVertex: public RegionGrower<RegionType> {
typedef typename RegionType RegionType;
typedef typename RegionType::MeshType MeshType;
typedef typename RegionType::ElemType ElemType;
typedef typename RegionType::ScalarType ScalarType;
typedef typename std::list<RegionType> ::iterator TriRegIterator;
typename MeshType:: template PerVertexAttributeHandle<RegionType*> attr_r;
typename MeshType:: template PerVertexAttributeHandle<RegionType*> attr_r_old;
RegionType * &B(ElemType *e){return (RegionType *)attr_r[e];};
RegionType * &B_old(ElemType *e){return (RegionType *)attr_r_old[e];}
void Init(MeshType & mesh, int n_seeds,int max_regions, float max_err){
RegionGrower<RegionType>::Init(mesh, n_seeds,max_regions,max_err);
attr_r = vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<RegionType*> (*m,"r");
if(!vcg::tri::Allocator<MeshType>::IsValidHandle(*m,attr_r))
attr_r = vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<RegionType*> (*m,"r");
attr_r_old = vcg::tri::Allocator<MeshType>::template GetPerVertexAttribute<RegionType*> (*m,"r_old");
if(!vcg::tri::Allocator<MeshType>::IsValidHandle(*m,attr_r_old))
attr_r_old = vcg::tri::Allocator<MeshType>::template AddPerVertexAttribute<RegionType*> (*m,"r_old");
for(int i = 0; i < m->vert.size(); ++i){
attr_r[i] = NULL;
attr_r_old[i] = NULL;
if( (i%(m->vn/n_seeds))==0)
CreateRegion(&m->vert[i]);
}
}
void Candidates(RegionType & r, std::vector< typename RegionType::ElemType*> & c){
typename RegionType::VertexIterator vi;
for(vi = r.elems.begin(); vi!= r.elems.end(); ++vi)
for(int i =0; i < (*vi)->SizeNeigh();++i)
if((attr_r[(*vi)->Neigh(i)] == NULL) )
c.push_back((*vi)->Neigh(i));
}
///Export the regions as a collection of meshes
void Export(std::vector<MeshType*> & submeshes){
MeshType::PerVertexAttributeHandle<char> sel = Allocator<MeshType>::AddPerVertexAttribute<char>(*m,"sel");
assert(Allocator<MeshType>::IsValidHandle<char>(*m,sel));
for( MeshType::VertexIterator vi = m->vert.begin();vi!= m->vert.end(); ++vi)
sel[(*vi)] = (*vi).IsS();
for(TriRegIterator ti = regions.begin(); ti != regions.end(); ++ti){
submeshes.push_back(new MeshType());
for(std::vector<MeshType::VertexPointer>::iterator vi = (*ti).elems.begin();vi!=(*ti).elems.end(); ++vi)
{(*vi)->SetS(); }
Append<MeshType,MeshType>::Mesh( *submeshes.back(),*m,true);
for(std::vector<MeshType::VertexPointer>::iterator vi = (*ti).elems.begin();vi!=(*ti).elems.end(); ++vi)
{(*vi)->ClearS();}
}
for( MeshType::VertexIterator vi = m->vert.begin();vi!= m->vert.end(); ++vi)
if(sel[(*vi)]) (*vi).SetS(); else (*vi).ClearS();
}
typename MeshType::VertexType * RandomElem(){
int id = std::max<int>(0,std::min<int>((rand()/float(RAND_MAX))*m->vert.size()-1,m->vert.size()-1));
return &m->vert[id];
}
};
template < class RegionType>
struct RegionGrowerFace: public RegionGrower<RegionType> {
typedef typename RegionType RegionType;
typedef typename RegionType::MeshType MeshType;
typedef typename RegionType::ElemType ElemType;
typedef typename RegionType::ScalarType ScalarType;
typedef typename std::list<RegionType> ::iterator TriRegIterator;
typename MeshType:: template PerFaceAttributeHandle<RegionType*> attr_r;
typename MeshType:: template PerFaceAttributeHandle<RegionType*> attr_r_old;
RegionType * &B(ElemType *e){return (RegionType *)attr_r[e];};
RegionType * &B_old(ElemType *e){return (RegionType *)attr_r_old[e];}
void Init(MeshType & mesh, int n_seeds,int max_regions, float max_err){
RegionGrower<RegionType>::Init(mesh, n_seeds,max_regions,max_err);
attr_r = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<RegionType*> (*m,"r");
if(!vcg::tri::Allocator<MeshType>::IsValidHandle(*m,attr_r))
attr_r = vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<RegionType*> (*m,"r");
attr_r_old = vcg::tri::Allocator<MeshType>::template GetPerFaceAttribute<RegionType*> (*m,"r_old");
if(!vcg::tri::Allocator<MeshType>::IsValidHandle(*m,attr_r_old))
attr_r_old = vcg::tri::Allocator<MeshType>::template AddPerFaceAttribute<RegionType*> (*m,"r_old");
for(int i = 0; i < m->face.size(); ++i){
attr_r[i] = NULL;
attr_r_old[i] = NULL;
if( (i%(m->fn/n_seeds))==0)
CreateRegion(&m->face[i]);
}
}
void Candidates(RegionType & r, std::vector< typename RegionType::ElemType*> & c){
typename RegionType::FaceIterator fi;
for(fi = r.elems.begin(); fi!= r.elems.end(); ++fi)
for(int i =0; i < (*fi)->VN();++i)
if((attr_r[(*fi)->FFp(i)] == NULL) )
c.push_back((*fi)->FFp(i));
}
///Export the regions as a collection of meshes
void Export(std::vector<MeshType*> & submeshes){
MeshType::PerFaceAttributeHandle<char> sel = Allocator<MeshType>::AddPerFaceAttribute<char>(*m,"sel");
assert(Allocator<MeshType>::IsValidHandle<char>(*m,sel));
for( MeshType::FaceIterator fi = m->face.begin();fi!= m->face.end(); ++fi)
sel[(*fi)] = (*fi).IsS();
for(TriRegIterator ti = regions.begin(); ti != regions.end(); ++ti){
submeshes.push_back(new MeshType());
for(std::vector<MeshType::FacePointer>::iterator fi = (*ti).face.begin();fi!=(*ti).face.end(); ++fi)
{(*fi)->SetS(); for(unsigned int i = 0; i < 3; ++i) (*fi)->V(i)->SetS();}
Append<MeshType,MeshType>::Mesh( *submeshes.back(),*m,true);
for(std::vector<MeshType::FacePointer>::iterator fi = (*ti).face.begin();fi!=(*ti).face.end(); ++fi)
{(*fi)->ClearS();for(unsigned int i = 0; i < 3; ++i) (*fi)->V(i)->ClearS();}
}
for( MeshType::FaceIterator fi = m->face.begin();fi!= m->face.end(); ++fi)
if(sel[(*fi)]) (*fi).SetS(); else (*fi).ClearS();
}
typename MeshType::FaceType * RandomElem(){
int id = std::max<int>(0,std::min<int>((rand()/float(RAND_MAX))*m->face.size()-1,m->face.size()-1));
return &m->face[id];
}
};
/**
Planar region made of vertices
*/
template<class MeshType>
struct PlanarRegionVertex: public RegionBase<MeshType,typename MeshType::VertexType> {
public:
typedef PlanarRegionVertex RegionType;
typedef typename MeshType::VertexType ElemType;
typedef typename std::vector<ElemType*>::iterator VertexIterator;
PlanarRegionVertex(): RegionBase<MeshType,typename MeshType::VertexType>(){};
// planes characterizing the region
vcg::Plane3<ScalarType,true> p;
CoordType center;
vcg::Color4b color;
ScalarType PlaneFittingError(std::vector<CoordType> & samples,vcg::Plane3<float> p);
// evaluate the gain if the triangle is added to the region
ScalarType Evaluate( ElemType & f);
// refit the plane
void Fit(std::vector<ElemType*>& faces, vcg::Plane3f & plane, float & err);
void Refit();
void Init(ElemType * );
};
template<class MeshType>
typename PlanarRegionVertex<MeshType>::ScalarType
PlanarRegionVertex< MeshType>:: PlaneFittingError( std::vector<CoordType> & samples,vcg::Plane3<float> p){
typename PlanarRegionVertex<MeshType>::ScalarType err =0.0;
typename std::vector< CoordType>::iterator si;
for(si = samples.begin(); si != samples.end(); ++si)
err += fabs((*si)*p.Direction()-p.Offset());
return err/samples.size();
}
// evaluate the gain if the vertex is added to the region
template<class MeshType>
typename PlanarRegionVertex<MeshType>::ScalarType PlanarRegionVertex< MeshType>::Evaluate( ElemType & e){
//return fabs(vcg::Distance(p,e.P()));
return (1-fabs(e.N()*p.Direction()));
}
// refit the planes
template<class MeshType>
void PlanarRegionVertex< MeshType>::Fit(std::vector<ElemType*>& verts, vcg::Plane3f & plane, float & err){
VertexIterator vi;
std::vector<CoordType> samples;
if(verts.size()<3){
samples.push_back( (*verts.begin())->P());
for(int i = 0; i < (*verts.begin())->SizeNeigh(); ++i)
samples.push_back((*verts.begin())->Neigh(i)->P());
vcg::PlaneFittingPoints<typename MeshType::ScalarType >(samples,plane);
}else
if(verts.size()==3){
plane.Init(verts[0]->P(),verts[1]->P(),verts[2]->P());
}else{
for( vi = verts.begin(); vi != verts.end(); ++vi)
samples.push_back( (*(*vi)).P());
vcg::PlaneFittingPoints<typename MeshType::ScalarType >(samples,plane);
}
err = PlaneFittingError(samples,plane);
}
template<class MeshType>
void PlanarRegionVertex< MeshType>::Refit(){
Fit(elems,p,this->approx_err);
}
template<class MeshType>
void PlanarRegionVertex< MeshType>::Init(ElemType * b){
std::vector<CoordType> samples;
samples.push_back(b->P());
for(int i = 0; i < b->SizeNeigh(); ++i)
samples.push_back(b->Neigh(i)->P());
vcg::PlaneFittingPoints<typename MeshType::ScalarType >(samples,this->p);
}
/**
Planar region made of vertices
*/
template<class MeshType>
struct PlanarRegionFace: public RegionBase<MeshType,typename MeshType::FaceType> {
public:
typedef PlanarRegionFace RegionType;
typedef typename MeshType::FaceType ElemType;
typedef typename std::vector<ElemType*>::iterator FaceIterator;
PlanarRegionFace(): RegionBase<MeshType,typename MeshType::FaceType>(){};
// planes characterizing the region
vcg::Plane3<ScalarType,true> p;
CoordType center;
vcg::Color4b color;
ScalarType PlaneFittingError(std::vector<CoordType> & samples,vcg::Plane3<float> p);
// evaluate the gain if the triangle is added to the region
ScalarType Evaluate( ElemType & f);
// refit the plane
void Fit(std::vector<ElemType*>& faces, vcg::Plane3f & plane, float & err);
void Refit();
void Init(ElemType * );
};
template<class MeshType>
typename PlanarRegionFace<MeshType>::ScalarType
PlanarRegionFace< MeshType>:: PlaneFittingError( std::vector<CoordType> & samples,vcg::Plane3<float> p){
typename PlanarRegionFace<MeshType>::ScalarType err =0.0;
typename std::vector< CoordType>::iterator si;
for(si = samples.begin(); si != samples.end(); ++si)
err += fabs((*si)*p.Direction()-p.Offset());
return err/samples.size();
}
// evaluate the gain if the triangle is added to the region
template<class MeshType>
typename PlanarRegionFace<MeshType>::ScalarType PlanarRegionFace< MeshType>::Evaluate( ElemType & f){
// return (f.N()-p.Direction()).SquaredNorm()*f.Q()*0.5;
return (vcg::NormalizedNormal(f)-p.Direction()).SquaredNorm()*vcg::DoubleArea(f)*0.5;
// return vcg::Distance(vcg::Barycenter(f),p)*vcg::DoubleArea(f)*0.5;
}
// refit the planes
template<class MeshType>
void PlanarRegionFace< MeshType>::Fit(std::vector<ElemType*>& faces, vcg::Plane3f & plane, float & err){
FaceIterator fi;
center = CoordType(0.0,0.0,0.0);
std::vector<CoordType> samples;
if(faces.size()<3){
samples.push_back((*faces.begin())->V(0)->P());
samples.push_back((*faces.begin())->V(1)->P());
samples.push_back((*faces.begin())->V(2)->P());
center+=vcg::Barycenter(*(*faces.begin()))*3;
}else
for( fi = faces.begin(); fi != faces.end(); ++fi)
{
//AddSamples(samples,*(*fi));
samples.push_back(vcg::Barycenter(*(*fi)));
center+=vcg::Barycenter(*(*fi));
}
center*=1/(typename MeshType::ScalarType)(samples.size());
if(samples.size()==3){
plane.SetDirection(vcg::Normal(vcg::Triangle3<typename PlanarRegionFace<MeshType>::ScalarType >(samples[0],samples[1],samples[2])));
typename MeshType::ScalarType off=samples[0]*plane.Direction();
plane.SetOffset(off);
}else
{ vcg::PlaneFittingPoints<typename MeshType::ScalarType >(samples,plane);
}
err = PlaneFittingError(samples,plane);
}
template<class MeshType>
void PlanarRegionFace< MeshType>::Refit(){
Fit(elems,p,this->approx_err);
}
template<class MeshType>
void PlanarRegionFace< MeshType>::Init(ElemType * b){
p.Init(vcg::Barycenter(*b),b->N());
}
} // namespace tri
}// namespace vcg
#endif