Removed useless include
This commit is contained in:
parent
998312b65e
commit
6f7e2872af
|
|
@ -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 <vcg/complex/complex.h>
|
||||
#include <vcg/space/point_matching.h>
|
||||
#include <vcg/complex/algorithms/closest.h>
|
||||
#include <vcg/complex/complex.h>
|
||||
#include <wrap/io_trimesh/export_ply.h>
|
||||
|
||||
// note: temporary (callback.h should be moved inside vcg)
|
||||
|
|
@ -150,191 +150,191 @@ private:
|
|||
void GetBBox(vcg::Box3<ScalarType> & b){b.Add(pos);}
|
||||
};
|
||||
|
||||
GridType *ugrid; // griglia
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType > ugridQ;
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType > ugridP;
|
||||
GridType *ugrid; // griglia
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType > ugridQ;
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType > 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<ScalarType> & 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<ScalarType> & 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<vcg::Matrix44f> 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<vcg::Matrix44f> 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 <class MeshType>
|
||||
void FourPCS<MeshType>:: 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<ScalarType > dists;
|
||||
std::vector<VertexType* > ress;
|
||||
vcg::tri::GetKClosestVertex<
|
||||
MeshType,
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType>,
|
||||
std::vector<VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
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<ScalarType > dists;
|
||||
std::vector<VertexType* > ress;
|
||||
vcg::tri::GetKClosestVertex<
|
||||
MeshType,
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType>,
|
||||
std::vector<VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
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 <class MeshType>
|
||||
bool
|
||||
FourPCS<MeshType>::SelectCoplanarBase(){
|
||||
|
||||
vcg::tri::UpdateBounding<MeshType>::Box(*P);
|
||||
vcg::tri::UpdateBounding<MeshType>::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<float>::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<float>::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<typename MeshType::VertexType*> closests;
|
||||
std::vector<ScalarType> distances;
|
||||
std::vector<CoordType> points;
|
||||
std::vector<typename MeshType::VertexType*> closests;
|
||||
std::vector<ScalarType> distances;
|
||||
std::vector<CoordType> points;
|
||||
|
||||
vcg::tri::GetInSphereVertex<
|
||||
MeshType,
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType >,
|
||||
std::vector<typename MeshType::VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
std::vector<CoordType>
|
||||
>(*P,ugridP,B4,radius,closests,distances,points);
|
||||
vcg::tri::GetInSphereVertex<
|
||||
MeshType,
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType >,
|
||||
std::vector<typename MeshType::VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
std::vector<CoordType>
|
||||
>(*P,ugridP,B4,radius,closests,distances,points);
|
||||
|
||||
if(closests.empty())
|
||||
return false;
|
||||
best = -1; bestv=std::numeric_limits<float>::max();
|
||||
for(i = 0; i <closests.size(); ++i){
|
||||
ScalarType angle = fabs((closests[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<float>::max();
|
||||
for(i = 0; i <closests.size(); ++i){
|
||||
ScalarType angle = fabs((closests[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<ScalarType > dists;
|
||||
radius =side*0.5;
|
||||
std::vector< CoordType > samples,d_samples;
|
||||
std::vector<ScalarType > dists;
|
||||
|
||||
for(int i = 0 ; i< 4; ++i){
|
||||
vcg::tri::GetKClosestVertex<
|
||||
MeshType,
|
||||
vcg::GridStaticPtr<typename MeshType::VertexType, ScalarType >,
|
||||
std::vector<VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
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<typename MeshType::VertexType, ScalarType >,
|
||||
std::vector<VertexType*>,
|
||||
std::vector<ScalarType>,
|
||||
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<MeshType>::IsTransfCongruent(FourPoints fp, vcg::Matrix44<ScalarTyp
|
|||
template <class MeshType>
|
||||
void
|
||||
FourPCS<MeshType>::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 <class MeshType>
|
||||
bool FourPCS<MeshType>::FindCongruent() { // of base B, on Q, with approximation delta
|
||||
bool done = false;
|
||||
std::vector<EPoint> 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<EPoint> 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<Couple>::iterator bR1,eR1,bR2,eR2,ite,cite;
|
||||
bR1 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d1-par.delta*2.0));
|
||||
eR1 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d1+par.delta*2.0));
|
||||
bR2 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d2-par.delta*2.0));
|
||||
eR2 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d2+par.delta*2.0));
|
||||
typename PMesh::VertexIterator vii;
|
||||
typename std::vector<Couple>::iterator bR1,eR1,bR2,eR2,ite,cite;
|
||||
bR1 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d1-par.delta*2.0));
|
||||
eR1 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d1+par.delta*2.0));
|
||||
bR2 = std::lower_bound<typename std::vector<Couple>::iterator,Couple>(R1.begin(),R1.end(),Couple(d2-par.delta*2.0));
|
||||
eR2 = std::lower_bound<typename std::vector<Couple>::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<PMesh>::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<PMesh>::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<int> id = vcg::tri::Allocator<PMesh>::template AddPerVertexAttribute<int>(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<PMesh>::Box(Invr);
|
||||
// printf("Invr size %d\n",Invr.vn);
|
||||
vcg::tri::UpdateBounding<PMesh>::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<typename PMesh::VertexType*> closests;
|
||||
std::vector<typename PMesh::VertexType*> closests;
|
||||
|
||||
// for each point in R2inv get all the points in R1 closer than par.delta
|
||||
vcg::Matrix44<ScalarType> 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<ScalarType> 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<PMesh,GridType,std::vector<typename PMesh::VertexType*> >
|
||||
(Invr,*ugrid,bb,closests);
|
||||
vcg::tri::GetInBoxVertex<PMesh,GridType,std::vector<typename PMesh::VertexType*> >
|
||||
(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<PMesh>::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<PMesh>::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<MeshType>::EvaluateSample(Candidate & fp, CoordType & tp, CoordType
|
|||
template <class MeshType>
|
||||
void
|
||||
FourPCS<MeshType>::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 <class MeshType>
|
||||
void
|
||||
FourPCS<MeshType>::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 <class MeshType>
|
|||
bool
|
||||
FourPCS<MeshType>:: 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<float>::max();
|
||||
iwinner = 0;
|
||||
bestv = -std::numeric_limits<float>::max();
|
||||
iwinner = 0;
|
||||
|
||||
for(int i = 0 ; i < U.size() ;++i)
|
||||
{
|
||||
|
|
@ -635,15 +635,15 @@ FourPCS<MeshType>:: 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
|
||||
|
|
|
|||
|
|
@ -3,7 +3,6 @@
|
|||
|
||||
#include <iostream>
|
||||
#include <list>
|
||||
#include <wrap/callback.h>
|
||||
#include <vcg/complex/algorithms/update/topology.h>
|
||||
#include <vcg/complex/algorithms/update/flag.h>
|
||||
|
||||
|
|
|
|||
|
|
@ -22,7 +22,6 @@
|
|||
****************************************************************************/
|
||||
#ifndef __VCG_TRIVIAL_WALKER
|
||||
#define __VCG_TRIVIAL_WALKER
|
||||
#include <wrap/callback.h>
|
||||
|
||||
namespace vcg {
|
||||
|
||||
|
|
@ -56,36 +55,36 @@ public:
|
|||
//else return numeric_limits<float>::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<class EXTRACTOR_TYPE>
|
||||
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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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<MeshType>::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
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -33,7 +33,6 @@
|
|||
#include <vcg/complex/algorithms/update/topology.h>
|
||||
#include <vcg/complex/algorithms/update/flag.h>
|
||||
#include <vcg/space/triangle3.h>
|
||||
#include <wrap/callback.h>
|
||||
|
||||
namespace vcg{
|
||||
namespace tri{
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -33,7 +33,6 @@
|
|||
#include <vcg/complex/algorithms/intersection.h>
|
||||
#include <vcg/complex/algorithms/inertia.h>
|
||||
#include <eigenlib/Eigen/Core>
|
||||
#include <wrap/callback.h>
|
||||
|
||||
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<FaceType> 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<FaceType> 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<Eigen::Matrix3d> 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<Eigen::Matrix3d> 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 <a href="http://www2.in.tu-clausthal.de/~hormann/papers/Dyn.2001.OTU.pdf"> <em> "Optimizing 3d triangulations using discrete curvature analysis" </em> </a>
|
||||
*/
|
||||
/**
|
||||
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 <a href="http://www2.in.tu-clausthal.de/~hormann/papers/Dyn.2001.OTU.pdf"> <em> "Optimizing 3d triangulations using discrete curvature analysis" </em> </a>
|
||||
*/
|
||||
|
||||
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<float>::epsilon()) {
|
||||
v->Kh() = 0;
|
||||
v->Kg() = 0;
|
||||
}
|
||||
else {
|
||||
v->Kh() /= A;
|
||||
v->Kg() /= A;
|
||||
}
|
||||
}
|
||||
if(norm) {
|
||||
if(A <= std::numeric_limits<float>::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<ScalarType> m33;m33.SetZero();
|
||||
face::JumpingPos<typename MeshType::FaceType> 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<ScalarType> m33;m33.SetZero();
|
||||
face::JumpingPos<typename MeshType::FaceType> 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<ScalarType> normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P();
|
||||
ScalarType edge_length = normalized_edge.Norm();
|
||||
normalized_edge/=edge_length;
|
||||
Point3<ScalarType> n1 = p.F()->cN();n1.Normalize();
|
||||
Point3<ScalarType> 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<ScalarType> normalized_edge = p.F()->V(p.F()->Next(p.VInd()))->cP() - (*vi).P();
|
||||
ScalarType edge_length = normalized_edge.Norm();
|
||||
normalized_edge/=edge_length;
|
||||
Point3<ScalarType> n1 = p.F()->cN();n1.Normalize();
|
||||
Point3<ScalarType> 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<Eigen::Matrix3d> eig(it);
|
||||
Eigen::Vector3d c_val = eig.eigenvalues();
|
||||
Eigen::Matrix3d c_vec = eig.eigenvectors();
|
||||
Eigen::Matrix3d it;
|
||||
m33.ToEigenMatrix(it);
|
||||
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eig(it);
|
||||
Eigen::Vector3d c_val = eig.eigenvalues();
|
||||
Eigen::Matrix3d c_vec = eig.eigenvectors();
|
||||
|
||||
Point3<ScalarType> lambda;
|
||||
Matrix33<ScalarType> vect;
|
||||
vect.FromEigenMatrix(c_vec);
|
||||
lambda.FromEigenVector(c_val);
|
||||
Point3<ScalarType> lambda;
|
||||
Matrix33<ScalarType> 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<ScalarType>*)(& vect[maxI][0]);
|
||||
(*vi).PD2() = *(Point3<ScalarType>*)(& vect[minI][0]);
|
||||
(*vi).K1() = lambda[2];
|
||||
(*vi).K2() = lambda[1];
|
||||
}
|
||||
}
|
||||
(*vi).PD1() = *(Point3<ScalarType>*)(& vect[maxI][0]);
|
||||
(*vi).PD2() = *(Point3<ScalarType>*)(& 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<m.vert.size();++i) {
|
||||
CoordType dd = m.vert[i].P()-c;
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue