diff --git a/vcg/complex/algorithms/variational_shape_approximation.h b/vcg/complex/algorithms/variational_shape_approximation.h new file mode 100644 index 00000000..d6420d56 --- /dev/null +++ b/vcg/complex/algorithms/variational_shape_approximation.h @@ -0,0 +1,921 @@ +#ifndef VCGLIB_REGION_GROWING_VERTEX_ +#define VCGLIB_REGION_GROWING_VERTEX_ + +#include +#include +#include + +#include +#include +#include +#include +//#include + +namespace vcg{ + namespace tri{ + +/** +Element and associated error +*/ +template +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 +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::iterator AdjIterator; + typedef ElemError ElemError; + typedef typename std::vector::iterator ElemIte; + + // elem belonging to the region + std::vector elems; + + // adjacent regions + std::list 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 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 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 ::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& 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 ::iterator TriRegIterator; + typedef typename std::list* >::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_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::max(); + return fabs(samples[std::min(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::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 samples; + std::deque boxes; + int ns; + + }; + + + RegionGrower():lastAdded(NULL),ave_error(-1),ave_var(-1),n_steps(0){} + + std::list regions; + std::vector toTunnel; + std::vector 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 elemheap; + std::vector 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::iterator & ri){std::remove(ri);} + + + void PushHeap(std::vector & candi, RegionType & r){ + typename std::vector::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::iterator ri; + std::vector 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::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 & 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::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 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::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 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 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::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 { + typedef typename RegionType RegionType; + typedef typename RegionType::MeshType MeshType; + + typedef typename RegionType::ElemType ElemType; + typedef typename RegionType::ScalarType ScalarType; + typedef typename std::list ::iterator TriRegIterator; + + typename MeshType:: template PerVertexAttributeHandle attr_r; + typename MeshType:: template PerVertexAttributeHandle 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::Init(mesh, n_seeds,max_regions,max_err); + + attr_r = vcg::tri::Allocator::template GetPerVertexAttribute (*m,"r"); + if(!vcg::tri::Allocator::IsValidHandle(*m,attr_r)) + attr_r = vcg::tri::Allocator::template AddPerVertexAttribute (*m,"r"); + + attr_r_old = vcg::tri::Allocator::template GetPerVertexAttribute (*m,"r_old"); + if(!vcg::tri::Allocator::IsValidHandle(*m,attr_r_old)) + attr_r_old = vcg::tri::Allocator::template AddPerVertexAttribute (*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 & submeshes){ + + MeshType::PerVertexAttributeHandle sel = Allocator::AddPerVertexAttribute(*m,"sel"); + assert(Allocator::IsValidHandle(*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::iterator vi = (*ti).elems.begin();vi!=(*ti).elems.end(); ++vi) + {(*vi)->SetS(); } + + Append::Mesh( *submeshes.back(),*m,true); + + for(std::vector::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(0,std::min((rand()/float(RAND_MAX))*m->vert.size()-1,m->vert.size()-1)); + return &m->vert[id]; + } + + }; + + + +template < class RegionType> +struct RegionGrowerFace: public RegionGrower { + typedef typename RegionType RegionType; + typedef typename RegionType::MeshType MeshType; + + typedef typename RegionType::ElemType ElemType; + typedef typename RegionType::ScalarType ScalarType; + typedef typename std::list ::iterator TriRegIterator; + + typename MeshType:: template PerFaceAttributeHandle attr_r; + typename MeshType:: template PerFaceAttributeHandle 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::Init(mesh, n_seeds,max_regions,max_err); + + attr_r = vcg::tri::Allocator::template GetPerFaceAttribute (*m,"r"); + if(!vcg::tri::Allocator::IsValidHandle(*m,attr_r)) + attr_r = vcg::tri::Allocator::template AddPerFaceAttribute (*m,"r"); + + attr_r_old = vcg::tri::Allocator::template GetPerFaceAttribute (*m,"r_old"); + if(!vcg::tri::Allocator::IsValidHandle(*m,attr_r_old)) + attr_r_old = vcg::tri::Allocator::template AddPerFaceAttribute (*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 & submeshes){ + MeshType::PerFaceAttributeHandle sel = Allocator::AddPerFaceAttribute(*m,"sel"); + assert(Allocator::IsValidHandle(*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::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::Mesh( *submeshes.back(),*m,true); + + for(std::vector::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(0,std::min((rand()/float(RAND_MAX))*m->face.size()-1,m->face.size()-1)); + return &m->face[id]; + } + + }; + + +/** + Planar region made of vertices +*/ +template +struct PlanarRegionVertex: public RegionBase { +public: + + typedef PlanarRegionVertex RegionType; + typedef typename MeshType::VertexType ElemType; + typedef typename std::vector::iterator VertexIterator; + + PlanarRegionVertex(): RegionBase(){}; + + + // planes characterizing the region + vcg::Plane3 p; + + CoordType center; + vcg::Color4b color; + + ScalarType PlaneFittingError(std::vector & samples,vcg::Plane3 p); + + // evaluate the gain if the triangle is added to the region + ScalarType Evaluate( ElemType & f); + + // refit the plane + void Fit(std::vector& faces, vcg::Plane3f & plane, float & err); + void Refit(); + void Init(ElemType * ); +}; + + + template +typename PlanarRegionVertex::ScalarType + PlanarRegionVertex< MeshType>:: PlaneFittingError( std::vector & samples,vcg::Plane3 p){ + + typename PlanarRegionVertex::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 +typename PlanarRegionVertex::ScalarType PlanarRegionVertex< MeshType>::Evaluate( ElemType & e){ + //return fabs(vcg::Distance(p,e.P())); + return (1-fabs(e.N()*p.Direction())); +} + + // refit the planes +template +void PlanarRegionVertex< MeshType>::Fit(std::vector& verts, vcg::Plane3f & plane, float & err){ + VertexIterator vi; + std::vector 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(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(samples,plane); + } + + err = PlaneFittingError(samples,plane); +} + +template +void PlanarRegionVertex< MeshType>::Refit(){ + + Fit(elems,p,this->approx_err); +} + +template +void PlanarRegionVertex< MeshType>::Init(ElemType * b){ + std::vector samples; + samples.push_back(b->P()); + for(int i = 0; i < b->SizeNeigh(); ++i) + samples.push_back(b->Neigh(i)->P()); + vcg::PlaneFittingPoints(samples,this->p); +} + + +/** + Planar region made of vertices +*/ +template +struct PlanarRegionFace: public RegionBase { +public: + + typedef PlanarRegionFace RegionType; + typedef typename MeshType::FaceType ElemType; + typedef typename std::vector::iterator FaceIterator; + + PlanarRegionFace(): RegionBase(){}; + + + // planes characterizing the region + vcg::Plane3 p; + + CoordType center; + vcg::Color4b color; + + ScalarType PlaneFittingError(std::vector & samples,vcg::Plane3 p); + + // evaluate the gain if the triangle is added to the region + ScalarType Evaluate( ElemType & f); + + // refit the plane + void Fit(std::vector& faces, vcg::Plane3f & plane, float & err); + void Refit(); + void Init(ElemType * ); +}; + + +template +typename PlanarRegionFace::ScalarType + PlanarRegionFace< MeshType>:: PlaneFittingError( std::vector & samples,vcg::Plane3 p){ + + typename PlanarRegionFace::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 +typename PlanarRegionFace::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 +void PlanarRegionFace< MeshType>::Fit(std::vector& faces, vcg::Plane3f & plane, float & err){ + FaceIterator fi; + center = CoordType(0.0,0.0,0.0); + std::vector 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::ScalarType >(samples[0],samples[1],samples[2]))); + typename MeshType::ScalarType off=samples[0]*plane.Direction(); + plane.SetOffset(off); + }else + { vcg::PlaneFittingPoints(samples,plane); + } + + err = PlaneFittingError(samples,plane); + } + + template + void PlanarRegionFace< MeshType>::Refit(){ + Fit(elems,p,this->approx_err); + } + + template + void PlanarRegionFace< MeshType>::Init(ElemType * b){ + p.Init(vcg::Barycenter(*b),b->N()); + } + + + } // namespace tri +}// namespace vcg + + +#endif