diff --git a/vcg/complex/algorithms/autoalign_4pcs.h b/vcg/complex/algorithms/autoalign_4pcs.h index a720983b..5da1907a 100644 --- a/vcg/complex/algorithms/autoalign_4pcs.h +++ b/vcg/complex/algorithms/autoalign_4pcs.h @@ -31,9 +31,9 @@ ps: the name of the variables are out of vcg standard but like the one used in the paper pseudocode. */ +#include #include #include -#include #include // note: temporary (callback.h should be moved inside vcg) @@ -150,191 +150,191 @@ private: void GetBBox(vcg::Box3 & b){b.Add(pos);} }; - GridType *ugrid; // griglia - vcg::GridStaticPtr ugridQ; - vcg::GridStaticPtr ugridP; + GridType *ugrid; // griglia + vcg::GridStaticPtr ugridQ; + vcg::GridStaticPtr ugridP; - bool SelectCoplanarBase(); // on P - bool FindCongruent() ; // of base B, on Q, with approximation delta - void ComputeR1R2(ScalarType d1,ScalarType d2); + bool SelectCoplanarBase(); // on P + bool FindCongruent() ; // of base B, on Q, with approximation delta + void ComputeR1R2(ScalarType d1,ScalarType d2); - bool IsTransfCongruent(FourPoints fp,vcg::Matrix44 & mat, float & trerr); - int EvaluateSample(Candidate & fp, CoordType & tp, CoordType & np, const float & angle); - void EvaluateAlignment(Candidate & fp); - void TestAlignment(Candidate & fp); + bool IsTransfCongruent(FourPoints fp,vcg::Matrix44 & mat, float & trerr); + int EvaluateSample(Candidate & fp, CoordType & tp, CoordType & np, const float & angle); + void EvaluateAlignment(Candidate & fp); + void TestAlignment(Candidate & fp); - /* debug tools */ + /* debug tools */ public: - std::vector allTr;// tutte le trasformazioni provate - FILE * db; - char namemesh1[255],namemesh2[255]; - int n_base; - void InitDebug(const char * name1, const char * name2){ - db = fopen("debugPCS.txt","w"); - sprintf(&namemesh1[0],"%s",name1); - sprintf(&namemesh2[0],"%s",name2); - n_base = 0; - } + std::vector allTr;// tutte le trasformazioni provate + FILE * db; + char namemesh1[255],namemesh2[255]; + int n_base; + void InitDebug(const char * name1, const char * name2){ + db = fopen("debugPCS.txt","w"); + sprintf(&namemesh1[0],"%s",name1); + sprintf(&namemesh2[0],"%s",name2); + n_base = 0; + } - void FinishDebug(){ - fclose(db); - } - //void SaveALN(char * name,vcg::Matrix44f mat ){ - // FILE * o = fopen(name,"w"); - // fprintf(o,"2\n%s\n#\n",namemesh1); - // for(int i = 0 ; i < 4; ++i) - // fprintf(o,"%f %f %f %f\n",mat[i][0],mat[i][1],mat[i][2],mat[i][3]); - // fprintf(o,"%s\n#\n",namemesh2); - // fprintf(o,"1.0 0.0 0.0 0.0 \n"); - // fprintf(o,"0.0 1.0 0.0 0.0 \n"); - // fprintf(o,"0.0 0.0 1.0 0.0 \n"); - // fprintf(o,"0.0 0.0 0.0 1.0 \n"); + void FinishDebug(){ + fclose(db); + } + //void SaveALN(char * name,vcg::Matrix44f mat ){ + // FILE * o = fopen(name,"w"); + // fprintf(o,"2\n%s\n#\n",namemesh1); + // for(int i = 0 ; i < 4; ++i) + // fprintf(o,"%f %f %f %f\n",mat[i][0],mat[i][1],mat[i][2],mat[i][3]); + // fprintf(o,"%s\n#\n",namemesh2); + // fprintf(o,"1.0 0.0 0.0 0.0 \n"); + // fprintf(o,"0.0 1.0 0.0 0.0 \n"); + // fprintf(o,"0.0 0.0 1.0 0.0 \n"); + // fprintf(o,"0.0 0.0 0.0 1.0 \n"); - // fclose(o); - //} + // fclose(o); + //} }; template void FourPCS:: Init(MeshType &_P,MeshType &_Q) { - P = &_P;Q=&_Q; - ugridQ.Set(Q->vert.begin(),Q->vert.end()); - ugridP.Set(P->vert.begin(),P->vert.end()); + P = &_P;Q=&_Q; + ugridQ.Set(Q->vert.begin(),Q->vert.end()); + ugridP.Set(P->vert.begin(),P->vert.end()); - float ratio = 800 / (float) Q->vert.size(); - for(int vi = 0; vi < Q->vert.size(); ++vi) - if(rand()/(float) RAND_MAX < ratio) - mapsub.push_back(vi); + float ratio = 800 / (float) Q->vert.size(); + for(int vi = 0; vi < Q->vert.size(); ++vi) + if(rand()/(float) RAND_MAX < ratio) + mapsub.push_back(vi); - for(int vi = 0; vi < P->vert.size(); ++vi) - if(rand()/(float) RAND_MAX < ratio) - subsetP.push_back(&P->vert[vi]); + for(int vi = 0; vi < P->vert.size(); ++vi) + if(rand()/(float) RAND_MAX < ratio) + subsetP.push_back(&P->vert[vi]); - // estimate neigh distance - float avD = 0.0; - for(int i = 0 ; i < 100; ++i){ - int ri = rand()/(float) RAND_MAX * Q->vert.size() -1; - std::vector< CoordType > samples,d_samples; - std::vector dists; - std::vector ress; - vcg::tri::GetKClosestVertex< - MeshType, - vcg::GridStaticPtr, - std::vector, - std::vector, - std::vector< CoordType > >(*Q,ugridQ,2,Q->vert[ri].cP(),Q->bbox.Diag(), ress,dists, samples); - assert(ress.size() == 2); - avD+=dists[1]; - } - avD /=100; // average vertex-vertex distance - avD /= sqrt(ratio); // take into account the ratio + // estimate neigh distance + float avD = 0.0; + for(int i = 0 ; i < 100; ++i){ + int ri = rand()/(float) RAND_MAX * Q->vert.size() -1; + std::vector< CoordType > samples,d_samples; + std::vector dists; + std::vector ress; + vcg::tri::GetKClosestVertex< + MeshType, + vcg::GridStaticPtr, + std::vector, + std::vector, + std::vector< CoordType > >(*Q,ugridQ,2,Q->vert[ri].cP(),Q->bbox.Diag(), ress,dists, samples); + assert(ress.size() == 2); + avD+=dists[1]; + } + avD /=100; // average vertex-vertex distance + avD /= sqrt(ratio); // take into account the ratio - par.delta = avD * par.delta; - side = P->bbox.Dim()[P->bbox.MaxDim()]*par.f; //rough implementation + par.delta = avD * par.delta; + side = P->bbox.Dim()[P->bbox.MaxDim()]*par.f; //rough implementation - } + } template bool FourPCS::SelectCoplanarBase(){ - vcg::tri::UpdateBounding::Box(*P); + vcg::tri::UpdateBounding::Box(*P); - // choose the inter point distance - ScalarType dtol = side*0.1; //rough implementation + // choose the inter point distance + ScalarType dtol = side*0.1; //rough implementation - //choose the first two points - int i = 0,ch; + //choose the first two points + int i = 0,ch; - // first point random - ch = (rand()/(float)RAND_MAX)*(P->vert.size()-2); - B[0] = P->vert[ch].P(); + // first point random + ch = (rand()/(float)RAND_MAX)*(P->vert.size()-2); + B[0] = P->vert[ch].P(); //printf("B[0] %d\n",ch); - // second a point at distance d+-dtol - for(i = 0; i < P->vert.size(); ++i){ - ScalarType dd = (P->vert[i].P() - B[0]).Norm(); - if( ( dd < side + dtol) && (dd > side - dtol)){ - B[1] = P->vert[i].P(); + // second a point at distance d+-dtol + for(i = 0; i < P->vert.size(); ++i){ + ScalarType dd = (P->vert[i].P() - B[0]).Norm(); + if( ( dd < side + dtol) && (dd > side - dtol)){ + B[1] = P->vert[i].P(); //printf("B[1] %d\n",i); - break; - } - } - if(i == P->vert.size()) - return false; + break; + } + } + if(i == P->vert.size()) + return false; - // third point at distance d from B[1] and forming a right angle - int best = -1; ScalarType bestv=std::numeric_limits::max(); - for(i = 0; i < P->vert.size(); ++i){ - int id = rand()/(float)RAND_MAX * (P->vert.size()-1); - ScalarType dd = (P->vert[id].P() - B[1]).Norm(); - if( ( dd < side + dtol) && (dd > side - dtol)){ - ScalarType angle = fabs( ( P->vert[id].P()-B[1]).normalized().dot((B[1]-B[0]).normalized())); - if( angle < bestv){ - bestv = angle; - best = id; - } - } - } - if(best == -1) - return false; - B[2] = P->vert[best].P(); + // third point at distance d from B[1] and forming a right angle + int best = -1; ScalarType bestv=std::numeric_limits::max(); + for(i = 0; i < P->vert.size(); ++i){ + int id = rand()/(float)RAND_MAX * (P->vert.size()-1); + ScalarType dd = (P->vert[id].P() - B[1]).Norm(); + if( ( dd < side + dtol) && (dd > side - dtol)){ + ScalarType angle = fabs( ( P->vert[id].P()-B[1]).normalized().dot((B[1]-B[0]).normalized())); + if( angle < bestv){ + bestv = angle; + best = id; + } + } + } + if(best == -1) + return false; + B[2] = P->vert[best].P(); //printf("B[2] %d\n",best); - CoordType n = ((B[0]-B[1]).normalized() ^ (B[2]-B[1]).normalized()).normalized(); - CoordType B4 = B[1] + (B[0]-B[1]) + (B[2]-B[1]); - VertexType * v =0; - ScalarType radius = dtol*4.0; + CoordType n = ((B[0]-B[1]).normalized() ^ (B[2]-B[1]).normalized()).normalized(); + CoordType B4 = B[1] + (B[0]-B[1]) + (B[2]-B[1]); + VertexType * v =0; + ScalarType radius = dtol*4.0; - std::vector closests; - std::vector distances; - std::vector points; + std::vector closests; + std::vector distances; + std::vector points; - vcg::tri::GetInSphereVertex< - MeshType, - vcg::GridStaticPtr, - std::vector, - std::vector, - std::vector - >(*P,ugridP,B4,radius,closests,distances,points); + vcg::tri::GetInSphereVertex< + MeshType, + vcg::GridStaticPtr, + std::vector, + std::vector, + std::vector + >(*P,ugridP,B4,radius,closests,distances,points); - if(closests.empty()) - return false; - best = -1; bestv=std::numeric_limits::max(); - for(i = 0; i P() - B[1]).normalized().dot(n)); - if( angle < bestv){ - bestv = angle; - best = i; - } - } - B[3] = closests[best]->P(); + if(closests.empty()) + return false; + best = -1; bestv=std::numeric_limits::max(); + for(i = 0; i P() - B[1]).normalized().dot(n)); + if( angle < bestv){ + bestv = angle; + best = i; + } + } + B[3] = closests[best]->P(); //printf("B[3] %d\n", (typename MeshType::VertexType*)closests[best] - &(*P->vert.begin())); - // compute r1 and r2 - CoordType x; - std::swap(B[1],B[2]); - IntersectionLineLine(B[0],B[1],B[2],B[3],x); + // compute r1 and r2 + CoordType x; + std::swap(B[1],B[2]); + IntersectionLineLine(B[0],B[1],B[2],B[3],x); - r1 = (x - B[0]).dot(B[1]-B[0]) / (B[1]-B[0]).SquaredNorm(); - r2 = (x - B[2]).dot(B[3]-B[2]) / (B[3]-B[2]).SquaredNorm(); + r1 = (x - B[0]).dot(B[1]-B[0]) / (B[1]-B[0]).SquaredNorm(); + r2 = (x - B[2]).dot(B[3]-B[2]) / (B[3]-B[2]).SquaredNorm(); - if( ((B[0]+(B[1]-B[0])*r1)-(B[2]+(B[3]-B[2])*r2)).Norm() > par.delta ) - return false; + if( ((B[0]+(B[1]-B[0])*r1)-(B[2]+(B[3]-B[2])*r2)).Norm() > par.delta ) + return false; - radius =side*0.5; - std::vector< CoordType > samples,d_samples; - std::vector dists; + radius =side*0.5; + std::vector< CoordType > samples,d_samples; + std::vector dists; - for(int i = 0 ; i< 4; ++i){ - vcg::tri::GetKClosestVertex< - MeshType, - vcg::GridStaticPtr, - std::vector, - std::vector, - std::vector< CoordType > >(*P,ugridP, par.feetsize ,B[i],radius, ExtB[i],dists, samples); - } + for(int i = 0 ; i< 4; ++i){ + vcg::tri::GetKClosestVertex< + MeshType, + vcg::GridStaticPtr, + std::vector, + std::vector, + std::vector< CoordType > >(*P,ugridP, par.feetsize ,B[i],radius, ExtB[i],dists, samples); + } //for(int i = 0 ; i< 4; ++i) // printf("%d ",ExtB[i].size()); @@ -372,144 +372,144 @@ bool FourPCS::IsTransfCongruent(FourPoints fp, vcg::Matrix44 void FourPCS::ComputeR1R2(ScalarType d1,ScalarType d2){ - int vi,vj; - R1.clear(); - //R2.clear(); - int start = clock(); - for(vi = 0; vi < mapsub.size(); ++vi) for(vj = vi; vj < mapsub.size(); ++vj){ - ScalarType d = ((Q->vert[mapsub[vi]]).P()-(Q->vert[mapsub[vj]]).P()).Norm(); - if( (d < d1+ side*0.5) && (d > d1-side*0.5)) - { - R1.push_back(Couple(mapsub[vi],mapsub[vj],d )); - R1.push_back(Couple(mapsub[vj],mapsub[vi],d)); - } - } - //for( vi = 0; vi < mapsub.size(); ++ vi ) for( vj = vi ; vj < mapsub.size(); ++ vj ){ - // ScalarType d = ((Q->vert[mapsub[vi]]).P()-(Q->vert[mapsub[vj]]).P()).Norm(); - // if( (d < d2+side*0.5) && (d > d2-side*0.5)) - // { - // R2.push_back(Couple(mapsub[vi],mapsub[vj],d)); - // R2.push_back(Couple(mapsub[vj],mapsub[vi],d)); - // } - //} + int vi,vj; + R1.clear(); + //R2.clear(); + int start = clock(); + for(vi = 0; vi < mapsub.size(); ++vi) for(vj = vi; vj < mapsub.size(); ++vj){ + ScalarType d = ((Q->vert[mapsub[vi]]).P()-(Q->vert[mapsub[vj]]).P()).Norm(); + if( (d < d1+ side*0.5) && (d > d1-side*0.5)) + { + R1.push_back(Couple(mapsub[vi],mapsub[vj],d )); + R1.push_back(Couple(mapsub[vj],mapsub[vi],d)); + } + } + //for( vi = 0; vi < mapsub.size(); ++ vi ) for( vj = vi ; vj < mapsub.size(); ++ vj ){ + // ScalarType d = ((Q->vert[mapsub[vi]]).P()-(Q->vert[mapsub[vj]]).P()).Norm(); + // if( (d < d2+side*0.5) && (d > d2-side*0.5)) + // { + // R2.push_back(Couple(mapsub[vi],mapsub[vj],d)); + // R2.push_back(Couple(mapsub[vj],mapsub[vi],d)); + // } + //} - std::sort(R1.begin(),R1.end()); + std::sort(R1.begin(),R1.end()); // std::sort(R2.begin(),R2.end()); } template bool FourPCS::FindCongruent() { // of base B, on Q, with approximation delta - bool done = false; - std::vector R2inv; - int n_closests = 0, n_congr = 0; - int ac =0 ,acf = 0,tr = 0,trf =0; - ScalarType d1,d2; - d1 = (B[1]-B[0]).Norm(); - d2 = (B[3]-B[2]).Norm(); + bool done = false; + std::vector R2inv; + int n_closests = 0, n_congr = 0; + int ac =0 ,acf = 0,tr = 0,trf =0; + ScalarType d1,d2; + d1 = (B[1]-B[0]).Norm(); + d2 = (B[3]-B[2]).Norm(); - int start = clock(); - //int vi,vj; + int start = clock(); + //int vi,vj; - typename PMesh::VertexIterator vii; - typename std::vector::iterator bR1,eR1,bR2,eR2,ite,cite; - bR1 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d1-par.delta*2.0)); - eR1 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d1+par.delta*2.0)); - bR2 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d2-par.delta*2.0)); - eR2 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d2+par.delta*2.0)); + typename PMesh::VertexIterator vii; + typename std::vector::iterator bR1,eR1,bR2,eR2,ite,cite; + bR1 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d1-par.delta*2.0)); + eR1 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d1+par.delta*2.0)); + bR2 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d2-par.delta*2.0)); + eR2 = std::lower_bound::iterator,Couple>(R1.begin(),R1.end(),Couple(d2+par.delta*2.0)); - // in [bR1,eR1) there are all the pairs ad a distance d1 +- par.delta - // in [bR1,eR1) there are all the pairs ad a distance d2 +- par.delta + // in [bR1,eR1) there are all the pairs ad a distance d1 +- par.delta + // in [bR1,eR1) there are all the pairs ad a distance d2 +- par.delta - if(bR1 == R1.end()) return false;// if there are no such pairs return - if(bR2 == R1.end()) return false; // if there are no such pairs return + if(bR1 == R1.end()) return false;// if there are no such pairs return + if(bR2 == R1.end()) return false; // if there are no such pairs return - // put [bR1,eR1) in a mesh to have the search operator for free (lazy me) - Invr.Clear(); - int i = &(*bR1)-&(*R1.begin()); - for(ite = bR1; ite != eR1;++ite){ - vii = vcg::tri::Allocator::AddVertices(Invr,1); - (*vii).P() = Q->vert[R1[i][0]].P() + (Q->vert[R1[i][1]].P()-Q->vert[R1[i][0]].P()) * r1; - ++i; - } - if(Invr.vert.empty() ) return false; + // put [bR1,eR1) in a mesh to have the search operator for free (lazy me) + Invr.Clear(); + int i = &(*bR1)-&(*R1.begin()); + for(ite = bR1; ite != eR1;++ite){ + vii = vcg::tri::Allocator::AddVertices(Invr,1); + (*vii).P() = Q->vert[R1[i][0]].P() + (Q->vert[R1[i][1]].P()-Q->vert[R1[i][0]].P()) * r1; + ++i; + } + if(Invr.vert.empty() ) return false; // index remaps a vertex of Invr to its corresponding point in R1 typename PMesh::template PerVertexAttributeHandle id = vcg::tri::Allocator::template AddPerVertexAttribute(Invr,std::string("index")); i = &(*bR1)-&(*R1.begin()); for(vii = Invr.vert.begin(); vii != Invr.vert.end();++vii,++i) id[vii] = i; - vcg::tri::UpdateBounding::Box(Invr); - // printf("Invr size %d\n",Invr.vn); + vcg::tri::UpdateBounding::Box(Invr); + // printf("Invr size %d\n",Invr.vn); - ugrid = new GridType(); - ugrid->Set(Invr.vert.begin(),Invr.vert.end()); + ugrid = new GridType(); + ugrid->Set(Invr.vert.begin(),Invr.vert.end()); - i = &(*bR2)-&(*R1.begin()); - // R2inv contains all the points generated by the couples in R2 (with the reference to remap into R2) - for(ite = bR2; ite != eR2;++ite){ - R2inv.push_back( EPoint( Q->vert[R1[i][0]].P() + (Q->vert[R1[i][1]].P()-Q->vert[R1[i][0]].P()) * r2,i)); - ++i; - } + i = &(*bR2)-&(*R1.begin()); + // R2inv contains all the points generated by the couples in R2 (with the reference to remap into R2) + for(ite = bR2; ite != eR2;++ite){ + R2inv.push_back( EPoint( Q->vert[R1[i][0]].P() + (Q->vert[R1[i][1]].P()-Q->vert[R1[i][0]].P()) * r2,i)); + ++i; + } - n_closests = 0; n_congr = 0; ac =0 ; acf = 0; tr = 0; trf = 0; - printf("R2Inv.size = %d \n",R2inv.size()); - for(uint i = 0 ; i < R2inv.size() ; ++i){ + n_closests = 0; n_congr = 0; ac =0 ; acf = 0; tr = 0; trf = 0; + printf("R2Inv.size = %d \n",R2inv.size()); + for(uint i = 0 ; i < R2inv.size() ; ++i){ - std::vector closests; + std::vector closests; - // for each point in R2inv get all the points in R1 closer than par.delta - vcg::Matrix44 mat; - vcg::Box3f bb; - bb.Add(R2inv[i].pos+vcg::Point3f(par.delta * 0.1,par.delta * 0.1 , par.delta * 0.1 )); - bb.Add(R2inv[i].pos-vcg::Point3f(par.delta * 0.1,par.delta* 0.1 , par.delta* 0.1)); + // for each point in R2inv get all the points in R1 closer than par.delta + vcg::Matrix44 mat; + vcg::Box3f bb; + bb.Add(R2inv[i].pos+vcg::Point3f(par.delta * 0.1,par.delta * 0.1 , par.delta * 0.1 )); + bb.Add(R2inv[i].pos-vcg::Point3f(par.delta * 0.1,par.delta* 0.1 , par.delta* 0.1)); - vcg::tri::GetInBoxVertex > - (Invr,*ugrid,bb,closests); + vcg::tri::GetInBoxVertex > + (Invr,*ugrid,bb,closests); - n_closests+=closests.size(); - for(uint ip = 0; ip < closests.size(); ++ip){ - FourPoints p; - p[0] = Q->vert[R1[id[closests[ip]]][0]].P(); - p[1] = Q->vert[R1[id[closests[ip]]][1]].P(); - p[2] = Q->vert[R1[ R2inv[i].pi][0]].P(); - p[3] = Q->vert[R1[ R2inv[i].pi][1]].P(); + n_closests+=closests.size(); + for(uint ip = 0; ip < closests.size(); ++ip){ + FourPoints p; + p[0] = Q->vert[R1[id[closests[ip]]][0]].P(); + p[1] = Q->vert[R1[id[closests[ip]]][1]].P(); + p[2] = Q->vert[R1[ R2inv[i].pi][0]].P(); + p[3] = Q->vert[R1[ R2inv[i].pi][1]].P(); - float trerr; - n_base++; - if(!IsTransfCongruent(p,mat,trerr)) { - trf++; - //char name[255]; - //sprintf(name,"faileTR_%d_%f.aln",n_base,trerr); - //fprintf(db,"TransCongruent %s\n", name); - //SaveALN(name, mat); - } - else{ - tr++; - n_congr++; - U.push_back(Candidate(p,mat)); - EvaluateAlignment(U.back()); - U.back().base = bases.size()-1; + float trerr; + n_base++; + if(!IsTransfCongruent(p,mat,trerr)) { + trf++; + //char name[255]; + //sprintf(name,"faileTR_%d_%f.aln",n_base,trerr); + //fprintf(db,"TransCongruent %s\n", name); + //SaveALN(name, mat); + } + else{ + tr++; + n_congr++; + U.push_back(Candidate(p,mat)); + EvaluateAlignment(U.back()); + U.back().base = bases.size()-1; - if( U.back().score > par.scoreFeet){ - TestAlignment(U.back()); - if(U.back().score > par.scoreAln) - { - done = true; break; - } - } - //char name[255]; - //sprintf(name,"passed_score_%5d_%d.aln",U.back().score,n_base); - //fprintf(db,"OK TransCongruent %s, score: %d \n", name,U.back().score); - //SaveALN(name, mat); - } - } - } + if( U.back().score > par.scoreFeet){ + TestAlignment(U.back()); + if(U.back().score > par.scoreAln) + { + done = true; break; + } + } + //char name[255]; + //sprintf(name,"passed_score_%5d_%d.aln",U.back().score,n_base); + //fprintf(db,"OK TransCongruent %s, score: %d \n", name,U.back().score); + //SaveALN(name, mat); + } + } + } - delete ugrid; - vcg::tri::Allocator::DeletePerVertexAttribute(Invr,id); - printf("n_closests %5d = (An %5d ) + ( Tr %5d ) + (OK) %5d\n",n_closests,acf,trf,n_congr); + delete ugrid; + vcg::tri::Allocator::DeletePerVertexAttribute(Invr,id); + printf("n_closests %5d = (An %5d ) + ( Tr %5d ) + (OK) %5d\n",n_closests,acf,trf,n_congr); - return done; + return done; // printf("done n_closests %d congr %d in %f s\n ",n_closests,n_congr,(clock()-start)/(float)CLOCKS_PER_SEC); // printf("angle:%d %d, trasf %d %d\n",ac,acf,tr,trf); } @@ -550,28 +550,28 @@ int FourPCS::EvaluateSample(Candidate & fp, CoordType & tp, CoordType template void FourPCS::EvaluateAlignment(Candidate & fp){ - int n_delta_close = 0; - for(int i = 0 ; i< 4; ++i) { - for(uint j = 0; j < ExtB[i].size();++j){ - CoordType np = ExtB[i][j]->cN();; - CoordType tp = ExtB[i][j]->P(); - n_delta_close+=EvaluateSample(fp,tp,np,0.9); - } - } - fp.score = n_delta_close; + int n_delta_close = 0; + for(int i = 0 ; i< 4; ++i) { + for(uint j = 0; j < ExtB[i].size();++j){ + CoordType np = ExtB[i][j]->cN();; + CoordType tp = ExtB[i][j]->P(); + n_delta_close+=EvaluateSample(fp,tp,np,0.9); + } + } + fp.score = n_delta_close; } template void FourPCS::TestAlignment(Candidate & fp){ - radius = par.delta; - int n_delta_close = 0; - for(uint j = 0; j < subsetP.size();++j){ - CoordType np = subsetP[j]->N(); - CoordType tp = subsetP[j]->P(); - n_delta_close+=EvaluateSample(fp,tp,np,0.6); - } - fp.score = n_delta_close; + radius = par.delta; + int n_delta_close = 0; + for(uint j = 0; j < subsetP.size();++j){ + CoordType np = subsetP[j]->N(); + CoordType tp = subsetP[j]->P(); + n_delta_close+=EvaluateSample(fp,tp,np,0.6); + } + fp.score = n_delta_close; } @@ -579,50 +579,50 @@ template bool FourPCS:: Align( int L, vcg::Matrix44f & result, vcg::CallBackPos * cb ){ // main loop - int bestv = 0; - bool found; - int n_tries = 0; - U.clear(); + int bestv = 0; + bool found; + int n_tries = 0; + U.clear(); - if(L==0) - { - L = (log(1.0-0.9999) / log(1.0-pow((float)par.f,3.f)))+1; - printf("using %d bases\n",L); - } + if(L==0) + { + L = (log(1.0-0.9999) / log(1.0-pow((float)par.f,3.f)))+1; + printf("using %d bases\n",L); + } - ComputeR1R2(side*1.4,side*1.4); + ComputeR1R2(side*1.4,side*1.4); - for(int t = 0; t < L; ++t ){ - do{ - n_tries = 0; - do{ - n_tries++; - found = SelectCoplanarBase(); - } - while(!found && (n_tries <50)); - if(!found) { - par.f*=0.98; - side = P->bbox.Dim()[P->bbox.MaxDim()]*par.f; //rough implementation - ComputeR1R2(side*1.4,side*1.4); - } - } while (!found && (par.f >0.1)); + for(int t = 0; t < L; ++t ){ + do{ + n_tries = 0; + do{ + n_tries++; + found = SelectCoplanarBase(); + } + while(!found && (n_tries <50)); + if(!found) { + par.f*=0.98; + side = P->bbox.Dim()[P->bbox.MaxDim()]*par.f; //rough implementation + ComputeR1R2(side*1.4,side*1.4); + } + } while (!found && (par.f >0.1)); - if(par.f <0.1) { - printf("FAILED"); - return false; - } - bases.push_back(B); - if(cb) cb(t*100/L,"trying bases"); - if(FindCongruent()) - break; - } + if(par.f <0.1) { + printf("FAILED"); + return false; + } + bases.push_back(B); + if(cb) cb(t*100/L,"trying bases"); + if(FindCongruent()) + break; + } - if(U.empty()) return false; + if(U.empty()) return false; - std::sort(U.begin(),U.end()); + std::sort(U.begin(),U.end()); - bestv = -std::numeric_limits::max(); - iwinner = 0; + bestv = -std::numeric_limits::max(); + iwinner = 0; for(int i = 0 ; i < U.size() ;++i) { @@ -635,15 +635,15 @@ FourPCS:: Align( int L, vcg::Matrix44f & result, vcg::CallBackPos * printf("Best score: %d \n", bestv); - winner = U[iwinner]; - result = winner.T; + winner = U[iwinner]; + result = winner.T; - // deallocations - Invr.Clear(); + // deallocations + Invr.Clear(); - return true; + return true; } - } // namespace tri + } // namespace tri } // namespace vcg #endif diff --git a/vcg/complex/algorithms/create/advancing_front.h b/vcg/complex/algorithms/create/advancing_front.h index f4a97c8b..cfbda931 100644 --- a/vcg/complex/algorithms/create/advancing_front.h +++ b/vcg/complex/algorithms/create/advancing_front.h @@ -3,7 +3,6 @@ #include #include -#include #include #include diff --git a/vcg/complex/algorithms/create/mc_trivial_walker.h b/vcg/complex/algorithms/create/mc_trivial_walker.h index 775ff0a9..ef19d331 100644 --- a/vcg/complex/algorithms/create/mc_trivial_walker.h +++ b/vcg/complex/algorithms/create/mc_trivial_walker.h @@ -22,7 +22,6 @@ ****************************************************************************/ #ifndef __VCG_TRIVIAL_WALKER #define __VCG_TRIVIAL_WALKER -#include namespace vcg { @@ -56,36 +55,36 @@ public: //else return numeric_limits::quiet_NaN( ); } - VOX_TYPE &V(const int &x,const int &y,const int &z) { - return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; - } + VOX_TYPE &V(const int &x,const int &y,const int &z) { + return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; + } - VOX_TYPE &V(const Point3i &pi) { - return Vol[ pi[0] + pi[1]*sz[0] + pi[2]*sz[0]*sz[1] ]; - } + VOX_TYPE &V(const Point3i &pi) { + return Vol[ pi[0] + pi[1]*sz[0] + pi[2]*sz[0]*sz[1] ]; + } - const VOX_TYPE &cV(const int &x,const int &y,const int &z) const { - return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; - } + const VOX_TYPE &cV(const int &x,const int &y,const int &z) const { + return Vol[x+y*sz[0]+z*sz[0]*sz[1]]; + } - typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis; + typedef enum { XAxis=0,YAxis=1,ZAxis=2} VolumeAxis; - template < class VertexPointerType, VolumeAxis AxisVal > - void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) - { - float f1 = V(p1).V()-thr; - float f2 = V(p2).V()-thr; - float u = (float) f1/(f1-f2); - if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X(); - else v->P().X() = (float) p1.X(); - if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y(); - else v->P().Y() = (float) p1.Y(); - if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z(); - else v->P().Z() = (float) p1.Z(); + template < class VertexPointerType, VolumeAxis AxisVal > + void GetIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) + { + float f1 = V(p1).V()-thr; + float f2 = V(p2).V()-thr; + float u = (float) f1/(f1-f2); + if(AxisVal==XAxis) v->P().X() = (float) p1.X()*(1-u) + u*p2.X(); + else v->P().X() = (float) p1.X(); + if(AxisVal==YAxis) v->P().Y() = (float) p1.Y()*(1-u) + u*p2.Y(); + else v->P().Y() = (float) p1.Y(); + if(AxisVal==ZAxis) v->P().Z() = (float) p1.Z()*(1-u) + u*p2.Z(); + else v->P().Z() = (float) p1.Z(); - if(VoxelType::HasNormal()) v->N() = V(p1).N()*(1-u) + V(p2).N()*u; - } + if(VoxelType::HasNormal()) v->N() = V(p1).N()*(1-u) + V(p2).N()*u; + } template < class VertexPointerType > void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointerType &v, const float thr) @@ -157,16 +156,16 @@ private: _bbox = Box3i(Point3i(0,0,0),volume.ISize()); _slice_dimension = _bbox.DimX()*_bbox.DimZ(); - _x_cs = new VertexIndex[ _slice_dimension ]; - _y_cs = new VertexIndex[ _slice_dimension ]; - _z_cs = new VertexIndex[ _slice_dimension ]; - _x_ns = new VertexIndex[ _slice_dimension ]; - _z_ns = new VertexIndex[ _slice_dimension ]; + _x_cs = new VertexIndex[ _slice_dimension ]; + _y_cs = new VertexIndex[ _slice_dimension ]; + _z_cs = new VertexIndex[ _slice_dimension ]; + _x_ns = new VertexIndex[ _slice_dimension ]; + _z_ns = new VertexIndex[ _slice_dimension ]; - }; + }; - ~TrivialWalker() - {_thr=0;} + ~TrivialWalker() + {_thr=0;} template void BuildMesh(MeshType &mesh, VolumeType &volume, EXTRACTOR_TYPE &extractor, const float threshold, vcg::CallBackPos * cb=0) @@ -178,148 +177,148 @@ private: _thr=threshold; vcg::Point3i p1, p2; - Begin(); - extractor.Initialize(); - for (int j=_bbox.min.Y(); j<(_bbox.max.Y()-1)-1; j+=1) - { + Begin(); + extractor.Initialize(); + for (int j=_bbox.min.Y(); j<(_bbox.max.Y()-1)-1; j+=1) + { - if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume"); + if(cb && ((j%10)==0) ) cb(j*_bbox.DimY()/100.0,"Marching volume"); - for (int i=_bbox.min.X(); i<(_bbox.max.X()-1)-1; i+=1) - { - for (int k=_bbox.min.Z(); k<(_bbox.max.Z()-1)-1; k+=1) - { - p1.X()=i; p1.Y()=j; p1.Z()=k; - p2.X()=i+1; p2.Y()=j+1; p2.Z()=k+1; - extractor.ProcessCell(p1, p2); - } - } - NextSlice(); - } - extractor.Finalize(); - _volume = NULL; - _mesh = NULL; - }; + for (int i=_bbox.min.X(); i<(_bbox.max.X()-1)-1; i+=1) + { + for (int k=_bbox.min.Z(); k<(_bbox.max.Z()-1)-1; k+=1) + { + p1.X()=i; p1.Y()=j; p1.Z()=k; + p2.X()=i+1; p2.Y()=j+1; p2.Z()=k+1; + extractor.ProcessCell(p1, p2); + } + } + NextSlice(); + } + extractor.Finalize(); + _volume = NULL; + _mesh = NULL; + }; - float V(int pi, int pj, int pk) - { - return _volume->Val(pi, pj, pk)-_thr; - } + float V(int pi, int pj, int pk) + { + return _volume->Val(pi, pj, pk)-_thr; + } - bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v) - { - int pos = p0.X()+p0.Z()*_bbox.max.X(); - int vidx; + bool Exist(const vcg::Point3i &p0, const vcg::Point3i &p1, VertexPointer &v) + { + int pos = p0.X()+p0.Z()*_bbox.max.X(); + int vidx; - if (p0.X()!=p1.X()) // punti allineati lungo l'asse X - vidx = (p0.Y()==_current_slice) ? _x_cs[pos] : _x_ns[pos]; - else if (p0.Y()!=p1.Y()) // punti allineati lungo l'asse Y - vidx = _y_cs[pos]; - else if (p0.Z()!=p1.Z()) // punti allineati lungo l'asse Z - vidx = (p0.Y()==_current_slice)? _z_cs[pos] : _z_ns[pos]; - else - assert(false); + if (p0.X()!=p1.X()) // punti allineati lungo l'asse X + vidx = (p0.Y()==_current_slice) ? _x_cs[pos] : _x_ns[pos]; + else if (p0.Y()!=p1.Y()) // punti allineati lungo l'asse Y + vidx = _y_cs[pos]; + else if (p0.Z()!=p1.Z()) // punti allineati lungo l'asse Z + vidx = (p0.Y()==_current_slice)? _z_cs[pos] : _z_ns[pos]; + else + assert(false); - v = (vidx!=-1)? &_mesh->vert[vidx] : NULL; - return v!=NULL; - } + v = (vidx!=-1)? &_mesh->vert[vidx] : NULL; + return v!=NULL; + } - void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) - { - int i = p1.X() - _bbox.min.X(); - int z = p1.Z() - _bbox.min.Z(); - VertexIndex index = i+z*_bbox.max.X(); - VertexIndex pos; - if (p1.Y()==_current_slice) - { - if ((pos=_x_cs[index])==-1) - { - _x_cs[index] = (VertexIndex) _mesh->vert.size(); - pos = _x_cs[index]; - Allocator::AddVertices( *_mesh, 1 ); - v = &_mesh->vert[pos]; - _volume->GetXIntercept(p1, p2, v, _thr); - return; - } - } - if (p1.Y()==_current_slice+1) - { - if ((pos=_x_ns[index])==-1) - { - _x_ns[index] = (VertexIndex) _mesh->vert.size(); - pos = _x_ns[index]; - Allocator::AddVertices( *_mesh, 1 ); - v = &_mesh->vert[pos]; - _volume->GetXIntercept(p1, p2, v,_thr); - return; - } - } - assert(pos >=0 && size_t(pos)< _mesh->vert.size()); - v = &_mesh->vert[pos]; - } - void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) - { - int i = p1.X() - _bbox.min.X(); - int z = p1.Z() - _bbox.min.Z(); - VertexIndex index = i+z*_bbox.max.X(); - VertexIndex pos; - if ((pos=_y_cs[index])==-1) - { - _y_cs[index] = (VertexIndex) _mesh->vert.size(); - pos = _y_cs[index]; - Allocator::AddVertices( *_mesh, 1); - v = &_mesh->vert[ pos ]; - _volume->GetYIntercept(p1, p2, v,_thr); - } - v = &_mesh->vert[pos]; - } - void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) - { - int i = p1.X() - _bbox.min.X(); - int z = p1.Z() - _bbox.min.Z(); - VertexIndex index = i+z*_bbox.max.X(); - VertexIndex pos; - if (p1.Y()==_current_slice) - { - if ((pos=_z_cs[index])==-1) - { - _z_cs[index] = (VertexIndex) _mesh->vert.size(); - pos = _z_cs[index]; - Allocator::AddVertices( *_mesh, 1 ); - v = &_mesh->vert[pos]; - _volume->GetZIntercept(p1, p2, v,_thr); - return; - } - } - if (p1.Y()==_current_slice+1) - { - if ((pos=_z_ns[index])==-1) - { - _z_ns[index] = (VertexIndex) _mesh->vert.size(); - pos = _z_ns[index]; - Allocator::AddVertices( *_mesh, 1 ); - v = &_mesh->vert[pos]; - _volume->GetZIntercept(p1, p2, v,_thr); - return; - } - } - v = &_mesh->vert[pos]; - } + void GetXIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + { + int i = p1.X() - _bbox.min.X(); + int z = p1.Z() - _bbox.min.Z(); + VertexIndex index = i+z*_bbox.max.X(); + VertexIndex pos; + if (p1.Y()==_current_slice) + { + if ((pos=_x_cs[index])==-1) + { + _x_cs[index] = (VertexIndex) _mesh->vert.size(); + pos = _x_cs[index]; + Allocator::AddVertices( *_mesh, 1 ); + v = &_mesh->vert[pos]; + _volume->GetXIntercept(p1, p2, v, _thr); + return; + } + } + if (p1.Y()==_current_slice+1) + { + if ((pos=_x_ns[index])==-1) + { + _x_ns[index] = (VertexIndex) _mesh->vert.size(); + pos = _x_ns[index]; + Allocator::AddVertices( *_mesh, 1 ); + v = &_mesh->vert[pos]; + _volume->GetXIntercept(p1, p2, v,_thr); + return; + } + } + assert(pos >=0 && size_t(pos)< _mesh->vert.size()); + v = &_mesh->vert[pos]; + } + void GetYIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + { + int i = p1.X() - _bbox.min.X(); + int z = p1.Z() - _bbox.min.Z(); + VertexIndex index = i+z*_bbox.max.X(); + VertexIndex pos; + if ((pos=_y_cs[index])==-1) + { + _y_cs[index] = (VertexIndex) _mesh->vert.size(); + pos = _y_cs[index]; + Allocator::AddVertices( *_mesh, 1); + v = &_mesh->vert[ pos ]; + _volume->GetYIntercept(p1, p2, v,_thr); + } + v = &_mesh->vert[pos]; + } + void GetZIntercept(const vcg::Point3i &p1, const vcg::Point3i &p2, VertexPointer &v) + { + int i = p1.X() - _bbox.min.X(); + int z = p1.Z() - _bbox.min.Z(); + VertexIndex index = i+z*_bbox.max.X(); + VertexIndex pos; + if (p1.Y()==_current_slice) + { + if ((pos=_z_cs[index])==-1) + { + _z_cs[index] = (VertexIndex) _mesh->vert.size(); + pos = _z_cs[index]; + Allocator::AddVertices( *_mesh, 1 ); + v = &_mesh->vert[pos]; + _volume->GetZIntercept(p1, p2, v,_thr); + return; + } + } + if (p1.Y()==_current_slice+1) + { + if ((pos=_z_ns[index])==-1) + { + _z_ns[index] = (VertexIndex) _mesh->vert.size(); + pos = _z_ns[index]; + Allocator::AddVertices( *_mesh, 1 ); + v = &_mesh->vert[pos]; + _volume->GetZIntercept(p1, p2, v,_thr); + return; + } + } + v = &_mesh->vert[pos]; + } protected: - Box3i _bbox; + Box3i _bbox; - int _slice_dimension; - int _current_slice; + int _slice_dimension; + int _current_slice; - VertexIndex *_x_cs; // indici dell'intersezioni della superficie lungo gli Xedge della fetta corrente - VertexIndex *_y_cs; // indici dell'intersezioni della superficie lungo gli Yedge della fetta corrente - VertexIndex *_z_cs; // indici dell'intersezioni della superficie lungo gli Zedge della fetta corrente - VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta - VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta + VertexIndex *_x_cs; // indici dell'intersezioni della superficie lungo gli Xedge della fetta corrente + VertexIndex *_y_cs; // indici dell'intersezioni della superficie lungo gli Yedge della fetta corrente + VertexIndex *_z_cs; // indici dell'intersezioni della superficie lungo gli Zedge della fetta corrente + VertexIndex *_x_ns; // indici dell'intersezioni della superficie lungo gli Xedge della prossima fetta + VertexIndex *_z_ns; // indici dell'intersezioni della superficie lungo gli Zedge della prossima fetta - MeshType *_mesh; - VolumeType *_volume; + MeshType *_mesh; + VolumeType *_volume; float _thr; void NextSlice() @@ -328,23 +327,23 @@ protected: memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); - std::swap(_x_cs, _x_ns); - std::swap(_z_cs, _z_ns); + std::swap(_x_cs, _x_ns); + std::swap(_z_cs, _z_ns); - _current_slice += 1; - } + _current_slice += 1; + } - void Begin() - { - _current_slice = _bbox.min.Y(); + void Begin() + { + _current_slice = _bbox.min.Y(); - memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex)); - memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_x_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_y_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_cs, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_x_ns, -1, _slice_dimension*sizeof(VertexIndex)); + memset(_z_ns, -1, _slice_dimension*sizeof(VertexIndex)); - } + } }; } // end namespace } // end namespace diff --git a/vcg/complex/algorithms/hole.h b/vcg/complex/algorithms/hole.h index 059599c4..630a0e93 100644 --- a/vcg/complex/algorithms/hole.h +++ b/vcg/complex/algorithms/hole.h @@ -8,7 +8,7 @@ * \ * * All rights reserved. * * * -* This program is free software; you can redistribute it and/or modify * +* 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. * @@ -23,10 +23,7 @@ #ifndef __VCG_TRI_UPDATE_HOLE #define __VCG_TRI_UPDATE_HOLE -#include -#include #include -#include // This file contains three Ear Classes // - TrivialEar @@ -37,30 +34,30 @@ namespace vcg { - namespace tri { + namespace tri { - /* - An ear is identified by TWO pos. - The Three vertexes of an Ear are: - e0.VFlip().v - e0.v - e1.v - Invariants: - e1 == e0.NextB(); - e1.FlipV() == e0; + /* + An ear is identified by TWO pos. + The Three vertexes of an Ear are: + e0.VFlip().v + e0.v + e1.v + Invariants: + e1 == e0.NextB(); + e1.FlipV() == e0; - Situazioni ear non manifold, e degeneri (buco triangolare) + Situazioni ear non manifold, e degeneri (buco triangolare) - T XXXXXXXXXXXXX A /XXXXX B en/XXXXX - /XXXXXXXXXXXXXXX /XXXXXX /XXXXXX - XXXXXXep==en XXX ep\ /en XXXX /e1 XXXX - XXXXXX ----/| XX ------ ----/| XX ------ ----/|XXX - XXXXXX| /e1 XX XXXXXX| /e1 XX XXXXXX| o/e0 XX - XXXXXX| /XXXXXX XXXXXX| /XXXXXX XXXXXX| /XXXXXX - XXX e0|o/XXXXXXX XXX e0|o/XXXXXXX XXX ep| /XXXXXXX - XXX \|/XXXXXXXX XXX \|/XXXXXXXX XXX \|/XXXXXXXX - XXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXX - */ + T XXXXXXXXXXXXX A /XXXXX B en/XXXXX + /XXXXXXXXXXXXXXX /XXXXXX /XXXXXX + XXXXXXep==en XXX ep\ /en XXXX /e1 XXXX + XXXXXX ----/| XX ------ ----/| XX ------ ----/|XXX + XXXXXX| /e1 XX XXXXXX| /e1 XX XXXXXX| o/e0 XX + XXXXXX| /XXXXXX XXXXXX| /XXXXXX XXXXXX| /XXXXXX + XXX e0|o/XXXXXXX XXX e0|o/XXXXXXX XXX ep| /XXXXXXX + XXX \|/XXXXXXXX XXX \|/XXXXXXXX XXX \|/XXXXXXXX + XXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXX + */ template class TrivialEar { public: @@ -328,315 +325,315 @@ template class Hole { public: - typedef typename MESH::VertexType VertexType; - typedef typename MESH::VertexPointer VertexPointer; - typedef typename MESH::ScalarType ScalarType; - typedef typename MESH::FaceType FaceType; - typedef typename MESH::FacePointer FacePointer; - typedef typename MESH::FaceIterator FaceIterator; - typedef typename MESH::CoordType CoordType; + typedef typename MESH::VertexType VertexType; + typedef typename MESH::VertexPointer VertexPointer; + typedef typename MESH::ScalarType ScalarType; + typedef typename MESH::FaceType FaceType; + typedef typename MESH::FacePointer FacePointer; + typedef typename MESH::FaceIterator FaceIterator; + typedef typename MESH::CoordType CoordType; typedef typename vcg::Box3 Box3Type; - typedef typename face::Pos PosType; + typedef typename face::Pos PosType; public: - class Info - { - public: - Info(){} - Info(PosType const &pHole, int const pHoleSize, Box3 &pHoleBB) - { - p=pHole; - size=pHoleSize; - bb=pHoleBB; - } + class Info + { + public: + Info(){} + Info(PosType const &pHole, int const pHoleSize, Box3 &pHoleBB) + { + p=pHole; + size=pHoleSize; + bb=pHoleBB; + } - PosType p; - int size; - Box3Type bb; - - bool operator < (const Info & hh) const {return size < hh.size;} + PosType p; + int size; + Box3Type bb; - ScalarType Perimeter() - { - ScalarType sum=0; - PosType ip = p; - do - { - sum+=Distance(ip.v->cP(),ip.VFlip()->cP()); - ip.NextB(); - } - while (ip != p); - return sum; - } + bool operator < (const Info & hh) const {return size < hh.size;} + + ScalarType Perimeter() + { + ScalarType sum=0; + PosType ip = p; + do + { + sum+=Distance(ip.v->cP(),ip.VFlip()->cP()); + ip.NextB(); + } + while (ip != p); + return sum; + } // Support function to test the validity of a single hole loop // for now it test only that all the edges are border; - // The real test should check if all non manifold vertices + // The real test should check if all non manifold vertices // are touched only by edges belonging to this hole loop. bool CheckValidity() { - if(!p.IsBorder()) + if(!p.IsBorder()) return false; PosType ip=p;ip.NextB(); for(;ip!=p;ip.NextB()) { - if(!ip.IsBorder()) + if(!ip.IsBorder()) return false; - } + } return true; } - }; + }; - class EdgeToBeAvoided - { - VertexPointer v0,v1; - EdgeToBeAvoided(VertexPointer _v0, VertexPointer _v1):v0(_v0),v1(_v1) - { - if(v0>v1) swap(v0,v1); - } - bool operator < (const EdgeToBeAvoided &e) - { - if(this->v0!=e.v0) return this->v0v1v1) swap(v0,v1); + } + bool operator < (const EdgeToBeAvoided &e) + { + if(this->v0!=e.v0) return this->v0v1 - static void FillHoleEar(MESH &m, // The mesh to be filled - Info &h, // the particular hole to be filled - std::vector &facePointersToBeUpdated) - { - //Aggiungo le facce e aggiorno il puntatore alla faccia! - FaceIterator f = tri::Allocator::AddFaces(m, h.size-2, facePointersToBeUpdated); + static void FillHoleEar(MESH &m, // The mesh to be filled + Info &h, // the particular hole to be filled + std::vector &facePointersToBeUpdated) + { + //Aggiungo le facce e aggiorno il puntatore alla faccia! + FaceIterator f = tri::Allocator::AddFaces(m, h.size-2, facePointersToBeUpdated); - assert(h.p.f >= &*m.face.begin()); - assert(h.p.f <= &m.face.back()); - assert(h.p.IsBorder()); + assert(h.p.f >= &*m.face.begin()); + assert(h.p.f <= &m.face.back()); + assert(h.p.IsBorder()); - std::vector< EAR > EarHeap; - EarHeap.reserve(h.size); - int nmBit= VertexType::NewBitFlag(); // non manifoldness bit + std::vector< EAR > EarHeap; + EarHeap.reserve(h.size); + int nmBit= VertexType::NewBitFlag(); // non manifoldness bit - //First loops around the hole to mark non manifold vertices. - PosType ip = h.p; // Pos iterator - do{ - ip.V()->ClearUserBit(nmBit); - ip.V()->ClearV(); - ip.NextB(); - } while(ip!=h.p); + //First loops around the hole to mark non manifold vertices. + PosType ip = h.p; // Pos iterator + do{ + ip.V()->ClearUserBit(nmBit); + ip.V()->ClearV(); + ip.NextB(); + } while(ip!=h.p); - ip = h.p; // Re init the pos iterator for another loop (useless if everithing is ok!!) - do{ - if(!ip.V()->IsV()) - ip.V()->SetV(); // All the vertexes that are visited more than once are non manifold - else ip.V()->SetUserBit(nmBit); - ip.NextB(); - } while(ip!=h.p); + ip = h.p; // Re init the pos iterator for another loop (useless if everithing is ok!!) + do{ + if(!ip.V()->IsV()) + ip.V()->SetV(); // All the vertexes that are visited more than once are non manifold + else ip.V()->SetUserBit(nmBit); + ip.NextB(); + } while(ip!=h.p); - PosType fp = h.p; - do{ - EAR appEar = EAR(fp); - EarHeap.push_back( appEar ); - //printf("Adding ear %s ",app.Dump()); - fp.NextB(); - assert(fp.IsBorder()); - }while(fp!=h.p); + PosType fp = h.p; + do{ + EAR appEar = EAR(fp); + EarHeap.push_back( appEar ); + //printf("Adding ear %s ",app.Dump()); + fp.NextB(); + assert(fp.IsBorder()); + }while(fp!=h.p); - int cnt=h.size; + int cnt=h.size; - make_heap(EarHeap.begin(), EarHeap.end()); + make_heap(EarHeap.begin(), EarHeap.end()); - //finche' il buco non e' chiuso o non ci sono piu' orecchie da analizzare. - while( cnt > 2 && !EarHeap.empty() ) - { - //printf("Front of the heap is %s", H.front().Dump()); - pop_heap(EarHeap.begin(), EarHeap.end()); // retrieve the MAXIMUM value and put in the back; - EAR BestEar=EarHeap.back(); - EarHeap.pop_back(); + //finche' il buco non e' chiuso o non ci sono piu' orecchie da analizzare. + while( cnt > 2 && !EarHeap.empty() ) + { + //printf("Front of the heap is %s", H.front().Dump()); + pop_heap(EarHeap.begin(), EarHeap.end()); // retrieve the MAXIMUM value and put in the back; + EAR BestEar=EarHeap.back(); + EarHeap.pop_back(); - if(BestEar.IsUpToDate() && !BestEar.IsDegen(nmBit)) - { - if((*f).HasPolyInfo()) (*f).Alloc(3); - PosType ep0,ep1; - if(BestEar.Close(ep0,ep1,&*f)) - { - if(!ep0.IsNull()){ - EarHeap.push_back(EAR(ep0)); - push_heap( EarHeap.begin(), EarHeap.end()); - } - if(!ep1.IsNull()){ - EarHeap.push_back(EAR(ep1)); - push_heap( EarHeap.begin(), EarHeap.end()); - } - --cnt; - ++f; - } - }//is update() - }//fine del while principale. + if(BestEar.IsUpToDate() && !BestEar.IsDegen(nmBit)) + { + if((*f).HasPolyInfo()) (*f).Alloc(3); + PosType ep0,ep1; + if(BestEar.Close(ep0,ep1,&*f)) + { + if(!ep0.IsNull()){ + EarHeap.push_back(EAR(ep0)); + push_heap( EarHeap.begin(), EarHeap.end()); + } + if(!ep1.IsNull()){ + EarHeap.push_back(EAR(ep1)); + push_heap( EarHeap.begin(), EarHeap.end()); + } + --cnt; + ++f; + } + }//is update() + }//fine del while principale. - while(f!=m.face.end()){ - tri::Allocator::DeleteFace(m,*f); - f++; - } + while(f!=m.face.end()){ + tri::Allocator::DeleteFace(m,*f); + f++; + } - VertexType::DeleteBitFlag(nmBit); // non manifoldness bit - } + VertexType::DeleteBitFlag(nmBit); // non manifoldness bit + } - template - static int EarCuttingFill(MESH &m, int sizeHole, bool Selected = false, CallBackPos *cb=0) - { - std::vector< Info > vinfo; - GetInfo(m, Selected,vinfo); + template + static int EarCuttingFill(MESH &m, int sizeHole, bool Selected = false, CallBackPos *cb=0) + { + std::vector< Info > vinfo; + GetInfo(m, Selected,vinfo); - typename std::vector::iterator ith; - int indCb=0; - int holeCnt=0; - std::vector facePtrToBeUpdated; - for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) - facePtrToBeUpdated.push_back( &(*ith).p.f ); + typename std::vector::iterator ith; + int indCb=0; + int holeCnt=0; + std::vector facePtrToBeUpdated; + for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) + facePtrToBeUpdated.push_back( &(*ith).p.f ); - for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) - { - indCb++; - if(cb) (*cb)(indCb*10/vinfo.size(),"Closing Holes"); - if((*ith).size < sizeHole){ - holeCnt++; - FillHoleEar< EAR >(m, *ith,facePtrToBeUpdated); - } - } - return holeCnt; - } + for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) + { + indCb++; + if(cb) (*cb)(indCb*10/vinfo.size(),"Closing Holes"); + if((*ith).size < sizeHole){ + holeCnt++; + FillHoleEar< EAR >(m, *ith,facePtrToBeUpdated); + } + } + return holeCnt; + } /// Main Hole Filling function. /// Given a mesh search for all the holes smaller than a given size and fill them /// It returns the number of filled holes. template - static int EarCuttingIntersectionFill(MESH &m, const int maxSizeHole, bool Selected, CallBackPos *cb=0) - { - std::vector vinfo; - GetInfo(m, Selected,vinfo); - typename std::vector::iterator ith; + static int EarCuttingIntersectionFill(MESH &m, const int maxSizeHole, bool Selected, CallBackPos *cb=0) + { + std::vector vinfo; + GetInfo(m, Selected,vinfo); + typename std::vector::iterator ith; - // collect the face pointer that has to be updated by the various addfaces - std::vector vfpOrig; - for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) - vfpOrig.push_back( &(*ith).p.f ); + // collect the face pointer that has to be updated by the various addfaces + std::vector vfpOrig; + for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) + vfpOrig.push_back( &(*ith).p.f ); - int indCb=0; - int holeCnt=0; - for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) - { - indCb++; - if(cb) (*cb)(indCb*10/vinfo.size(),"Closing Holes"); - if((*ith).size < maxSizeHole){ - std::vector facePtrToBeUpdated; - holeCnt++; - facePtrToBeUpdated=vfpOrig; - EAR::AdjacencyRing().clear(); - //Loops around the hole to collect the faces that have to be tested for intersection. - PosType ip = (*ith).p; - do - { - PosType inp = ip; - do - { - inp.FlipE(); - inp.FlipF(); - EAR::AdjacencyRing().push_back(inp.f); - } while(!inp.IsBorder()); - ip.NextB(); - }while(ip != (*ith).p); + int indCb=0; + int holeCnt=0; + for(ith = vinfo.begin(); ith!= vinfo.end(); ++ith) + { + indCb++; + if(cb) (*cb)(indCb*10/vinfo.size(),"Closing Holes"); + if((*ith).size < maxSizeHole){ + std::vector facePtrToBeUpdated; + holeCnt++; + facePtrToBeUpdated=vfpOrig; + EAR::AdjacencyRing().clear(); + //Loops around the hole to collect the faces that have to be tested for intersection. + PosType ip = (*ith).p; + do + { + PosType inp = ip; + do + { + inp.FlipE(); + inp.FlipF(); + EAR::AdjacencyRing().push_back(inp.f); + } while(!inp.IsBorder()); + ip.NextB(); + }while(ip != (*ith).p); - typename std::vector::iterator fpi; - for(fpi=EAR::AdjacencyRing().begin();fpi!=EAR::AdjacencyRing().end();++fpi) - facePtrToBeUpdated.push_back( &*fpi ); + typename std::vector::iterator fpi; + for(fpi=EAR::AdjacencyRing().begin();fpi!=EAR::AdjacencyRing().end();++fpi) + facePtrToBeUpdated.push_back( &*fpi ); - FillHoleEar(m, *ith,facePtrToBeUpdated); - EAR::AdjacencyRing().clear(); - } - } - return holeCnt; - } + FillHoleEar(m, *ith,facePtrToBeUpdated); + EAR::AdjacencyRing().clear(); + } + } + return holeCnt; + } - static void GetInfo(MESH &m, bool Selected ,std::vector& VHI) - { - tri::UpdateFlags::FaceClearV(m); - for(FaceIterator fi = m.face.begin(); fi!=m.face.end(); ++fi) - { - if(!(*fi).IsD()) - { - if(Selected && !(*fi).IsS()) - { - //se devo considerare solo i triangoli selezionati e - //quello che sto considerando non lo e' lo marchio e vado avanti - (*fi).SetV(); - } - else - { - for(int j =0; j<3 ; ++j) - { - if( face::IsBorder(*fi,j) && !(*fi).IsV() ) - {//Trovato una faccia di bordo non ancora visitata. - (*fi).SetV(); - PosType sp(&*fi, j, (*fi).V(j)); - PosType fp=sp; - int holesize=0; + static void GetInfo(MESH &m, bool Selected ,std::vector& VHI) + { + tri::UpdateFlags::FaceClearV(m); + for(FaceIterator fi = m.face.begin(); fi!=m.face.end(); ++fi) + { + if(!(*fi).IsD()) + { + if(Selected && !(*fi).IsS()) + { + //se devo considerare solo i triangoli selezionati e + //quello che sto considerando non lo e' lo marchio e vado avanti + (*fi).SetV(); + } + else + { + for(int j =0; j<3 ; ++j) + { + if( face::IsBorder(*fi,j) && !(*fi).IsV() ) + {//Trovato una faccia di bordo non ancora visitata. + (*fi).SetV(); + PosType sp(&*fi, j, (*fi).V(j)); + PosType fp=sp; + int holesize=0; - Box3Type hbox; - hbox.Add(sp.v->cP()); + Box3Type hbox; + hbox.Add(sp.v->cP()); //printf("Looping %i : (face %i edge %i) \n", VHI.size(),sp.f-&*m.face.begin(),sp.z); sp.f->SetV(); - do - { - sp.f->SetV(); - hbox.Add(sp.v->cP()); - ++holesize; - sp.NextB(); - sp.f->SetV(); - assert(sp.IsBorder()); - }while(sp != fp); + do + { + sp.f->SetV(); + hbox.Add(sp.v->cP()); + ++holesize; + sp.NextB(); + sp.f->SetV(); + assert(sp.IsBorder()); + }while(sp != fp); - //ho recuperato l'inofrmazione su tutto il buco - VHI.push_back( Info(sp,holesize,hbox) ); - } - }//for sugli edge del triangolo - }//S & !S - }//!IsD() - }//for principale!!! - } + //ho recuperato l'inofrmazione su tutto il buco + VHI.push_back( Info(sp,holesize,hbox) ); + } + }//for sugli edge del triangolo + }//S & !S + }//!IsD() + }//for principale!!! + } - //Minimum Weight Algorithm - class Weight - { - public: + //Minimum Weight Algorithm + class Weight + { + public: - Weight() { ang = 180; ar = FLT_MAX ;} - Weight( float An, float Ar ) { ang=An ; ar= Ar;} - ~Weight() {} + Weight() { ang = 180; ar = FLT_MAX ;} + Weight( float An, float Ar ) { ang=An ; ar= Ar;} + ~Weight() {} - float angle() const { return ang; } - float area() const { return ar; } + float angle() const { return ang; } + float area() const { return ar; } - Weight operator+( const Weight & other ) const {return Weight( std::max( angle(), other.angle() ), area() + other.area());} - bool operator<( const Weight & rhs ) const {return ( angle() < rhs.angle() ||(angle() == rhs.angle() && area() < rhs.area())); } + Weight operator+( const Weight & other ) const {return Weight( std::max( angle(), other.angle() ), area() + other.area());} + bool operator<( const Weight & rhs ) const {return ( angle() < rhs.angle() ||(angle() == rhs.angle() && area() < rhs.area())); } - private: - float ang; - float ar; - }; + private: + float ang; + float ar; + }; - /* - \ / \/ - v1*---------*v4 + /* + \ / \/ + v1*---------*v4 / \ / / \ / / \ / @@ -644,235 +641,235 @@ template *---------*- | v3 v2\ */ - - static float ComputeDihedralAngle(CoordType p1,CoordType p2,CoordType p3,CoordType p4) - { - CoordType n1 = NormalizedNormal(p1,p3,p2); - CoordType n2 = NormalizedNormal(p1,p2,p4); - return math::ToDeg(AngleN(n1,n2)); - } + + static float ComputeDihedralAngle(CoordType p1,CoordType p2,CoordType p3,CoordType p4) + { + CoordType n1 = NormalizedNormal(p1,p3,p2); + CoordType n2 = NormalizedNormal(p1,p2,p4); + return math::ToDeg(AngleN(n1,n2)); + } static bool existEdge(PosType pi,PosType pf) - { - PosType app = pi; - PosType appF = pi; - PosType tmp; - assert(pi.IsBorder()); - appF.NextB(); - appF.FlipV(); - do - { - tmp = app; - tmp.FlipV(); - if(tmp.v == pf.v) - return true; - app.FlipE(); - app.FlipF(); + { + PosType app = pi; + PosType appF = pi; + PosType tmp; + assert(pi.IsBorder()); + appF.NextB(); + appF.FlipV(); + do + { + tmp = app; + tmp.FlipV(); + if(tmp.v == pf.v) + return true; + app.FlipE(); + app.FlipF(); - if(app == pi)return false; - }while(app != appF); - return false; - } + if(app == pi)return false; + }while(app != appF); + return false; + } - static Weight computeWeight( int i, int j, int k, - std::vector pv, - std::vector< std::vector< int > > v) - { - PosType pi = pv[i]; - PosType pj = pv[j]; - PosType pk = pv[k]; + static Weight computeWeight( int i, int j, int k, + std::vector pv, + std::vector< std::vector< int > > v) + { + PosType pi = pv[i]; + PosType pj = pv[j]; + PosType pk = pv[k]; - //test complex edge - if(existEdge(pi,pj) || existEdge(pj,pk)|| existEdge(pk,pi) ) - { - return Weight(); - } - // Return an infinite weight, if one of the neighboring patches - // could not be created. - if(v[i][j] == -1){return Weight();} - if(v[j][k] == -1){return Weight();} + //test complex edge + if(existEdge(pi,pj) || existEdge(pj,pk)|| existEdge(pk,pi) ) + { + return Weight(); + } + // Return an infinite weight, if one of the neighboring patches + // could not be created. + if(v[i][j] == -1){return Weight();} + if(v[j][k] == -1){return Weight();} - //calcolo il massimo angolo diedrale, se esiste. - float angle = 0.0f; - PosType px; - if(i + 1 == j) - { - px = pj; - px.FlipE(); px.FlipV(); - angle = std::max(angle , ComputeDihedralAngle(pi.v->P(), pj.v->P(), pk.v->P(), px.v->P()) ); - } - else - { - angle = std::max( angle, ComputeDihedralAngle(pi.v->P(),pj.v->P(), pk.v->P(), pv[ v[i][j] ].v->P())); - } + //calcolo il massimo angolo diedrale, se esiste. + float angle = 0.0f; + PosType px; + if(i + 1 == j) + { + px = pj; + px.FlipE(); px.FlipV(); + angle = std::max(angle , ComputeDihedralAngle(pi.v->P(), pj.v->P(), pk.v->P(), px.v->P()) ); + } + else + { + angle = std::max( angle, ComputeDihedralAngle(pi.v->P(),pj.v->P(), pk.v->P(), pv[ v[i][j] ].v->P())); + } - if(j + 1 == k) - { - px = pk; - px.FlipE(); px.FlipV(); - angle = std::max(angle , ComputeDihedralAngle(pj.v->P(), pk.v->P(), pi.v->P(), px.v->P()) ); - } - else - { - angle = std::max( angle, ComputeDihedralAngle(pj.v->P(),pk.v->P(), pi.v->P(), pv[ v[j][k] ].v->P())); - } + if(j + 1 == k) + { + px = pk; + px.FlipE(); px.FlipV(); + angle = std::max(angle , ComputeDihedralAngle(pj.v->P(), pk.v->P(), pi.v->P(), px.v->P()) ); + } + else + { + angle = std::max( angle, ComputeDihedralAngle(pj.v->P(),pk.v->P(), pi.v->P(), pv[ v[j][k] ].v->P())); + } - if( i == 0 && k == (int)v.size() - 1) - { - px = pi; - px.FlipE(); px.FlipV(); - angle = std::max(angle , ComputeDihedralAngle(pk.v->P(), pi.v->P(), pj.v->P(),px.v->P() ) ); - } + if( i == 0 && k == (int)v.size() - 1) + { + px = pi; + px.FlipE(); px.FlipV(); + angle = std::max(angle , ComputeDihedralAngle(pk.v->P(), pi.v->P(), pj.v->P(),px.v->P() ) ); + } - ScalarType area = ( (pj.v->P() - pi.v->P()) ^ (pk.v->P() - pi.v->P()) ).Norm() * 0.5; + ScalarType area = ( (pj.v->P() - pi.v->P()) ^ (pk.v->P() - pi.v->P()) ).Norm() * 0.5; - return Weight(angle, area); - } + return Weight(angle, area); + } - static void calculateMinimumWeightTriangulation(MESH &m, FaceIterator f,std::vector vv ) - { - std::vector< std::vector< Weight > > w; //matrice dei pesi minimali di ogni orecchio preso in conzideraione - std::vector< std::vector< int > > vi;//memorizza l'indice del terzo vertice del triangolo + static void calculateMinimumWeightTriangulation(MESH &m, FaceIterator f,std::vector vv ) + { + std::vector< std::vector< Weight > > w; //matrice dei pesi minimali di ogni orecchio preso in conzideraione + std::vector< std::vector< int > > vi;//memorizza l'indice del terzo vertice del triangolo - //hole size - int nv = vv.size(); - - w.clear(); - w.resize( nv, std::vector( nv, Weight() ) ); + //hole size + int nv = vv.size(); - vi.resize( nv, std::vector( nv, 0 ) ); + w.clear(); + w.resize( nv, std::vector( nv, Weight() ) ); - //inizializzo tutti i pesi possibili del buco - for ( int i = 0; i < nv-1; ++i ) - w[i][i+1] = Weight( 0, 0 ); + vi.resize( nv, std::vector( nv, 0 ) ); - //doppio ciclo for per calcolare di tutti i possibili triangoli i loro pesi. - for ( int j = 2; j < nv; ++j ) - { - for ( int i = 0; i + j < nv; ++i ) - { - //per ogni triangolazione mi mantengo il minimo valore del peso tra i triangoli possibili - Weight minval; + //inizializzo tutti i pesi possibili del buco + for ( int i = 0; i < nv-1; ++i ) + w[i][i+1] = Weight( 0, 0 ); - //indice del vertice che da il peso minimo nella triangolazione corrente - int minIndex = -1; + //doppio ciclo for per calcolare di tutti i possibili triangoli i loro pesi. + for ( int j = 2; j < nv; ++j ) + { + for ( int i = 0; i + j < nv; ++i ) + { + //per ogni triangolazione mi mantengo il minimo valore del peso tra i triangoli possibili + Weight minval; - //ciclo tra i vertici in mezzo a i due prefissati - for ( int m = i + 1; m < i + j; ++m ) - { - Weight a = w[i][m]; - Weight b = w[m][i+j]; - Weight newval = a + b + computeWeight( i, m, i+j, vv, vi); - if ( newval < minval ) - { - minval = newval; - minIndex = m; - } - } - w[i][i+j] = minval; - vi[i][i+j] = minIndex; - } - } + //indice del vertice che da il peso minimo nella triangolazione corrente + int minIndex = -1; - //Triangulate - int i, j; - i=0; j=nv-1; - - triangulate(m,f, i, j, vi, vv); + //ciclo tra i vertici in mezzo a i due prefissati + for ( int m = i + 1; m < i + j; ++m ) + { + Weight a = w[i][m]; + Weight b = w[m][i+j]; + Weight newval = a + b + computeWeight( i, m, i+j, vv, vi); + if ( newval < minval ) + { + minval = newval; + minIndex = m; + } + } + w[i][i+j] = minval; + vi[i][i+j] = minIndex; + } + } - while(f!=m.face.end()) - { - (*f).SetD(); - ++f; - m.fn--; - } - } + //Triangulate + int i, j; + i=0; j=nv-1; + + triangulate(m,f, i, j, vi, vv); + + while(f!=m.face.end()) + { + (*f).SetD(); + ++f; + m.fn--; + } + } - static void triangulate(MESH &m, FaceIterator &f,int i, int j, + static void triangulate(MESH &m, FaceIterator &f,int i, int j, std::vector< std::vector > vi, std::vector vv) - { - if(i + 1 == j){return;} - if(i==j)return; + { + if(i + 1 == j){return;} + if(i==j)return; - int k = vi[i][j]; + int k = vi[i][j]; - if(k == -1) return; + if(k == -1) return; - //Setto i vertici - f->V(0) = vv[i].v; - f->V(1) = vv[k].v; - f->V(2) = vv[j].v; - - f++; - triangulate(m,f,i,k,vi,vv); - triangulate(m,f,k,j,vi,vv); - } + //Setto i vertici + f->V(0) = vv[i].v; + f->V(1) = vv[k].v; + f->V(2) = vv[j].v; + + f++; + triangulate(m,f,i,k,vi,vv); + triangulate(m,f,k,j,vi,vv); + } static void MinimumWeightFill(MESH &m, int holeSize, bool Selected) - { - std::vector vvi; - std::vector vfp; + { + std::vector vvi; + std::vector vfp; - std::vector vinfo; - typename std::vector::iterator VIT; - GetInfo(m, Selected,vinfo); + std::vector vinfo; + typename std::vector::iterator VIT; + GetInfo(m, Selected,vinfo); - for(VIT = vinfo.begin(); VIT != vinfo.end();++VIT) - { - vvi.push_back(VIT->p); - } + for(VIT = vinfo.begin(); VIT != vinfo.end();++VIT) + { + vvi.push_back(VIT->p); + } - typename std::vector::iterator ith; - typename std::vector::iterator ithn; - typename std::vector::iterator itf; + typename std::vector::iterator ith; + typename std::vector::iterator ithn; + typename std::vector::iterator itf; - std::vector app; - PosType ps; - std::vector tr; - std::vector vf; + std::vector app; + PosType ps; + std::vector tr; + std::vector vf; - for(ith = vvi.begin(); ith!= vvi.end(); ++ith) - { - tr.clear(); - vf.clear(); - app.clear(); - vfp.clear(); + for(ith = vvi.begin(); ith!= vvi.end(); ++ith) + { + tr.clear(); + vf.clear(); + app.clear(); + vfp.clear(); - ps = *ith; - getBoundHole(ps,app); + ps = *ith; + getBoundHole(ps,app); if(app.size() <= size_t(holeSize) ) - { - typename std::vector::iterator itP; - std::vector vfp; + { + typename std::vector::iterator itP; + std::vector vfp; - for(ithn = vvi.begin(); ithn!= vvi.end(); ++ithn) - vfp.push_back(&(ithn->f)); + for(ithn = vvi.begin(); ithn!= vvi.end(); ++ithn) + vfp.push_back(&(ithn->f)); - for(itP = app.begin (); itP != app.end ();++itP) - vfp.push_back( &(*itP).f ); + for(itP = app.begin (); itP != app.end ();++itP) + vfp.push_back( &(*itP).f ); - //aggiungo le facce - FaceIterator f = tri::Allocator::AddFaces(m, (app.size()-2) , vfp); + //aggiungo le facce + FaceIterator f = tri::Allocator::AddFaces(m, (app.size()-2) , vfp); - calculateMinimumWeightTriangulation(m,f, app); - } - } + calculateMinimumWeightTriangulation(m,f, app); + } + } - } + } - static void getBoundHole (PosType sp,std::vector&ret) - { - PosType fp = sp; - //take vertex around the hole - do - { - assert(fp.IsBorder()); - ret.push_back(fp); - fp.NextB(); - }while(sp != fp); - } + static void getBoundHole (PosType sp,std::vector&ret) + { + PosType fp = sp; + //take vertex around the hole + do + { + assert(fp.IsBorder()); + ret.push_back(fp); + fp.NextB(); + }while(sp != fp); + } };//close class Hole diff --git a/vcg/complex/algorithms/refine.h b/vcg/complex/algorithms/refine.h index 0138b7c6..2d333d1a 100644 --- a/vcg/complex/algorithms/refine.h +++ b/vcg/complex/algorithms/refine.h @@ -33,7 +33,6 @@ #include #include #include -#include namespace vcg{ namespace tri{ diff --git a/vcg/complex/algorithms/smooth.h b/vcg/complex/algorithms/smooth.h index 96123182..85889218 100644 --- a/vcg/complex/algorithms/smooth.h +++ b/vcg/complex/algorithms/smooth.h @@ -26,10 +26,7 @@ #define __VCGLIB__SMOOTH #include -#include -#include #include -#include #include #include #include @@ -50,24 +47,24 @@ class Smooth { public: - typedef SmoothMeshType MeshType; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::VertexType::CoordType CoordType; - typedef typename MeshType::VertexPointer VertexPointer; - typedef typename MeshType::VertexIterator VertexIterator; - typedef typename MeshType::ScalarType ScalarType; - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::FacePointer FacePointer; - typedef typename MeshType::FaceIterator FaceIterator; - typedef typename MeshType::FaceContainer FaceContainer; - typedef typename vcg::Box3 Box3Type; - typedef typename vcg::face::VFIterator VFLocalIterator; + typedef SmoothMeshType MeshType; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexType::CoordType CoordType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::ScalarType ScalarType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::FaceContainer FaceContainer; + typedef typename vcg::Box3 Box3Type; + typedef typename vcg::face::VFIterator VFLocalIterator; class ScaleLaplacianInfo { public: - CoordType PntSum; - ScalarType LenSum; + CoordType PntSum; + ScalarType LenSum; }; // This is precisely what curvature flow does. @@ -83,40 +80,40 @@ void VertexCoordLaplacianCurvatureFlow(MeshType &/*m*/, int /*step*/, ScalarType static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType delta) { - ScaleLaplacianInfo lpz; - lpz.PntSum=CoordType(0,0,0); - lpz.LenSum=0; - SimpleTempData TD(m.vert,lpz); - FaceIterator fi; - for(int i=0;iP() + (*fi).V(1)->P() + (*fi).V(2)->P())/3.0; - CoordType e0=((*fi).V(0)->P() - (*fi).V(1)->P()).Normalize(); - CoordType e1=((*fi).V(1)->P() - (*fi).V(2)->P()).Normalize(); - CoordType e2=((*fi).V(2)->P() - (*fi).V(0)->P()).Normalize(); + ScaleLaplacianInfo lpz; + lpz.PntSum=CoordType(0,0,0); + lpz.LenSum=0; + SimpleTempData TD(m.vert,lpz); + FaceIterator fi; + for(int i=0;iP() + (*fi).V(1)->P() + (*fi).V(2)->P())/3.0; + CoordType e0=((*fi).V(0)->P() - (*fi).V(1)->P()).Normalize(); + CoordType e1=((*fi).V(1)->P() - (*fi).V(2)->P()).Normalize(); + CoordType e2=((*fi).V(2)->P() - (*fi).V(0)->P()).Normalize(); - a[0]=AngleN(-e0,e2); - a[1]=AngleN(-e1,e0); - a[2]=AngleN(-e2,e1); - //assert(fabs(M_PI -a[0] -a[1] -a[2])<0.0000001); + a[0]=AngleN(-e0,e2); + a[1]=AngleN(-e1,e0); + a[2]=AngleN(-e2,e1); + //assert(fabs(M_PI -a[0] -a[1] -a[2])<0.0000001); - for(int j=0;j<3;++j){ - CoordType dir= (mp-(*fi).V(j)->P()).Normalize(); - TD[(*fi).V(j)].PntSum+=dir*a[j]; - TD[(*fi).V(j)].LenSum+=a[j]; // well, it should be named angleSum - } - } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].LenSum>0 ) - (*vi).P() = (*vi).P() + (TD[*vi].PntSum/TD[*vi].LenSum ) * delta; + for(int j=0;j<3;++j){ + CoordType dir= (mp-(*fi).V(j)->P()).Normalize(); + TD[(*fi).V(j)].PntSum+=dir*a[j]; + TD[(*fi).V(j)].LenSum+=a[j]; // well, it should be named angleSum + } + } + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].LenSum>0 ) + (*vi).P() = (*vi).P() + (TD[*vi].PntSum/TD[*vi].LenSum ) * delta; - } + } }; // Scale dependent laplacian smoothing [Fujiwara 95] @@ -131,72 +128,72 @@ static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType static void VertexCoordScaleDependentLaplacian_Fujiwara(MeshType &m, int step, ScalarType delta) { - SimpleTempData TD(m.vert); - ScaleLaplacianInfo lpz; - lpz.PntSum=CoordType(0,0,0); - lpz.LenSum=0; - FaceIterator fi; - for(int i=0;i TD(m.vert); + ScaleLaplacianInfo lpz; + lpz.PntSum=CoordType(0,0,0); + lpz.LenSum=0; + FaceIterator fi; + for(int i=0;iP() -(*fi).V(j)->P(); - ScalarType len=Norm(edge); - edge/=len; - TD[(*fi).V(j)].PntSum+=edge; - TD[(*fi).V1(j)].PntSum-=edge; - TD[(*fi).V(j)].LenSum+=len; - TD[(*fi).V1(j)].LenSum+=len; - } + for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) { + CoordType edge= (*fi).V1(j)->P() -(*fi).V(j)->P(); + ScalarType len=Norm(edge); + edge/=len; + TD[(*fi).V(j)].PntSum+=edge; + TD[(*fi).V1(j)].PntSum-=edge; + TD[(*fi).V(j)].LenSum+=len; + TD[(*fi).V1(j)].LenSum+=len; + } - for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD()) - for(int j=0;j<3;++j) - // se l'edge j e' di bordo si riazzera tutto e si riparte - if((*fi).IsB(j)) { - TD[(*fi).V(j)].PntSum=CoordType(0,0,0); - TD[(*fi).V1(j)].PntSum=CoordType(0,0,0); - TD[(*fi).V(j)].LenSum=0; - TD[(*fi).V1(j)].LenSum=0; - } + for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD()) + for(int j=0;j<3;++j) + // se l'edge j e' di bordo si riazzera tutto e si riparte + if((*fi).IsB(j)) { + TD[(*fi).V(j)].PntSum=CoordType(0,0,0); + TD[(*fi).V1(j)].PntSum=CoordType(0,0,0); + TD[(*fi).V(j)].LenSum=0; + TD[(*fi).V1(j)].LenSum=0; + } - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - CoordType edge= (*fi).V1(j)->P() -(*fi).V(j)->P(); - ScalarType len=Norm(edge); - edge/=len; - TD[(*fi).V(j)].PntSum+=edge; - TD[(*fi).V1(j)].PntSum-=edge; - TD[(*fi).V(j)].LenSum+=len; - TD[(*fi).V1(j)].LenSum+=len; - } + for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + CoordType edge= (*fi).V1(j)->P() -(*fi).V(j)->P(); + ScalarType len=Norm(edge); + edge/=len; + TD[(*fi).V(j)].PntSum+=edge; + TD[(*fi).V1(j)].PntSum-=edge; + TD[(*fi).V(j)].LenSum+=len; + TD[(*fi).V1(j)].LenSum+=len; + } // The fundamental part: // We move the new point of a quantity // // L(M) = 1/Sum(edgelen) * Sum(Normalized edges) // - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].LenSum>0 ) - (*vi).P() = (*vi).P() + (TD[*vi].PntSum/TD[*vi].LenSum)*delta; - } + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].LenSum>0 ) + (*vi).P() = (*vi).P() + (TD[*vi].PntSum/TD[*vi].LenSum)*delta; + } }; class LaplacianInfo { public: - LaplacianInfo(const CoordType &_p, const int _n):sum(_p),cnt(_n) {} - LaplacianInfo() {} - CoordType sum; - ScalarType cnt; + LaplacianInfo(const CoordType &_p, const int _n):sum(_p),cnt(_n) {} + LaplacianInfo() {} + CoordType sum; + ScalarType cnt; }; // Classical Laplacian Smoothing. Each vertex can be moved onto the average of the adjacent vertices. @@ -225,39 +222,39 @@ static void AccumulateLaplacianInfo(MeshType &m, SimpleTempDataP(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->P(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } - } + // se l'edge j e' di bordo si deve mediare solo con gli adiacenti + for(fi=m.face.begin();fi!=m.face.end();++fi) + { + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->P(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->P(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } + } } static void VertexCoordLaplacian(MeshType &m, int step, bool SmoothSelected=false, bool cotangentWeight=false, vcg::CallBackPos * cb=0) @@ -298,60 +295,60 @@ static void VertexCoordPlanarLaplacian(MeshType &m, int step, float AngleThrRad TD[*vi].sum = ( (*vi).P() + TD[*vi].sum)/(TD[*vi].cnt+1); } - for(fi=m.face.begin();fi!=m.face.end();++fi){ - if(!(*fi).IsD()){ - for (int j = 0; j < 3; ++j) { - if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, (*fi).P1(j), (*fi).P2(j) ), - NormalizedNormal( (*fi).P0(j) , (*fi).P1(j), (*fi).P2(j) ) ) > AngleThrRad ) - TD[(*fi).V0(j)].sum = (*fi).P0(j); - } - } - } - for(fi=m.face.begin();fi!=m.face.end();++fi){ - if(!(*fi).IsD()){ - for (int j = 0; j < 3; ++j) { - if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, TD[(*fi).V1(j)].sum, (*fi).P2(j) ), - NormalizedNormal( (*fi).P0(j) , (*fi).P1(j), (*fi).P2(j) ) ) > AngleThrRad ) - { - TD[(*fi).V0(j)].sum = (*fi).P0(j); - TD[(*fi).V1(j)].sum = (*fi).P1(j); - } - } - } - } + for(fi=m.face.begin();fi!=m.face.end();++fi){ + if(!(*fi).IsD()){ + for (int j = 0; j < 3; ++j) { + if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, (*fi).P1(j), (*fi).P2(j) ), + NormalizedNormal( (*fi).P0(j) , (*fi).P1(j), (*fi).P2(j) ) ) > AngleThrRad ) + TD[(*fi).V0(j)].sum = (*fi).P0(j); + } + } + } + for(fi=m.face.begin();fi!=m.face.end();++fi){ + if(!(*fi).IsD()){ + for (int j = 0; j < 3; ++j) { + if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, TD[(*fi).V1(j)].sum, (*fi).P2(j) ), + NormalizedNormal( (*fi).P0(j) , (*fi).P1(j), (*fi).P2(j) ) ) > AngleThrRad ) + { + TD[(*fi).V0(j)].sum = (*fi).P0(j); + TD[(*fi).V1(j)].sum = (*fi).P1(j); + } + } + } + } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - (*vi).P()= TD[*vi].sum; + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + (*vi).P()= TD[*vi].sum; - }// end step + }// end step } static void VertexCoordLaplacianBlend(MeshType &m, int step, float alpha, bool SmoothSelected=false) { - VertexIterator vi; - LaplacianInfo lpz(CoordType(0,0,0),0); - assert (alpha<= 1.0); - SimpleTempData TD(m.vert); + VertexIterator vi; + LaplacianInfo lpz(CoordType(0,0,0),0); + assert (alpha<= 1.0); + SimpleTempData TD(m.vert); - for(int i=0;i0 ) - { - if(!SmoothSelected || (*vi).IsS()) - { - CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); - (*vi).P() = (*vi).P() + Delta*alpha; - } - } - } + for(int i=0;i0 ) + { + if(!SmoothSelected || (*vi).IsS()) + { + CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); + (*vi).P() = (*vi).P() + Delta*alpha; + } + } + } } /* a couple of notes about the lambda mu values @@ -359,7 +356,7 @@ We assume that 0 < lambda , and mu is a negative scale factor such that mu < - Holds mu+lambda < 0 (e.g in absolute value mu is greater) let kpb be the pass-band frequency, taubin says that: - kpb = 1/lambda + 1/mu >0 + kpb = 1/lambda + 1/mu >0 Values of kpb from 0.01 to 0.1 produce good results according to the original paper. @@ -374,35 +371,35 @@ So if static void VertexCoordTaubin(MeshType &m, int step, float lambda, float mu, bool SmoothSelected=false, vcg::CallBackPos * cb=0) { - LaplacianInfo lpz(CoordType(0,0,0),0); - SimpleTempData TD(m.vert,lpz); - VertexIterator vi; - for(int i=0;i0 ) - { - if(!SmoothSelected || (*vi).IsS()) - { - CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); - (*vi).P() = (*vi).P() + Delta*lambda ; - } - } - TD.Init(lpz); - AccumulateLaplacianInfo(m,TD); - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - { - if(!SmoothSelected || (*vi).IsS()) - { - CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); - (*vi).P() = (*vi).P() + Delta*mu ; - } - } - } // end for step + LaplacianInfo lpz(CoordType(0,0,0),0); + SimpleTempData TD(m.vert,lpz); + VertexIterator vi; + for(int i=0;i0 ) + { + if(!SmoothSelected || (*vi).IsS()) + { + CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); + (*vi).P() = (*vi).P() + Delta*lambda ; + } + } + TD.Init(lpz); + AccumulateLaplacianInfo(m,TD); + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + { + if(!SmoothSelected || (*vi).IsS()) + { + CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P(); + (*vi).P() = (*vi).P() + Delta*mu ; + } + } + } // end for step } @@ -433,9 +430,9 @@ static void VertexCoordLaplacianQuality(MeshType &m, int step, bool SmoothSelect class HCSmoothInfo { public: - CoordType dif; - CoordType sum; - int cnt; + CoordType dif; + CoordType sum; + int cnt; }; static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=false ) @@ -472,29 +469,29 @@ static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=fa for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD()) TD[*vi].sum/=(float)TD[*vi].cnt; - // Second Loop compute average difference - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - for(int j=0;j<3;++j) - { - TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P(); - TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P(); - // se l'edge j e' di bordo si deve sommare due volte - if((*fi).IsB(j)) - { - TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P(); - TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P(); - } - } - } + // Second Loop compute average difference + for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) + { + for(int j=0;j<3;++j) + { + TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P(); + TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P(); + // se l'edge j e' di bordo si deve sommare due volte + if((*fi).IsB(j)) + { + TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P(); + TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P(); + } + } + } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - { - TD[*vi].dif/=(float)TD[*vi].cnt; - if(!SmoothSelected || (*vi).IsS()) - (*vi).P()= TD[*vi].sum - (TD[*vi].sum - (*vi).P())*beta + (TD[*vi].dif)*(1.f-beta); - } - } // end for step + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + { + TD[*vi].dif/=(float)TD[*vi].cnt; + if(!SmoothSelected || (*vi).IsS()) + (*vi).P()= TD[*vi].sum - (TD[*vi].sum - (*vi).P())*beta + (TD[*vi].dif)*(1.f-beta); + } + } // end for step }; // Laplacian smooth of the quality. @@ -503,124 +500,124 @@ static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=fa class ColorSmoothInfo { public: - unsigned int r; - unsigned int g; - unsigned int b; - unsigned int a; - int cnt; + unsigned int r; + unsigned int g; + unsigned int b; + unsigned int a; + int cnt; }; static void VertexColorLaplacian(MeshType &m, int step, bool SmoothSelected=false, vcg::CallBackPos * cb=0) { - ColorSmoothInfo csi; - csi.r=0; csi.g=0; csi.b=0; csi.cnt=0; - SimpleTempData TD(m.vert,csi); + ColorSmoothInfo csi; + csi.r=0; csi.g=0; csi.b=0; csi.cnt=0; + SimpleTempData TD(m.vert,csi); - for(int i=0;iC()[0]; - TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1]; - TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2]; - TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3]; + FaceIterator fi; + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) + { + TD[(*fi).V(j)].r+=(*fi).V1(j)->C()[0]; + TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1]; + TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2]; + TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3]; - TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0]; - TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1]; - TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2]; - TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3]; + TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0]; + TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1]; + TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2]; + TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3]; - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - // si azzaera i dati per i vertici di bordo - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)]=csi; - TD[(*fi).V1(j)]=csi; - } + // si azzaera i dati per i vertici di bordo + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)]=csi; + TD[(*fi).V1(j)]=csi; + } - // se l'edge j e' di bordo si deve mediare solo con gli adiacenti - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)].r+=(*fi).V1(j)->C()[0]; - TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1]; - TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2]; - TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3]; + // se l'edge j e' di bordo si deve mediare solo con gli adiacenti + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)].r+=(*fi).V1(j)->C()[0]; + TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1]; + TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2]; + TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3]; - TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0]; - TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1]; - TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2]; - TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3]; + TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0]; + TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1]; + TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2]; + TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3]; - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - if(!SmoothSelected || (*vi).IsS()) - { - (*vi).C()[0] = (unsigned int) ceil((double) (TD[*vi].r / TD[*vi].cnt)); - (*vi).C()[1] = (unsigned int) ceil((double) (TD[*vi].g / TD[*vi].cnt)); - (*vi).C()[2] = (unsigned int) ceil((double) (TD[*vi].b / TD[*vi].cnt)); - (*vi).C()[3] = (unsigned int) ceil((double) (TD[*vi].a / TD[*vi].cnt)); - } - } // end for step + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + if(!SmoothSelected || (*vi).IsS()) + { + (*vi).C()[0] = (unsigned int) ceil((double) (TD[*vi].r / TD[*vi].cnt)); + (*vi).C()[1] = (unsigned int) ceil((double) (TD[*vi].g / TD[*vi].cnt)); + (*vi).C()[2] = (unsigned int) ceil((double) (TD[*vi].b / TD[*vi].cnt)); + (*vi).C()[3] = (unsigned int) ceil((double) (TD[*vi].a / TD[*vi].cnt)); + } + } // end for step }; static void FaceColorLaplacian(MeshType &m, int step, bool SmoothSelected=false, vcg::CallBackPos * cb=0) { - ColorSmoothInfo csi; - csi.r=0; csi.g=0; csi.b=0; csi.cnt=0; - SimpleTempData TD(m.face,csi); + ColorSmoothInfo csi; + csi.r=0; csi.g=0; csi.b=0; csi.cnt=0; + SimpleTempData TD(m.face,csi); - for(int i=0;iC()[0]; - TD[*fi].g+=(*fi).FFp(j)->C()[1]; - TD[*fi].b+=(*fi).FFp(j)->C()[2]; - TD[*fi].a+=(*fi).FFp(j)->C()[3]; - ++TD[*fi].cnt; - } - } - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD() && TD[*fi].cnt>0 ) - if(!SmoothSelected || (*fi).IsS()) - { - (*fi).C()[0] = (unsigned int) ceil((float) (TD[*fi].r / TD[*fi].cnt)); - (*fi).C()[1] = (unsigned int) ceil((float) (TD[*fi].g / TD[*fi].cnt)); - (*fi).C()[2] = (unsigned int) ceil((float) (TD[*fi].b / TD[*fi].cnt)); - (*fi).C()[3] = (unsigned int) ceil((float) (TD[*fi].a / TD[*fi].cnt)); - } - } // end for step + for(fi=m.face.begin();fi!=m.face.end();++fi) + { + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) + { + TD[*fi].r+=(*fi).FFp(j)->C()[0]; + TD[*fi].g+=(*fi).FFp(j)->C()[1]; + TD[*fi].b+=(*fi).FFp(j)->C()[2]; + TD[*fi].a+=(*fi).FFp(j)->C()[3]; + ++TD[*fi].cnt; + } + } + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD() && TD[*fi].cnt>0 ) + if(!SmoothSelected || (*fi).IsS()) + { + (*fi).C()[0] = (unsigned int) ceil((float) (TD[*fi].r / TD[*fi].cnt)); + (*fi).C()[1] = (unsigned int) ceil((float) (TD[*fi].g / TD[*fi].cnt)); + (*fi).C()[2] = (unsigned int) ceil((float) (TD[*fi].b / TD[*fi].cnt)); + (*fi).C()[3] = (unsigned int) ceil((float) (TD[*fi].a / TD[*fi].cnt)); + } + } // end for step }; // Laplacian smooth of the quality. @@ -628,8 +625,8 @@ static void FaceColorLaplacian(MeshType &m, int step, bool SmoothSelected=false, class QualitySmoothInfo { public: - ScalarType sum; - int cnt; + ScalarType sum; + int cnt; }; @@ -646,48 +643,48 @@ static void VertexQualityLaplacian(MeshType &m, int step=1, bool SmoothSelected= for(vi=m.vert.begin();vi!=m.vert.end();++vi) TD[*vi]=lpz; - FaceIterator fi; - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if(!(*fi).IsB(j)) - { - TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + FaceIterator fi; + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - // si azzaera i dati per i vertici di bordo - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)]=lpz; - TD[(*fi).V1(j)]=lpz; - } + // si azzaera i dati per i vertici di bordo + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)]=lpz; + TD[(*fi).V1(j)]=lpz; + } - // se l'edge j e' di bordo si deve mediare solo con gli adiacenti - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + // se l'edge j e' di bordo si deve mediare solo con gli adiacenti + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - //VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - if(!SmoothSelected || (*vi).IsS()) - (*vi).Q()=TD[*vi].sum/TD[*vi].cnt; - } + //VertexIterator vi; + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + if(!SmoothSelected || (*vi).IsS()) + (*vi).Q()=TD[*vi].sum/TD[*vi].cnt; + } - //TD.Stop(); + //TD.Stop(); }; static void VertexNormalLaplacian(MeshType &m, int step,bool SmoothSelected=false) @@ -702,46 +699,46 @@ static void VertexNormalLaplacian(MeshType &m, int step,bool SmoothSelected=fals for(vi=m.vert.begin();vi!=m.vert.end();++vi) TD[*vi]=lpz; - FaceIterator fi; - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if(!(*fi).IsB(j)) - { - TD[(*fi).V(j)].sum+=(*fi).V1(j)->N(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->N(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + FaceIterator fi; + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->N(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->N(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - // si azzaera i dati per i vertici di bordo - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)]=lpz; - TD[(*fi).V1(j)]=lpz; - } + // si azzaera i dati per i vertici di bordo + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)]=lpz; + TD[(*fi).V1(j)]=lpz; + } - // se l'edge j e' di bordo si deve mediare solo con gli adiacenti - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)].sum+=(*fi).V1(j)->N(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->N(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + // se l'edge j e' di bordo si deve mediare solo con gli adiacenti + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->N(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->N(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - //VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - if(!SmoothSelected || (*vi).IsS()) - (*vi).N()=TD[*vi].sum/TD[*vi].cnt; - } + //VertexIterator vi; + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + if(!SmoothSelected || (*vi).IsS()) + (*vi).N()=TD[*vi].sum/TD[*vi].cnt; + } }; @@ -749,41 +746,41 @@ static void VertexNormalLaplacian(MeshType &m, int step,bool SmoothSelected=fals // alpha e' compreso fra 0(no smoot) e 1 (tutto smoot) // Nota che se smootare il bordo puo far fare bandierine. static void VertexCoordViewDepth(MeshType &m, - const CoordType & viewpoint, - const ScalarType alpha, - int step, bool SmoothBorder=false ) + const CoordType & viewpoint, + const ScalarType alpha, + int step, bool SmoothBorder=false ) { - LaplacianInfo lpz; - lpz.sum=CoordType(0,0,0); - lpz.cnt=0; - SimpleTempData TD(m.vert,lpz); - for(int i=0;i TD(m.vert,lpz); + for(int i=0;icP(); - TD[(*fi).V1(j)].sum+=(*fi).V(j)->cP(); - ++TD[(*fi).V(j)].cnt; - ++TD[(*fi).V1(j)].cnt; - } + FaceIterator fi; + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if(!(*fi).IsB(j)) + { + TD[(*fi).V(j)].sum+=(*fi).V1(j)->cP(); + TD[(*fi).V1(j)].sum+=(*fi).V(j)->cP(); + ++TD[(*fi).V(j)].cnt; + ++TD[(*fi).V1(j)].cnt; + } - // si azzaera i dati per i vertici di bordo - for(fi=m.face.begin();fi!=m.face.end();++fi) - if(!(*fi).IsD()) - for(int j=0;j<3;++j) - if((*fi).IsB(j)) - { - TD[(*fi).V(j)]=lpz; - TD[(*fi).V1(j)]=lpz; - } + // si azzaera i dati per i vertici di bordo + for(fi=m.face.begin();fi!=m.face.end();++fi) + if(!(*fi).IsD()) + for(int j=0;j<3;++j) + if((*fi).IsB(j)) + { + TD[(*fi).V(j)]=lpz; + TD[(*fi).V1(j)]=lpz; + } // se l'edge j e' di bordo si deve mediare solo con gli adiacenti if(SmoothBorder) @@ -798,15 +795,15 @@ static void VertexCoordViewDepth(MeshType &m, ++TD[(*fi).V1(j)].cnt; } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if(!(*vi).IsD() && TD[*vi].cnt>0 ) - { - CoordType np = TD[*vi].sum/TD[*vi].cnt; - CoordType d = (*vi).cP() - viewpoint; d.Normalize(); - ScalarType s = d .dot ( np - (*vi).cP() ); - (*vi).P() += d * (s*alpha); - } - } + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if(!(*vi).IsD() && TD[*vi].cnt>0 ) + { + CoordType np = TD[*vi].sum/TD[*vi].cnt; + CoordType d = (*vi).cP() - viewpoint; d.Normalize(); + ScalarType s = d .dot ( np - (*vi).cP() ); + (*vi).P() += d * (s*alpha); + } + } } @@ -827,14 +824,14 @@ static void VertexCoordViewDepth(MeshType &m, class PDVertInfo { public: - CoordType np; + CoordType np; }; class PDFaceInfo { public: - CoordType m; + CoordType m; }; /***************************************************************************/ // Paso Doble Step 1 compute the smoothed normals @@ -848,26 +845,26 @@ public: void FaceNormalFuzzyVectorSB(MeshType &m, - SimpleTempData &TD, - ScalarType sigma) + SimpleTempData &TD, + ScalarType sigma) { - int i; + int i; - FaceIterator fi; + FaceIterator fi; - for(fi=m.face.begin();fi!=m.face.end();++fi) - { - CoordType bc=(*fi).Barycenter(); - // 1) Clear all the visited flag of faces that are vertex-adjacent to fi - for(i=0;i<3;++i) - { - vcg::face::VFIterator ep(&*fi,i); - while (!ep.End()) - { - ep.f->ClearV(); - ++ep; - } - } + for(fi=m.face.begin();fi!=m.face.end();++fi) + { + CoordType bc=(*fi).Barycenter(); + // 1) Clear all the visited flag of faces that are vertex-adjacent to fi + for(i=0;i<3;++i) + { + vcg::face::VFIterator ep(&*fi,i); + while (!ep.End()) + { + ep.f->ClearV(); + ++ep; + } + } // 1) Effectively average the normals weighting them with (*fi).SetV(); @@ -904,54 +901,54 @@ void FaceNormalFuzzyVectorSB(MeshType &m, static void FaceNormalLaplacianVF(MeshType &m) { - SimpleTempData TDF(m.face); + SimpleTempData TDF(m.face); - PDFaceInfo lpzf; - lpzf.m=CoordType(0,0,0); + PDFaceInfo lpzf; + lpzf.m=CoordType(0,0,0); - assert(tri::HasPerVertexVFAdjacency(m) && tri::HasPerFaceVFAdjacency(m) ); - TDF.Start(lpzf); - int i; + assert(tri::HasPerVertexVFAdjacency(m) && tri::HasPerFaceVFAdjacency(m) ); + TDF.Start(lpzf); + int i; FaceIterator fi; tri::UpdateNormal::AreaNormalizeFace(m); - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - CoordType bc=Barycenter(*fi); - // 1) Clear all the visited flag of faces that are vertex-adjacent to fi - for(i=0;i<3;++i) - { - VFLocalIterator ep(&*fi,i); - for (;!ep.End();++ep) - ep.f->ClearV(); - } + for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) + { + CoordType bc=Barycenter(*fi); + // 1) Clear all the visited flag of faces that are vertex-adjacent to fi + for(i=0;i<3;++i) + { + VFLocalIterator ep(&*fi,i); + for (;!ep.End();++ep) + ep.f->ClearV(); + } - // 2) Effectively average the normals - CoordType normalSum=(*fi).N(); + // 2) Effectively average the normals + CoordType normalSum=(*fi).N(); - for(i=0;i<3;++i) - { - VFLocalIterator ep(&*fi,i); - for (;!ep.End();++ep) - { - if(! (*ep.f).IsV() ) - { - normalSum += ep.f->N(); - (*ep.f).SetV(); - } - } - } - normalSum.Normalize(); - TDF[*fi].m=normalSum; - } - for(fi=m.face.begin();fi!=m.face.end();++fi) - (*fi).N()=TDF[*fi].m; + for(i=0;i<3;++i) + { + VFLocalIterator ep(&*fi,i); + for (;!ep.End();++ep) + { + if(! (*ep.f).IsV() ) + { + normalSum += ep.f->N(); + (*ep.f).SetV(); + } + } + } + normalSum.Normalize(); + TDF[*fi].m=normalSum; + } + for(fi=m.face.begin();fi!=m.face.end();++fi) + (*fi).N()=TDF[*fi].m; - tri::UpdateNormal::NormalizePerFace(m); + tri::UpdateNormal::NormalizePerFace(m); - TDF.Stop(); + TDF.Stop(); } // Replace the normal of the face with the average of normals of the face adijacent faces. @@ -1005,55 +1002,55 @@ static void FaceNormalLaplacianFF(MeshType &m, int step=1, bool SmoothSelected=f static void FaceNormalAngleThreshold(MeshType &m, - SimpleTempData &TD, - ScalarType sigma) + SimpleTempData &TD, + ScalarType sigma) { - int i; + int i; FaceIterator fi; - for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) - { - CoordType bc=Barycenter(*fi); - // 1) Clear all the visited flag of faces that are vertex-adjacent to fi - for(i=0;i<3;++i) - { - VFLocalIterator ep(&*fi,i); - for (;!ep.End();++ep) - ep.f->ClearV(); - } + for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD()) + { + CoordType bc=Barycenter(*fi); + // 1) Clear all the visited flag of faces that are vertex-adjacent to fi + for(i=0;i<3;++i) + { + VFLocalIterator ep(&*fi,i); + for (;!ep.End();++ep) + ep.f->ClearV(); + } - // 1) Effectively average the normals weighting them with the squared difference of the angle similarity - // sigma is the cosine of a threshold angle. sigma \in 0..1 - // sigma == 0 All the normals are averaged - // sigma == 1 Nothing is averaged. - // The averaging is weighted with the difference between normals. more similar the normal more important they are. + // 1) Effectively average the normals weighting them with the squared difference of the angle similarity + // sigma is the cosine of a threshold angle. sigma \in 0..1 + // sigma == 0 All the normals are averaged + // sigma == 1 Nothing is averaged. + // The averaging is weighted with the difference between normals. more similar the normal more important they are. - CoordType normalSum=CoordType(0,0,0); - for(i=0;i<3;++i) - { - VFLocalIterator ep(&*fi,i); - for (;!ep.End();++ep) - { - if(! (*ep.f).IsV() ) - { - ScalarType cosang=ep.f->N().dot((*fi).N()); - // Note that if two faces form an angle larger than 90 deg, their contribution should be very very small. - // Without this clamping - cosang = math::Clamp(cosang,0.0001f,1.f); - if(cosang >= sigma) - { - ScalarType w = cosang-sigma; - normalSum += ep.f->N()*(w*w); // similar normals have a cosang very close to 1 so cosang - sigma is maximized - } - (*ep.f).SetV(); - } - } - } - normalSum.Normalize(); - TD[*fi].m=normalSum; - } + CoordType normalSum=CoordType(0,0,0); + for(i=0;i<3;++i) + { + VFLocalIterator ep(&*fi,i); + for (;!ep.End();++ep) + { + if(! (*ep.f).IsV() ) + { + ScalarType cosang=ep.f->N().dot((*fi).N()); + // Note that if two faces form an angle larger than 90 deg, their contribution should be very very small. + // Without this clamping + cosang = math::Clamp(cosang,0.0001f,1.f); + if(cosang >= sigma) + { + ScalarType w = cosang-sigma; + normalSum += ep.f->N()*(w*w); // similar normals have a cosang very close to 1 so cosang - sigma is maximized + } + (*ep.f).SetV(); + } + } + } + normalSum.Normalize(); + TD[*fi].m=normalSum; + } for(fi=m.face.begin();fi!=m.face.end();++fi) (*fi).N()=TD[*fi].m; @@ -1066,35 +1063,35 @@ static void FaceNormalAngleThreshold(MeshType &m, static CoordType TriAreaGradient(CoordType &p,CoordType &p0,CoordType &p1) { - CoordType dd = p1-p0; - CoordType d0 = p-p0; - CoordType d1 = p-p1; - CoordType grad; + CoordType dd = p1-p0; + CoordType d0 = p-p0; + CoordType d1 = p-p1; + CoordType grad; - ScalarType t16 = d0[1]* d1[2] - d0[2]* d1[1]; - ScalarType t5 = -d0[2]* d1[0] + d0[0]* d1[2]; - ScalarType t4 = -d0[0]* d1[1] + d0[1]* d1[0]; + ScalarType t16 = d0[1]* d1[2] - d0[2]* d1[1]; + ScalarType t5 = -d0[2]* d1[0] + d0[0]* d1[2]; + ScalarType t4 = -d0[0]* d1[1] + d0[1]* d1[0]; - ScalarType delta= sqrtf(t4*t4 + t5*t5 +t16*t16); + ScalarType delta= sqrtf(t4*t4 + t5*t5 +t16*t16); - grad[0]= (t5 * (-dd[2]) + t4 * ( dd[1]))/delta; - grad[1]= (t16 * (-dd[2]) + t4 * (-dd[0]))/delta; - grad[2]= (t16 * ( dd[1]) + t5 * ( dd[0]))/delta; + grad[0]= (t5 * (-dd[2]) + t4 * ( dd[1]))/delta; + grad[1]= (t16 * (-dd[2]) + t4 * (-dd[0]))/delta; + grad[2]= (t16 * ( dd[1]) + t5 * ( dd[0]))/delta; - return grad; + return grad; } template static CoordType CrossProdGradient(CoordType &p, CoordType &p0, CoordType &p1, CoordType &m) { - CoordType grad; - CoordType p00=p0-p; - CoordType p01=p1-p; - grad[0] = (-p00[2] + p01[2])*m[1] + (-p01[1] + p00[1])*m[2]; - grad[1] = (-p01[2] + p00[2])*m[0] + (-p00[0] + p01[0])*m[2]; - grad[2] = (-p00[1] + p01[1])*m[0] + (-p01[0] + p00[0])*m[1]; + CoordType grad; + CoordType p00=p0-p; + CoordType p01=p1-p; + grad[0] = (-p00[2] + p01[2])*m[1] + (-p01[1] + p00[1])*m[2]; + grad[1] = (-p01[2] + p00[2])*m[0] + (-p00[0] + p01[0])*m[2]; + grad[2] = (-p00[1] + p01[1])*m[0] + (-p01[0] + p00[0])*m[1]; - return grad; + return grad; } /* @@ -1110,8 +1107,8 @@ A(...) (2-2nm) = static CoordType FaceErrorGrad(CoordType &p,CoordType &p0,CoordType &p1, CoordType &m) { - return TriAreaGradient(p,p0,p1) *2.0f - - CrossProdGradient(p,p0,p1,m) *2.0f ; + return TriAreaGradient(p,p0,p1) *2.0f + - CrossProdGradient(p,p0,p1,m) *2.0f ; } /***************************************************************************/ // Paso Doble Step 2 Fitta la mesh a un dato insieme di normali @@ -1119,29 +1116,29 @@ static CoordType FaceErrorGrad(CoordType &p,CoordType &p0,CoordType &p1, CoordTy static void FitMesh(MeshType &m, - SimpleTempData &TDV, - SimpleTempData &TDF, - float lambda) + SimpleTempData &TDV, + SimpleTempData &TDF, + float lambda) { - //vcg::face::Pos ep; - vcg::face::VFIterator ep; - VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - { - CoordType ErrGrad=CoordType(0,0,0); + //vcg::face::Pos ep; + vcg::face::VFIterator ep; + VertexIterator vi; + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + { + CoordType ErrGrad=CoordType(0,0,0); - ep.f=(*vi).VFp(); - ep.z=(*vi).VFi(); - while (!ep.End()) - { - ErrGrad+=FaceErrorGrad(ep.f->V(ep.z)->P(),ep.f->V1(ep.z)->P(),ep.f->V2(ep.z)->P(),TDF[ep.f].m); - ++ep; - } - TDV[*vi].np=(*vi).P()-ErrGrad*(ScalarType)lambda; - } + ep.f=(*vi).VFp(); + ep.z=(*vi).VFi(); + while (!ep.End()) + { + ErrGrad+=FaceErrorGrad(ep.f->V(ep.z)->P(),ep.f->V1(ep.z)->P(),ep.f->V2(ep.z)->P(),TDF[ep.f].m); + ++ep; + } + TDV[*vi].np=(*vi).P()-ErrGrad*(ScalarType)lambda; + } - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - (*vi).P()=TDV[*vi].np; + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + (*vi).P()=TDV[*vi].np; } /****************************************************************************************************************/ @@ -1153,12 +1150,12 @@ static void FastFitMesh(MeshType &m, //SimpleTempData &TDF, bool OnlySelected=false) { - //vcg::face::Pos ep; - vcg::face::VFIterator ep; - VertexIterator vi; + //vcg::face::Pos ep; + vcg::face::VFIterator ep; + VertexIterator vi; - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - { + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + { CoordType Sum(0,0,0); ScalarType cnt=0; VFLocalIterator ep(&*vi); @@ -1171,44 +1168,44 @@ static void FastFitMesh(MeshType &m, TDV[*vi].np=(*vi).P()+ Sum*(1.0/cnt); } - if(OnlySelected) - { - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - if((*vi).IsS()) (*vi).P()=TDV[*vi].np; - } - else - { - for(vi=m.vert.begin();vi!=m.vert.end();++vi) - (*vi).P()=TDV[*vi].np; - } + if(OnlySelected) + { + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + if((*vi).IsS()) (*vi).P()=TDV[*vi].np; + } + else + { + for(vi=m.vert.begin();vi!=m.vert.end();++vi) + (*vi).P()=TDV[*vi].np; + } } static void VertexCoordPasoDoble(MeshType &m, int step, typename MeshType::ScalarType Sigma=0, int FitStep=10, typename MeshType::ScalarType FitLambda=0.05) { - SimpleTempData< typename MeshType::VertContainer, PDVertInfo> TDV(m.vert); - SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face); - PDVertInfo lpzv; - lpzv.np=CoordType(0,0,0); - PDFaceInfo lpzf; - lpzf.m=CoordType(0,0,0); + SimpleTempData< typename MeshType::VertContainer, PDVertInfo> TDV(m.vert); + SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face); + PDVertInfo lpzv; + lpzv.np=CoordType(0,0,0); + PDFaceInfo lpzf; + lpzf.m=CoordType(0,0,0); - assert(m.HasVFTopology()); - m.HasVFTopology(); - TDV.Start(lpzv); - TDF.Start(lpzf); - for(int j=0;j::PerFace(m); - FaceNormalAngleThreshold(m,TDF,Sigma); - for(int k=0;k::PerFace(m); + FaceNormalAngleThreshold(m,TDF,Sigma); + for(int k=0;k TDV(m.vert,lpzv); - SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face,lpzf); + assert(HasPerVertexVFAdjacency(m) && HasPerFaceVFAdjacency(m)); + SimpleTempData< typename MeshType::VertContainer, PDVertInfo> TDV(m.vert,lpzv); + SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face,lpzf); for(int j=0;j #include #include -#include namespace vcg { namespace tri { @@ -52,16 +51,16 @@ class UpdateCurvature { public: - typedef typename MeshType::FaceType FaceType; - typedef typename MeshType::FacePointer FacePointer; - typedef typename MeshType::FaceIterator FaceIterator; - typedef typename MeshType::VertexIterator VertexIterator; - typedef typename MeshType::VertContainer VertContainer; - typedef typename MeshType::VertexType VertexType; - typedef typename MeshType::VertexPointer VertexPointer; - typedef vcg::face::VFIterator VFIteratorType; - typedef typename MeshType::CoordType CoordType; - typedef typename CoordType::ScalarType ScalarType; + typedef typename MeshType::FaceType FaceType; + typedef typename MeshType::FacePointer FacePointer; + typedef typename MeshType::FaceIterator FaceIterator; + typedef typename MeshType::VertexIterator VertexIterator; + typedef typename MeshType::VertContainer VertContainer; + typedef typename MeshType::VertexType VertexType; + typedef typename MeshType::VertexPointer VertexPointer; + typedef vcg::face::VFIterator VFIteratorType; + typedef typename MeshType::CoordType CoordType; + typedef typename CoordType::ScalarType ScalarType; private: @@ -74,11 +73,11 @@ private: public: - /// \brief Compute principal direction and magnitudo of curvature. + /// \brief Compute principal direction and magnitudo of curvature. /* - Compute principal direction and magniuto of curvature as describe in the paper: - @InProceedings{bb33922, + Compute principal direction and magniuto of curvature as describe in the paper: + @InProceedings{bb33922, author = "G. Taubin", title = "Estimating the Tensor of Curvature of a Surface from a Polyhedral Approximation", @@ -336,13 +335,13 @@ If pointVSfaceInt==false the covariance is computed by (analytic)integration ove // Jacobi(A, eigenvalues , eigenvectors, nrot); - Eigen::Matrix3d AA; - A.ToEigenMatrix(AA); - Eigen::SelfAdjointEigenSolver eig(AA); - Eigen::Vector3d c_val = eig.eigenvalues(); - Eigen::Matrix3d c_vec = eig.eigenvectors(); // eigenvector are stored as columns. - eigenvectors.FromEigenMatrix(c_vec); - eigenvalues.FromEigenVector(c_val); + Eigen::Matrix3d AA; + A.ToEigenMatrix(AA); + Eigen::SelfAdjointEigenSolver eig(AA); + Eigen::Vector3d c_val = eig.eigenvalues(); + Eigen::Matrix3d c_vec = eig.eigenvectors(); // eigenvector are stored as columns. + eigenvectors.FromEigenMatrix(c_vec); + eigenvalues.FromEigenVector(c_val); // EV.transposeInPlace(); // ev.FromEigenVector(c_val); @@ -520,174 +519,174 @@ static void MeanAndGaussian(MeshType & m) } - /// \brief Update the mean and the gaussian curvature of a vertex. + /// \brief Update the mean and the gaussian curvature of a vertex. - /** - The function uses the VF adiacency to walk around the vertex. - \return It will return the voronoi area around the vertex. If (norm == true) the mean and the gaussian curvature are normalized. - Based on the paper "Optimizing 3d triangulations using discrete curvature analysis" - */ + /** + The function uses the VF adiacency to walk around the vertex. + \return It will return the voronoi area around the vertex. If (norm == true) the mean and the gaussian curvature are normalized. + Based on the paper "Optimizing 3d triangulations using discrete curvature analysis" + */ - static float ComputeSingleVertexCurvature(VertexPointer v, bool norm = true) - { - VFIteratorType vfi(v); - float A = 0; + static float ComputeSingleVertexCurvature(VertexPointer v, bool norm = true) + { + VFIteratorType vfi(v); + float A = 0; - v->Kh() = 0; - v->Kg() = 2 * M_PI; + v->Kh() = 0; + v->Kg() = 2 * M_PI; - while (!vfi.End()) { - if (!vfi.F()->IsD()) { - FacePointer f = vfi.F(); - int i = vfi.I(); - VertexPointer v0 = f->V0(i), v1 = f->V1(i), v2 = f->V2(i); + while (!vfi.End()) { + if (!vfi.F()->IsD()) { + FacePointer f = vfi.F(); + int i = vfi.I(); + VertexPointer v0 = f->V0(i), v1 = f->V1(i), v2 = f->V2(i); - float ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P() )); - float ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P() )); - float ang2 = M_PI - ang0 - ang1; + float ang0 = math::Abs(Angle(v1->P() - v0->P(), v2->P() - v0->P() )); + float ang1 = math::Abs(Angle(v0->P() - v1->P(), v2->P() - v1->P() )); + float ang2 = M_PI - ang0 - ang1; - float s01 = SquaredDistance(v1->P(), v0->P()); - float s02 = SquaredDistance(v2->P(), v0->P()); + float s01 = SquaredDistance(v1->P(), v0->P()); + float s02 = SquaredDistance(v2->P(), v0->P()); - // voronoi cell of current vertex - if (ang0 >= M_PI/2) - A += (0.5f * DoubleArea(*f) - (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 ); - else if (ang1 >= M_PI/2) - A += (s01 * tan(ang0)) / 8.0; - else if (ang2 >= M_PI/2) - A += (s02 * tan(ang0)) / 8.0; - else // non obctuse triangle - A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0; + // voronoi cell of current vertex + if (ang0 >= M_PI/2) + A += (0.5f * DoubleArea(*f) - (s01 * tan(ang1) + s02 * tan(ang2)) / 8.0 ); + else if (ang1 >= M_PI/2) + A += (s01 * tan(ang0)) / 8.0; + else if (ang2 >= M_PI/2) + A += (s02 * tan(ang0)) / 8.0; + else // non obctuse triangle + A += ((s02 / tan(ang1)) + (s01 / tan(ang2))) / 8.0; - // gaussian curvature update - v->Kg() -= ang0; + // gaussian curvature update + v->Kg() -= ang0; - // mean curvature update - ang1 = math::Abs(Angle(f->N(), v1->N())); - ang2 = math::Abs(Angle(f->N(), v2->N())); - v->Kh() += ( (math::Sqrt(s01) / 2.0) * ang1 + - (math::Sqrt(s02) / 2.0) * ang2 ); - } + // mean curvature update + ang1 = math::Abs(Angle(f->N(), v1->N())); + ang2 = math::Abs(Angle(f->N(), v2->N())); + v->Kh() += ( (math::Sqrt(s01) / 2.0) * ang1 + + (math::Sqrt(s02) / 2.0) * ang2 ); + } - ++vfi; - } + ++vfi; + } - v->Kh() /= 4.0f; + v->Kh() /= 4.0f; - if(norm) { - if(A <= std::numeric_limits::epsilon()) { - v->Kh() = 0; - v->Kg() = 0; - } - else { - v->Kh() /= A; - v->Kg() /= A; - } - } + if(norm) { + if(A <= std::numeric_limits::epsilon()) { + v->Kh() = 0; + v->Kg() = 0; + } + else { + v->Kh() /= A; + v->Kg() /= A; + } + } - return A; - } + return A; + } - static void PerVertex(MeshType & m) - { - tri::RequireVFAdjacency(m); + static void PerVertex(MeshType & m) + { + tri::RequireVFAdjacency(m); - for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) - ComputeSingleVertexCurvature(&*vi,false); - } + for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) + ComputeSingleVertexCurvature(&*vi,false); + } /* - Compute principal curvature directions and value with normal cycle: - @inproceedings{CohMor03, - author = {Cohen-Steiner, David and Morvan, Jean-Marie }, - booktitle = {SCG '03: Proceedings of the nineteenth annual symposium on Computational geometry}, - title - {Restricted delaunay triangulations and normal cycle} - year = {2003} + Compute principal curvature directions and value with normal cycle: + @inproceedings{CohMor03, + author = {Cohen-Steiner, David and Morvan, Jean-Marie }, + booktitle = {SCG '03: Proceedings of the nineteenth annual symposium on Computational geometry}, + title - {Restricted delaunay triangulations and normal cycle} + year = {2003} } - */ + */ - static void PrincipalDirectionsNormalCycle(MeshType & m){ - tri::RequireVFAdjacency(m); - tri::RequireFFAdjacency(m); - tri::RequirePerFaceNormal(m); + static void PrincipalDirectionsNormalCycle(MeshType & m){ + tri::RequireVFAdjacency(m); + tri::RequireFFAdjacency(m); + tri::RequirePerFaceNormal(m); - typename MeshType::VertexIterator vi; + typename MeshType::VertexIterator vi; - for(vi = m.vert.begin(); vi != m.vert.end(); ++vi) - if(!((*vi).IsD())){ - vcg::Matrix33 m33;m33.SetZero(); - face::JumpingPos p((*vi).VFp(),&(*vi)); - p.FlipE(); - typename MeshType::VertexType * firstv = p.VFlip(); - assert(p.F()->V(p.VInd())==&(*vi)); + for(vi = m.vert.begin(); vi != m.vert.end(); ++vi) + if(!((*vi).IsD())){ + vcg::Matrix33 m33;m33.SetZero(); + face::JumpingPos p((*vi).VFp(),&(*vi)); + p.FlipE(); + typename MeshType::VertexType * firstv = p.VFlip(); + assert(p.F()->V(p.VInd())==&(*vi)); - do{ - if( p.F() != p.FFlip()){ - Point3 normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P(); - ScalarType edge_length = normalized_edge.Norm(); - normalized_edge/=edge_length; - Point3 n1 = p.F()->cN();n1.Normalize(); - Point3 n2 = p.FFlip()->cN();n2.Normalize(); - ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge); - n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0)); - ScalarType beta = math::Asin(n1n2); - m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0]; - m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0]; - m33[1][1] += beta*edge_length*normalized_edge[1]*normalized_edge[1]; - m33[0][2] += beta*edge_length*normalized_edge[2]*normalized_edge[0]; - m33[1][2] += beta*edge_length*normalized_edge[2]*normalized_edge[1]; - m33[2][2] += beta*edge_length*normalized_edge[2]*normalized_edge[2]; - } - p.NextFE(); - }while(firstv != p.VFlip()); + do{ + if( p.F() != p.FFlip()){ + Point3 normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P(); + ScalarType edge_length = normalized_edge.Norm(); + normalized_edge/=edge_length; + Point3 n1 = p.F()->cN();n1.Normalize(); + Point3 n2 = p.FFlip()->cN();n2.Normalize(); + ScalarType n1n2 = (n1 ^ n2).dot(normalized_edge); + n1n2 = std::max(std::min( ScalarType(1.0),n1n2),ScalarType(-1.0)); + ScalarType beta = math::Asin(n1n2); + m33[0][0] += beta*edge_length*normalized_edge[0]*normalized_edge[0]; + m33[0][1] += beta*edge_length*normalized_edge[1]*normalized_edge[0]; + m33[1][1] += beta*edge_length*normalized_edge[1]*normalized_edge[1]; + m33[0][2] += beta*edge_length*normalized_edge[2]*normalized_edge[0]; + m33[1][2] += beta*edge_length*normalized_edge[2]*normalized_edge[1]; + m33[2][2] += beta*edge_length*normalized_edge[2]*normalized_edge[2]; + } + p.NextFE(); + }while(firstv != p.VFlip()); - if(m33.Determinant()==0.0){ // degenerate case - (*vi).K1() = (*vi).K2() = 0.0; continue;} + if(m33.Determinant()==0.0){ // degenerate case + (*vi).K1() = (*vi).K2() = 0.0; continue;} - m33[1][0] = m33[0][1]; - m33[2][0] = m33[0][2]; - m33[2][1] = m33[1][2]; + m33[1][0] = m33[0][1]; + m33[2][0] = m33[0][2]; + m33[2][1] = m33[1][2]; - Eigen::Matrix3d it; - m33.ToEigenMatrix(it); - Eigen::SelfAdjointEigenSolver eig(it); - Eigen::Vector3d c_val = eig.eigenvalues(); - Eigen::Matrix3d c_vec = eig.eigenvectors(); + Eigen::Matrix3d it; + m33.ToEigenMatrix(it); + Eigen::SelfAdjointEigenSolver eig(it); + Eigen::Vector3d c_val = eig.eigenvalues(); + Eigen::Matrix3d c_vec = eig.eigenvectors(); - Point3 lambda; - Matrix33 vect; - vect.FromEigenMatrix(c_vec); - lambda.FromEigenVector(c_val); + Point3 lambda; + Matrix33 vect; + vect.FromEigenMatrix(c_vec); + lambda.FromEigenVector(c_val); - ScalarType bestNormal = 0; - int bestNormalIndex = -1; - for(int i = 0; i < 3; ++i) - { - float agreeWithNormal = fabs((*vi).N().Normalize().dot(vect.GetColumn(i))); - if( agreeWithNormal > bestNormal ) - { - bestNormal= agreeWithNormal; - bestNormalIndex = i; - } - } - int maxI = (bestNormalIndex+2)%3; - int minI = (bestNormalIndex+1)%3; - if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI); + ScalarType bestNormal = 0; + int bestNormalIndex = -1; + for(int i = 0; i < 3; ++i) + { + float agreeWithNormal = fabs((*vi).N().Normalize().dot(vect.GetColumn(i))); + if( agreeWithNormal > bestNormal ) + { + bestNormal= agreeWithNormal; + bestNormalIndex = i; + } + } + int maxI = (bestNormalIndex+2)%3; + int minI = (bestNormalIndex+1)%3; + if(fabs(lambda[maxI]) < fabs(lambda[minI])) std::swap(maxI,minI); - (*vi).PD1() = *(Point3*)(& vect[maxI][0]); - (*vi).PD2() = *(Point3*)(& vect[minI][0]); - (*vi).K1() = lambda[2]; - (*vi).K2() = lambda[1]; - } - } + (*vi).PD1() = *(Point3*)(& vect[maxI][0]); + (*vi).PD2() = *(Point3*)(& vect[minI][0]); + (*vi).K1() = lambda[2]; + (*vi).K2() = lambda[1]; + } + } - static void PerVertexBasicRadialCrossField(MeshType &m, float anisotropyRatio = 1.0 ) - { - tri::RequirePerVertexCurvatureDir(m); - CoordType c=m.bbox.Center(); - float maxRad = m.bbox.Diag()/2.0f; + static void PerVertexBasicRadialCrossField(MeshType &m, float anisotropyRatio = 1.0 ) + { + tri::RequirePerVertexCurvatureDir(m); + CoordType c=m.bbox.Center(); + float maxRad = m.bbox.Diag()/2.0f; for(int i=0;i -#include -#include -#include #ifdef __glut_h__ #include @@ -37,755 +34,754 @@ #include #include #include -#include namespace vcg { - /*! - * Given an object or an object pointer, return the reference to the object - */ - template - struct Dereferencer - { - static TYPE& Ref(TYPE &t) { return ( t); } - static TYPE& Ref(TYPE* &t) { return (*t); } - static const TYPE& Ref(const TYPE &t) { return ( t); } - static const TYPE& Ref(const TYPE* &t) { return (*t); } - }; - - - /*! - * Given a type, return the type - */ - template - class ReferenceType - { - public: - typedef T Type; - }; - - - /*! - * Given as type a pointer to type, return the type - */ - template - class ReferenceType - { - public: - typedef typename ReferenceType::Type Type; - }; + /*! + * Given an object or an object pointer, return the reference to the object + */ + template + struct Dereferencer + { + static TYPE& Ref(TYPE &t) { return ( t); } + static TYPE& Ref(TYPE* &t) { return (*t); } + static const TYPE& Ref(const TYPE &t) { return ( t); } + static const TYPE& Ref(const TYPE* &t) { return (*t); } + }; - - /*! - * The type of the octree voxels - */ - struct Voxel - { - Voxel() { count = begin = end = -1; } - - void SetRange(const int begin, const int end) - { - this->begin = begin; - this->end = end; - count = end-begin; - }; - - void AddRange(const Voxel *voxel) - { - assert(voxel->end>end); - - count += voxel->count; - end = voxel->end; - }; - - int begin; - int end; - int count; - }; - - - template < class OBJECT_TYPE, class SCALAR_TYPE> - class Octree : public vcg::OctreeTemplate< Voxel, SCALAR_TYPE >, public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE > - { - protected: - struct Neighbour; - - public: - typedef SCALAR_TYPE ScalarType; - typedef OBJECT_TYPE ObjectType; - typedef typename Octree::Leaf * LeafPointer; - typedef typename Octree::InnerNode * InnerNodePointer; - typedef typename ReferenceType::Type * ObjectPointer; - - typedef vcg::Voxel VoxelType; - typedef VoxelType * VoxelPointer; - - typedef vcg::OctreeTemplate< VoxelType, SCALAR_TYPE > TemplatedOctree; - typedef typename TemplatedOctree::ZOrderType ZOrderType; - - typedef typename TemplatedOctree::BoundingBoxType BoundingBoxType; - typedef typename TemplatedOctree::CenterType CenterType; - typedef typename TemplatedOctree::CoordinateType CoordType; - - typedef typename TemplatedOctree::NodeType NodeType; - typedef typename TemplatedOctree::NodePointer NodePointer; - typedef typename TemplatedOctree::NodeIndex NodeIndex; - - typedef typename std::vector< Neighbour >::iterator NeighbourIterator; - - protected: - /*********************************************** - * INNER DATA STRUCTURES AND PREDICATES * - ***********************************************/ - /*! - * Structure used during the sorting of the dataset - */ - template < typename LEAF_TYPE > - struct ObjectPlaceholder - { - typedef LEAF_TYPE* LeafPointer; - - ObjectPlaceholder() { z_order = object_index = -1, leaf_pointer = NULL;} - - ObjectPlaceholder(ZOrderType zOrder, void* leafPointer, unsigned int objectIndex) - { - z_order = zOrder; - leaf_pointer = leafPointer; - object_index = objectIndex; - } - - ZOrderType z_order; - LeafPointer leaf_pointer; - unsigned int object_index; - }; - - - /*! - * Predicate used during the sorting of the dataset - */ - template - struct ObjectSorter - { - inline bool operator()(const ObjectPlaceholder< LEAF_TYPE > &first, const ObjectPlaceholder< LEAF_TYPE > &second) - { - return (first.z_orderobject = NULL; - this->distance = -1.0f; - }; - - Neighbour(ObjectPointer &object, CoordType &point, ScalarType distance) - { - this->object = object; - this->point = point; - this->distance = distance; - } - - inline bool operator<(const Neighbour &n) const - { - return distance + class ReferenceType + { + public: + typedef T Type; + }; - ObjectPointer object; - CoordType point; - ScalarType distance; - }; + /*! + * Given as type a pointer to type, return the type + */ + template + class ReferenceType + { + public: + typedef typename ReferenceType::Type Type; + }; + + + + /*! + * The type of the octree voxels + */ + struct Voxel + { + Voxel() { count = begin = end = -1; } + + void SetRange(const int begin, const int end) + { + this->begin = begin; + this->end = end; + count = end-begin; + }; + + void AddRange(const Voxel *voxel) + { + assert(voxel->end>end); + + count += voxel->count; + end = voxel->end; + }; + + int begin; + int end; + int count; + }; + + + template < class OBJECT_TYPE, class SCALAR_TYPE> + class Octree : public vcg::OctreeTemplate< Voxel, SCALAR_TYPE >, public vcg::SpatialIndex< OBJECT_TYPE, SCALAR_TYPE > + { + protected: + struct Neighbour; + + public: + typedef SCALAR_TYPE ScalarType; + typedef OBJECT_TYPE ObjectType; + typedef typename Octree::Leaf * LeafPointer; + typedef typename Octree::InnerNode * InnerNodePointer; + typedef typename ReferenceType::Type * ObjectPointer; + + typedef vcg::Voxel VoxelType; + typedef VoxelType * VoxelPointer; + + typedef vcg::OctreeTemplate< VoxelType, SCALAR_TYPE > TemplatedOctree; + typedef typename TemplatedOctree::ZOrderType ZOrderType; + + typedef typename TemplatedOctree::BoundingBoxType BoundingBoxType; + typedef typename TemplatedOctree::CenterType CenterType; + typedef typename TemplatedOctree::CoordinateType CoordType; + + typedef typename TemplatedOctree::NodeType NodeType; + typedef typename TemplatedOctree::NodePointer NodePointer; + typedef typename TemplatedOctree::NodeIndex NodeIndex; + + typedef typename std::vector< Neighbour >::iterator NeighbourIterator; + + protected: + /*********************************************** + * INNER DATA STRUCTURES AND PREDICATES * + ***********************************************/ + /*! + * Structure used during the sorting of the dataset + */ + template < typename LEAF_TYPE > + struct ObjectPlaceholder + { + typedef LEAF_TYPE* LeafPointer; + + ObjectPlaceholder() { z_order = object_index = -1, leaf_pointer = NULL;} + + ObjectPlaceholder(ZOrderType zOrder, void* leafPointer, unsigned int objectIndex) + { + z_order = zOrder; + leaf_pointer = leafPointer; + object_index = objectIndex; + } + + ZOrderType z_order; + LeafPointer leaf_pointer; + unsigned int object_index; + }; + + + /*! + * Predicate used during the sorting of the dataset + */ + template + struct ObjectSorter + { + inline bool operator()(const ObjectPlaceholder< LEAF_TYPE > &first, const ObjectPlaceholder< LEAF_TYPE > &second) + { + return (first.z_orderobject = NULL; + this->distance = -1.0f; + }; + + Neighbour(ObjectPointer &object, CoordType &point, ScalarType distance) + { + this->object = object; + this->point = point; + this->distance = distance; + } + + inline bool operator<(const Neighbour &n) const + { + return distance - void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj /*, vcg::CallBackPos *callback=NULL*/) - { - // Compute the bounding-box enclosing the whole dataset - typedef Dereferencer::Type > DereferencerType; - BoundingBoxType bounding_box, obj_bb; - bounding_box.SetNull(); - for (OBJECT_ITERATOR iObj=bObj; iObj!=eObj; iObj++) - { - (*iObj).GetBBox(obj_bb); - bounding_box.Add(obj_bb); - } + /*! + * Populate the octree + */ + template < class OBJECT_ITERATOR > + void Set(const OBJECT_ITERATOR & bObj, const OBJECT_ITERATOR & eObj /*, vcg::CallBackPos *callback=NULL*/) + { + // Compute the bounding-box enclosing the whole dataset + typedef Dereferencer::Type > DereferencerType; + BoundingBoxType bounding_box, obj_bb; + bounding_box.SetNull(); + for (OBJECT_ITERATOR iObj=bObj; iObj!=eObj; iObj++) + { + (*iObj).GetBBox(obj_bb); + bounding_box.Add(obj_bb); + } - //...and expand it a bit more - BoundingBoxType resulting_bb(bounding_box); - CoordType offset = bounding_box.Dim()*Octree::EXPANSION_FACTOR; - CoordType center = bounding_box.Center(); - resulting_bb.Offset(offset); + //...and expand it a bit more + BoundingBoxType resulting_bb(bounding_box); + CoordType offset = bounding_box.Dim()*Octree::EXPANSION_FACTOR; + CoordType center = bounding_box.Center(); + resulting_bb.Offset(offset); ScalarType longest_side = vcg::math::Max( resulting_bb.DimX(), resulting_bb.DimY(), resulting_bb.DimZ())/2.0f; - resulting_bb.Set(center); - resulting_bb.Offset(longest_side); - TemplatedOctree::boundingBox = resulting_bb; + resulting_bb.Set(center); + resulting_bb.Offset(longest_side); + TemplatedOctree::boundingBox = resulting_bb; - // Try to find a reasonable octree depth - int dataset_dimension = int(std::distance(bObj, eObj)); + // Try to find a reasonable octree depth + int dataset_dimension = int(std::distance(bObj, eObj)); - int primitives_per_voxel; - int depth = 4; - do - { - int number_of_voxel = 1<<(3*depth); // i.e. 8^depth - float density = float(number_of_voxel)/float(depth); - primitives_per_voxel = int(float(dataset_dimension)/density); - depth++; - } - while (primitives_per_voxel>25 && depth<15); - TemplatedOctree::Initialize(++depth); + int primitives_per_voxel; + int depth = 4; + do + { + int number_of_voxel = 1<<(3*depth); // i.e. 8^depth + float density = float(number_of_voxel)/float(depth); + primitives_per_voxel = int(float(dataset_dimension)/density); + depth++; + } + while (primitives_per_voxel>25 && depth<15); + TemplatedOctree::Initialize(++depth); - // Sort the dataset (using the lebesgue space filling curve...) - std::string message("Indexing dataset..."); - NodePointer *route = new NodePointer[depth+1]; - OBJECT_ITERATOR iObj = bObj; + // Sort the dataset (using the lebesgue space filling curve...) + std::string message("Indexing dataset..."); + NodePointer *route = new NodePointer[depth+1]; + OBJECT_ITERATOR iObj = bObj; - //if (callback!=NULL) callback(int((i+1)*100/dataset_dimension), message.c_str()); + //if (callback!=NULL) callback(int((i+1)*100/dataset_dimension), message.c_str()); - std::vector< ObjectPlaceholder< NodeType > > placeholders/*(dataset_dimension)*/; - vcg::Box3 object_bb; - vcg::Point3 hit_leaf; - for (int i=0; i > placeholders/*(dataset_dimension)*/; + vcg::Box3 object_bb; + vcg::Point3 hit_leaf; + for (int i=0; i() ); + while (object_bb.IsIn(hit_leaf)) + { + int placeholder_index = int(placeholders.size()); + placeholders.push_back( ObjectPlaceholder< NodeType >() ); placeholders[placeholder_index].z_order = TemplatedOctree::BuildRoute(hit_leaf, route); - placeholders[placeholder_index].leaf_pointer = route[depth]; - placeholders[placeholder_index].object_index = i; - - hit_leaf.X() += TemplatedOctree::leafDimension.X(); - if (hit_leaf.X()>object_bb.max.X()) - { - hit_leaf.X() = object_bb.min.X(); - hit_leaf.Z()+= TemplatedOctree::leafDimension.Z(); - if (hit_leaf.Z()>object_bb.max.Z()) - { - hit_leaf.Z() = object_bb.min.Z(); - hit_leaf.Y()+= TemplatedOctree::leafDimension.Y(); - } - } - } - } - delete []route; + placeholders[placeholder_index].leaf_pointer = route[depth]; + placeholders[placeholder_index].object_index = i; - int placeholder_count = int(placeholders.size()); - - // Allocate the mark array - global_mark = 1; - marks = new unsigned char[placeholder_count]; - memset(&marks[0], 0, sizeof(unsigned char)*placeholder_count); - - std::sort(placeholders.begin(), placeholders.end(), ObjectSorter< NodeType >()); - std::vector< NodePointer > filled_leaves(placeholder_count); - sorted_dataset.resize( placeholder_count ); - for (int i=0; iobject_bb.max.X()) + { + hit_leaf.X() = object_bb.min.X(); + hit_leaf.Z()+= TemplatedOctree::leafDimension.Z(); + if (hit_leaf.Z()>object_bb.max.Z()) + { + hit_leaf.Z() = object_bb.min.Z(); + hit_leaf.Y()+= TemplatedOctree::leafDimension.Y(); + } + } + } + } + delete []route; - // The dataset is sorted and the octree is built, but the indexing information aren't stored yet in the octree: - // we assign to each leaf the range inside the sorted dataset of the primitives contained inside the leaf - int begin = -1; - NodePointer initial_leaf = NULL; - for (int end=0; endSetRange(begin, end); - } + // Allocate the mark array + global_mark = 1; + marks = new unsigned char[placeholder_count]; + memset(&marks[0], 0, sizeof(unsigned char)*placeholder_count); - // The octree is built, the dataset is sorted but only the leaves are indexed: - // we propagate the indexing information bottom-up to the root - IndexInnerNodes( TemplatedOctree::Root() ); - } //end of Set + std::sort(placeholders.begin(), placeholders.end(), ObjectSorter< NodeType >()); + std::vector< NodePointer > filled_leaves(placeholder_count); + sorted_dataset.resize( placeholder_count ); + for (int i=0; i - ObjectPointer GetClosest - ( - OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, - OBJECT_MARKER & /*marker*/, - const CoordType & query_point, - const ScalarType & max_distance, - ScalarType & distance, - CoordType & point, - bool allow_zero_distance = true - ) - { - BoundingBoxType query_bb; - ScalarType sphere_radius; - if (!GuessInitialBoundingBox(query_point, max_distance, sphere_radius, query_bb)) - return NULL; - - std::vector< NodePointer > leaves; + // The dataset is sorted and the octree is built, but the indexing information aren't stored yet in the octree: + // we assign to each leaf the range inside the sorted dataset of the primitives contained inside the leaf + int begin = -1; + NodePointer initial_leaf = NULL; + for (int end=0; endSetRange(begin, end); + } - IncrementMark(); - AdjustBoundingBox(query_bb, sphere_radius, max_distance, leaves, 1); + // The octree is built, the dataset is sorted but only the leaves are indexed: + // we propagate the indexing information bottom-up to the root + IndexInnerNodes( TemplatedOctree::Root() ); + } //end of Set - if (sphere_radius>max_distance) - return NULL; + /*! + * Finds the closest object to a given point. + */ + template + ObjectPointer GetClosest + ( + OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, + OBJECT_MARKER & /*marker*/, + const CoordType & query_point, + const ScalarType & max_distance, + ScalarType & distance, + CoordType & point, + bool allow_zero_distance = true + ) + { + BoundingBoxType query_bb; + ScalarType sphere_radius; + if (!GuessInitialBoundingBox(query_point, max_distance, sphere_radius, query_bb)) + return NULL; - std::vector< Neighbour > neighbors; - RetrieveContainedObjects(query_point, distance_functor, max_distance, allow_zero_distance, leaves, neighbors); - - typename std::vector< Neighbour >::iterator first = neighbors.begin(); - typename std::vector< Neighbour >::iterator last = neighbors.end(); - std::partial_sort(first, first+1, last); - - distance = neighbors[0].distance; - point = neighbors[0].point; - return neighbors[0].object; - }; //end of GetClosest + std::vector< NodePointer > leaves; - /*! - * Retrieve the k closest element to the query point - */ - template - unsigned int GetKClosest - ( - OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, - OBJECT_MARKER & /*marker*/, - unsigned int k, - const CoordType & query_point, - const ScalarType & max_distance, - OBJECT_POINTER_CONTAINER & objects, - DISTANCE_CONTAINER & distances, - POINT_CONTAINER & points, - bool sort_per_distance = true, - bool allow_zero_distance = true - ) - { - BoundingBoxType query_bb; - ScalarType sphere_radius; - if (!GuessInitialBoundingBox(query_point, max_distance, sphere_radius, query_bb)) - return 0; + //unsigned int object_count; + //int leaves_count; - std::vector< NodePointer > leaves; - std::vector< Neighbour > neighbors; + IncrementMark(); + AdjustBoundingBox(query_bb, sphere_radius, max_distance, leaves, 1); - unsigned int object_count; - float k_distance; + if (sphere_radius>max_distance) + return NULL; + + std::vector< Neighbour > neighbors; + RetrieveContainedObjects(query_point, distance_functor, max_distance, allow_zero_distance, leaves, neighbors); + + typename std::vector< Neighbour >::iterator first = neighbors.begin(); + typename std::vector< Neighbour >::iterator last = neighbors.end(); + std::partial_sort(first, first+1, last); + + distance = neighbors[0].distance; + point = neighbors[0].point; + return neighbors[0].object; + }; //end of GetClosest + + /*! + * Retrieve the k closest element to the query point + */ + template + unsigned int GetKClosest + ( + OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, + OBJECT_MARKER & /*marker*/, + unsigned int k, + const CoordType & query_point, + const ScalarType & max_distance, + OBJECT_POINTER_CONTAINER & objects, + DISTANCE_CONTAINER & distances, + POINT_CONTAINER & points, + bool sort_per_distance = true, + bool allow_zero_distance = true + ) + { + BoundingBoxType query_bb; + ScalarType sphere_radius; + if (!GuessInitialBoundingBox(query_point, max_distance, sphere_radius, query_bb)) + return 0; + + std::vector< NodePointer > leaves; + std::vector< Neighbour > neighbors; + + unsigned int object_count; + float k_distance; OBJECT_RETRIEVER: - IncrementMark(); - AdjustBoundingBox(query_bb, sphere_radius, max_distance, leaves, k); - object_count = RetrieveContainedObjects(query_point, distance_functor, max_distance, allow_zero_distance, leaves, neighbors); + IncrementMark(); + AdjustBoundingBox(query_bb, sphere_radius, max_distance, leaves, k); + object_count = RetrieveContainedObjects(query_point, distance_functor, max_distance, allow_zero_distance, leaves, neighbors); - if (sphere_radius(first, first+object_count, last ); - else std::nth_element < NeighbourIterator >(first, first+object_count, last ); - - k_distance = neighbors[object_count-1].distance; - if (k_distance>sphere_radius && sphere_radius(first, first+object_count, last ); + else std::nth_element < NeighbourIterator >(first, first+object_count, last ); - return CopyQueryResults(neighbors, object_count, objects, distances, points); - }; //end of GetKClosest + k_distance = neighbors[object_count-1].distance; + if (k_distance>sphere_radius && sphere_radius(neighbors, object_count, objects, distances, points); + }; //end of GetKClosest - /*! - * Returns all the objects contained inside a specified sphere - */ - template - unsigned int GetInSphere - ( - OBJECT_POINT_DISTANCE_FUNCTOR &distance_functor, - OBJECT_MARKER &/*marker*/, - const CoordType &sphere_center, - const ScalarType &sphere_radius, - OBJECT_POINTER_CONTAINER &objects, - DISTANCE_CONTAINER &distances, - POINT_CONTAINER &points, - bool sort_per_distance = false, - bool allow_zero_distance = false - ) - { - // Define the minimum bounding-box containing the sphere - BoundingBoxType query_bb(sphere_center, sphere_radius); - - // If that bounding-box don't collide with the octree bounding-box, simply return 0 - if (!TemplatedOctree::boundingBox.Collide(query_bb)) - return 0; + /*! + * Returns all the objects contained inside a specified sphere + */ + template + unsigned int GetInSphere + ( + OBJECT_POINT_DISTANCE_FUNCTOR &distance_functor, + OBJECT_MARKER &/*marker*/, + const CoordType &sphere_center, + const ScalarType &sphere_radius, + OBJECT_POINTER_CONTAINER &objects, + DISTANCE_CONTAINER &distances, + POINT_CONTAINER &points, + bool sort_per_distance = false, + bool allow_zero_distance = false + ) + { + // Define the minimum bounding-box containing the sphere + BoundingBoxType query_bb(sphere_center, sphere_radius); - std::vector< NodePointer > leaves; - std::vector< Neighbour > neighbors; + // If that bounding-box don't collide with the octree bounding-box, simply return 0 + if (!TemplatedOctree::boundingBox.Collide(query_bb)) + return 0; - IncrementMark(); - ContainedLeaves(query_bb, leaves, TemplatedOctree::Root(), TemplatedOctree::boundingBox); + std::vector< NodePointer > leaves; + std::vector< Neighbour > neighbors; - int leaves_count = int(leaves.size()); - if (leaves_count==0) - return 0; + IncrementMark(); + ContainedLeaves(query_bb, leaves, TemplatedOctree::Root(), TemplatedOctree::boundingBox); - int object_count = RetrieveContainedObjects(sphere_center, distance_functor, sphere_radius, allow_zero_distance, leaves, neighbors); + int leaves_count = int(leaves.size()); + if (leaves_count==0) + return 0; - NeighbourIterator first = neighbors.begin(); - NeighbourIterator last = neighbors.end(); - if (sort_per_distance) std::partial_sort< NeighbourIterator >(first, first+object_count, last ); - else std::nth_element < NeighbourIterator >(first, first+object_count, last ); - - return CopyQueryResults(neighbors, object_count, objects, distances, points); - };//end of GetInSphere + int object_count = RetrieveContainedObjects(sphere_center, distance_functor, sphere_radius, allow_zero_distance, leaves, neighbors); - /*! - * Returns all the objects lying inside the specified bbox - */ - template - unsigned int GetInBox - ( - OBJECT_MARKER &/*marker*/, - const BoundingBoxType &query_bounding_box, - OBJECT_POINTER_CONTAINER &objects - ) - { - //if the query bounding-box don't collide with the octree bounding-box, simply return 0 - if (!query_bounding_box.Collide()) - { - objects.clear(); - return 0; - } - - //otherwise, retrieve the leaves and fill the container with the objects contained - std::vector< NodePointer > leaves; - unsigned int object_count; - int leaves_count; - - TemplatedOctree::ContainedLeaves(query_bounding_box, leaves, TemplatedOctree::Root(), TemplatedOctree::boundingBox); - leaves_count = int(leaves.size()); - if (leaves_count==0) - return 0; - - IncrementMark(); - for (int i=0; ibegin; - int end = voxel->end; - for ( ; beginpObject); - } //end of for ( ; begin(first, first+object_count, last ); + else std::nth_element < NeighbourIterator >(first, first+object_count, last ); - return int(objects.size()); - }; //end of GetInBox + return CopyQueryResults(neighbors, object_count, objects, distances, points); + };//end of GetInSphere - protected: - /*! - * Contains pointer to the objects in the dataset. - * The pointers are sorted so that the object pointed to result ordered in the space - */ - std::vector< ObjectReference > sorted_dataset; + /*! + * Returns all the objects lying inside the specified bbox + */ + template + unsigned int GetInBox + ( + OBJECT_MARKER &/*marker*/, + const BoundingBoxType &query_bounding_box, + OBJECT_POINTER_CONTAINER &objects + ) + { + //if the query bounding-box don't collide with the octree bounding-box, simply return 0 + if (!query_bounding_box.Collide()) + { + objects.clear(); + return 0; + } - /*! - * Markers used to avoid duplication of the same result during a query - */ - unsigned char *marks; - unsigned char global_mark; + //otherwise, retrieve the leaves and fill the container with the objects contained + std::vector< NodePointer > leaves; + unsigned int object_count; + int leaves_count; - /*! - * The expansion factor used to solve the spatial queries - * The current expansion factor is computed on the basis of the last expansion factor - * and on the history of these values, through the following heuristic: - * current_expansion_factor = alpha*last_expansion_factor + (1.0f-alpha)*mean_expansion_factor - * where alpha = 1.0/3.0; - */ - //float last_expansion_factor; - //float mean_expansion_factor; - //float ALPHA; - //float ONE_MINUS_ALPHA; + TemplatedOctree::ContainedLeaves(query_bounding_box, leaves, TemplatedOctree::Root(), TemplatedOctree::boundingBox); + leaves_count = int(leaves.size()); + if (leaves_count==0) + return 0; - protected: - /*! - */ - inline void IncrementMark() - { - // update the marks - global_mark = (global_mark+1)%255; - if (global_mark == 0) - { - memset(&marks[0], 0, sizeof(unsigned char)*int(sorted_dataset.size())); - global_mark++; - } - };//end of IncrementMark + IncrementMark(); + for (int i=0; ibegin; + int end = voxel->end; + for ( ; beginpMark == global_mark; }; + Mark(ref); + objects.push_back(ref->pObject); + } //end of for ( ; beginpMark = global_mark;}; + return int(objects.size()); + }; //end of GetInBox - /*! - * Guess an initial bounding-box from which starting the research of the closest point(s). - * \return true iff it's possibile to find a sphere, centered in query_point and having radius max_distance at most, which collide the octree bounding-box. - */ - inline bool GuessInitialBoundingBox(const CoordType &query_point, const ScalarType max_distance, ScalarType &sphere_radius, BoundingBoxType &query_bb) - { - // costruisco una bounging box centrata in query di dimensione pari a quella di una foglia. - // e controllo se in tale bounging box sono contenute un numero di elementi >= a k. - // Altrimenti espando il bounding box. - query_bb.Set(query_point); + protected: + /*! + * Contains pointer to the objects in the dataset. + * The pointers are sorted so that the object pointed to result ordered in the space + */ + std::vector< ObjectReference > sorted_dataset; - // the radius of the sphere centered in query - sphere_radius = 0.0f; + /*! + * Markers used to avoid duplication of the same result during a query + */ + unsigned char *marks; + unsigned char global_mark; - // if the bounding-box doesn't intersect the bounding-box of the octree, then it must be immediately expanded - if (!query_bb.IsIn(query_point)) - { - do - { - query_bb.Offset(TemplatedOctree::leafDiagonal); - sphere_radius += TemplatedOctree::leafDiagonal; - } - while ( !TemplatedOctree::boundingBox.Collide(query_bb) || sphere_radius>max_distance); - } - return (sphere_radius<=max_distance); - }; + /*! + * The expansion factor used to solve the spatial queries + * The current expansion factor is computed on the basis of the last expansion factor + * and on the history of these values, through the following heuristic: + * current_expansion_factor = alpha*last_expansion_factor + (1.0f-alpha)*mean_expansion_factor + * where alpha = 1.0/3.0; + */ + //float last_expansion_factor; + //float mean_expansion_factor; + //float ALPHA; + //float ONE_MINUS_ALPHA; - /*! - * Modify the bounding-box used during the query until either at least k points - * are contained inside the box or the box radius became greater than the threshold distance - * Return the number of leaves contained inside the bounding-box - */ - inline int AdjustBoundingBox - ( - BoundingBoxType & query_bb, - ScalarType & sphere_radius, - const ScalarType max_allowed_distance, - std::vector< NodePointer > & leaves, - const int required_object_count - ) - { - int leaves_count; - int object_count; - do - { - leaves.clear(); + protected: + /*! + */ + inline void IncrementMark() + { + // update the marks + global_mark = (global_mark+1)%255; + if (global_mark == 0) + { + memset(&marks[0], 0, sizeof(unsigned char)*int(sorted_dataset.size())); + global_mark++; + } + };//end of IncrementMark - query_bb.Offset(TemplatedOctree::leafDiagonal); - sphere_radius+= TemplatedOctree::leafDiagonal; + /* + */ + inline bool IsMarked(const ObjectReference *ref) const + { return *ref->pMark == global_mark; }; + + /* + */ + inline void Mark(const ObjectReference *ref) + { *ref->pMark = global_mark;}; + + /*! + * Guess an initial bounding-box from which starting the research of the closest point(s). + * \return true iff it's possibile to find a sphere, centered in query_point and having radius max_distance at most, which collide the octree bounding-box. + */ + inline bool GuessInitialBoundingBox(const CoordType &query_point, const ScalarType max_distance, ScalarType &sphere_radius, BoundingBoxType &query_bb) + { + // costruisco una bounging box centrata in query di dimensione pari a quella di una foglia. + // e controllo se in tale bounging box sono contenute un numero di elementi >= a k. + // Altrimenti espando il bounding box. + query_bb.Set(query_point); + + // the radius of the sphere centered in query + sphere_radius = 0.0f; + + // if the bounding-box doesn't intersect the bounding-box of the octree, then it must be immediately expanded + if (!query_bb.IsIn(query_point)) + { + do + { + query_bb.Offset(TemplatedOctree::leafDiagonal); + sphere_radius += TemplatedOctree::leafDiagonal; + } + while ( !TemplatedOctree::boundingBox.Collide(query_bb) || sphere_radius>max_distance); + } + return (sphere_radius<=max_distance); + }; + + /*! + * Modify the bounding-box used during the query until either at least k points + * are contained inside the box or the box radius became greater than the threshold distance + * Return the number of leaves contained inside the bounding-box + */ + inline int AdjustBoundingBox + ( + BoundingBoxType & query_bb, + ScalarType & sphere_radius, + const ScalarType max_allowed_distance, + std::vector< NodePointer > & leaves, + const int required_object_count + ) + { + int leaves_count; + int object_count; + do + { + leaves.clear(); + + query_bb.Offset(TemplatedOctree::leafDiagonal); + sphere_radius+= TemplatedOctree::leafDiagonal; TemplatedOctree::ContainedLeaves(query_bb, leaves, TemplatedOctree::Root(), TemplatedOctree::boundingBox); - leaves_count = int(leaves.size()); - object_count = 0; - for (int i=0; icount; - } - while (object_countcount; + } + while (object_count - inline int RetrieveContainedObjects - ( - const CoordType query_point, - OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, - const ScalarType max_allowed_distance, - bool allow_zero_distance, - std::vector< NodePointer > & leaves, - std::vector< Neighbour > & neighbors - ) - { - CoordType closest_point; - neighbors.clear(); - for (int i=0, leaves_count=int(leaves.size()); ibegin; - int end = voxel->end; - for ( ; begin + inline int RetrieveContainedObjects + ( + const CoordType query_point, + OBJECT_POINT_DISTANCE_FUNCTOR & distance_functor, + const ScalarType max_allowed_distance, + bool allow_zero_distance, + std::vector< NodePointer > & leaves, + std::vector< Neighbour > & neighbors + ) + { + CoordType closest_point; + neighbors.clear(); + for (int i=0, leaves_count=int(leaves.size()); ibegin; + int end = voxel->end; + for ( ; beginpObject, query_point, distance, closest_point)) - continue; + ScalarType distance = max_allowed_distance; + if (!distance_functor(*ref->pObject, query_point, distance, closest_point)) + continue; - Mark(ref); - if ((distance!=ScalarType(0.0) || allow_zero_distance)) - neighbors.push_back( Neighbour(ref->pObject, closest_point, distance) ); - } //end of for ( ; beginpObject, closest_point, distance) ); + } //end of for ( ; begin - inline int CopyQueryResults - ( - std::vector< Neighbour > &neighbors, - const unsigned int object_count, - OBJECT_POINTER_CONTAINER &objects, - DISTANCE_CONTAINER &distances, - POINT_CONTAINER &points - ) - { - // copy the nearest object into - points.resize( object_count ); - distances.resize( object_count ); - objects.resize( object_count ); + /*! + * Copy the results of a query + */ + template + inline int CopyQueryResults + ( + std::vector< Neighbour > &neighbors, + const unsigned int object_count, + OBJECT_POINTER_CONTAINER &objects, + DISTANCE_CONTAINER &distances, + POINT_CONTAINER &points + ) + { + // copy the nearest object into + points.resize( object_count ); + distances.resize( object_count ); + objects.resize( object_count ); - typename POINT_CONTAINER::iterator iPoint = points.begin(); - typename DISTANCE_CONTAINER::iterator iDistance = distances.begin(); - typename OBJECT_POINTER_CONTAINER::iterator iObject = objects.begin(); - for (unsigned int n=0; nAddRange( son_voxel ); - } - } - }; // end of IndexInnerNodes - }; + son_voxel = TemplatedOctree::Voxel(son_index); + current_voxel->AddRange( son_voxel ); + } + } + }; // end of IndexInnerNodes + }; #ifdef __glut_h__ - /************************************************************************/ - /* Rendering */ - /************************************************************************/ + /************************************************************************/ + /* Rendering */ + /************************************************************************/ protected: - /*! - * Structure which holds the rendering settings - */ - struct OcreeRenderingSetting - { - OcreeRenderingSetting() - { - color = vcg::Color4b(155, 155, 155, 255); - isVisible = false; - minVisibleDepth = 1; - maxVisibleDepth = 4; - }; + /*! + * Structure which holds the rendering settings + */ + struct OcreeRenderingSetting + { + OcreeRenderingSetting() + { + color = vcg::Color4b(155, 155, 155, 255); + isVisible = false; + minVisibleDepth = 1; + maxVisibleDepth = 4; + }; - int isVisible; - int minVisibleDepth; - int maxVisibleDepth; - vcg::Color4b color; - }; + int isVisible; + int minVisibleDepth; + int maxVisibleDepth; + vcg::Color4b color; + }; public: - /* - * Draw the octree in a valid OpenGL context according to the rendering settings - */ - void DrawOctree(vcg::Box3f &boundingBox, NodePointer n) - { - char level = TemplatedOctree::Level(n); - NodePointer son; - if (rendering_settings.minVisibleDepth>level) - { - for (int s=0; s<8; s++) - if ((son=Son(n, s))!=0) - DrawOctree(TemplatedOctree::SubBox(boundingBox, s), son); - } - else - { - vcg::glBoxWire(boundingBox); - if (levellevel) + { + for (int s=0; s<8; s++) + if ((son=Son(n, s))!=0) + DrawOctree(TemplatedOctree::SubBox(boundingBox, s), son); + } + else + { + vcg::glBoxWire(boundingBox); + if (level