manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Geodesic: Added possibility of saving also the implict tree of the shortest path. For each vertex you can give an attribute where the function will save the 'parent' e.g. the previous vertex in the shortest path to the closest source.

This commit is contained in:
Paolo Cignoni 2012-04-12 23:54:23 +00:00
parent 1c5f2c2264
commit ffdc2f2b28
1 changed files with 272 additions and 285 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

557
vcg/complex/algorithms/geodesic.h
View File

@ -47,342 +47,329 @@ g.FarthestVertex(m,seed,far,d);
#define __VCGLIB_GEODESIC #define __VCGLIB_GEODESIC
namespace vcg{ namespace vcg{
namespace tri{ namespace tri{
template <class MeshType> template <class MeshType>
struct EuclideanDistance{ struct EuclideanDistance{
typedef typename MeshType::VertexType VertexType; typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::ScalarType ScalarType;
EuclideanDistance(){} EuclideanDistance(){}
ScalarType operator()(const VertexType * v0, const VertexType * v1) const ScalarType operator()(const VertexType * v0, const VertexType * v1) const
{return vcg::Distance(v0->cP(),v1->cP());} {return vcg::Distance(v0->cP(),v1->cP());}
}; };
template <class MeshType, class DistanceFunctor = EuclideanDistance<MeshType> > template <class MeshType, class DistanceFunctor = EuclideanDistance<MeshType> >
class Geo{ class Geo{
public: public:
typedef typename MeshType::VertexType VertexType; typedef typename MeshType::VertexType VertexType;
typedef typename MeshType::VertexIterator VertexIterator; typedef typename MeshType::VertexIterator VertexIterator;
typedef typename MeshType::VertexPointer VertexPointer; typedef typename MeshType::VertexPointer VertexPointer;
typedef typename MeshType::FaceType FaceType; typedef typename MeshType::FaceType FaceType;
typedef typename MeshType::CoordType CoordType; typedef typename MeshType::CoordType CoordType;
typedef typename MeshType::ScalarType ScalarType; typedef typename MeshType::ScalarType ScalarType;
/* Auxiliary class for keeping the heap of vertices to visit and their estimated distance /* Auxiliary class for keeping the heap of vertices to visit and their estimated distance */
*/ struct VertDist{
struct VertDist{ VertDist(){}
VertDist(){} VertDist(VertexPointer _v, ScalarType _d):v(_v),d(_d){}
VertDist(VertexPointer _v, ScalarType _d):v(_v),d(_d){}
VertexPointer v; VertexPointer v;
ScalarType d; ScalarType d;
}; };
/* Temporary data to associate to all the vertices: estimated distance and boolean flag /* Temporary data to associate to all the vertices: estimated distance and boolean flag */
*/ struct TempData{
struct TempData{ TempData(){}
TempData(){} TempData(const ScalarType & _d):d(_d),source(0),parent(0){}
TempData(const ScalarType & d_){d=d_;source = NULL;}
ScalarType d;
VertexPointer source;//closest source
}; ScalarType d;
VertexPointer source;//closest source
VertexPointer parent;
};
typedef SimpleTempData<std::vector<VertexType>, TempData > TempDataType; typedef SimpleTempData<std::vector<VertexType>, TempData > TempDataType;
//TempDataType * TD;
struct pred: public std::binary_function<VertDist,VertDist,bool>{ struct pred: public std::binary_function<VertDist,VertDist,bool>{
pred(){}; pred(){}
bool operator()(const VertDist& v0, const VertDist& v1) const bool operator()(const VertDist& v0, const VertDist& v1) const
{return (v0.d > v1.d);} {return (v0.d > v1.d);}
}; };
struct pred_addr: public std::binary_function<VertDist,VertDist,bool>{
pred_addr(){};
bool operator()(const VertDist& v0, const VertDist& v1) const
{return (v0.v > v1.v);}
};
//************** calcolo della distanza di pw in base alle distanze note di pw1 e curr struct pred_addr: public std::binary_function<VertDist,VertDist,bool>{
//************** sapendo che (curr,pw,pw1) e'una faccia della mesh pred_addr(){}
//************** (vedi figura in file distance.gif) bool operator()(const VertDist& v0, const VertDist& v1) const
static ScalarType Distance(const VertexPointer &pw, {return (v0.v > v1.v);}
const VertexPointer &pw1, };
const VertexPointer &curr,
const ScalarType &d_pw1,
const ScalarType &d_curr)
{
ScalarType curr_d=0;
ScalarType ew_c = DistanceFunctor()(pw,curr); //************** calcolo della distanza di pw in base alle distanze note di pw1 e curr
ScalarType ew_w1 = DistanceFunctor()(pw,pw1); //************** sapendo che (curr,pw,pw1) e'una faccia della mesh
ScalarType ec_w1 = DistanceFunctor()(pw1,curr); //************** (vedi figura in file distance.gif)
CoordType w_c = (pw->cP()-curr->cP()).Normalize() * ew_c; static ScalarType Distance(const VertexPointer &pw,
CoordType w_w1 = (pw->cP() - pw1->cP()).Normalize() * ew_w1; const VertexPointer &pw1,
CoordType w1_c = (pw1->cP() - curr->cP()).Normalize() * ec_w1; const VertexPointer &curr,
const ScalarType &d_pw1,
const ScalarType &d_curr)
{
ScalarType curr_d=0;
ScalarType alpha,alpha_, beta,beta_,theta,h,delta,s,a,b; ScalarType ew_c = DistanceFunctor()(pw,curr);
ScalarType ew_w1 = DistanceFunctor()(pw,pw1);
ScalarType ec_w1 = DistanceFunctor()(pw1,curr);
CoordType w_c = (pw->cP()-curr->cP()).Normalize() * ew_c;
CoordType w_w1 = (pw->cP() - pw1->cP()).Normalize() * ew_w1;
CoordType w1_c = (pw1->cP() - curr->cP()).Normalize() * ec_w1;
alpha = acos((w_c.dot(w1_c))/(ew_c*ec_w1)); ScalarType alpha,alpha_, beta,beta_,theta,h,delta,s,a,b;
s = (d_curr + d_pw1+ec_w1)/2;
a = s/ec_w1;
b = a*s;
alpha_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_pw1)/d_curr)));
if ( alpha+alpha_ > M_PI){ alpha = acos((w_c.dot(w1_c))/(ew_c*ec_w1));
curr_d = d_curr + ew_c; s = (d_curr + d_pw1+ec_w1)/2;
}else a = s/ec_w1;
{ b = a*s;
beta_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_curr)/d_pw1))); alpha_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_pw1)/d_curr)));
beta = acos((w_w1).dot(-w1_c)/(ew_w1*ec_w1));
if ( beta+beta_ > M_PI) if ( alpha+alpha_ > M_PI){
curr_d = d_pw1 + ew_w1; curr_d = d_curr + ew_c;
else }else
{ {
theta = ScalarType(M_PI)-alpha-alpha_; beta_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_curr)/d_pw1)));
delta = cos(theta)* ew_c; beta = acos((w_w1).dot(-w1_c)/(ew_w1*ec_w1));
h = sin(theta)* ew_c;
curr_d = sqrt( pow(h,2)+ pow(d_curr + delta,2));
}
}
return (curr_d);
}
/* if ( beta+beta_ > M_PI)
This is the low level version of the geodesic computation framework. curr_d = d_pw1 + ew_w1;
Starting from the seeds, it assign a distance value to each vertex. The distance of a vertex is its else
approximated geodesic distance to the closest seeds. {
This is function is not meant to be called (although is not prevented). Instead, it is invoked by theta = ScalarType(M_PI)-alpha-alpha_;
wrapping function. delta = cos(theta)* ew_c;
*/ h = sin(theta)* ew_c;
static VertexPointer Visit( curr_d = sqrt( pow(h,2)+ pow(d_curr + delta,2));
MeshType & m, }
std::vector<VertDist> & seedVec, }
bool farthestOnBorder = false, return (curr_d);
ScalarType distance_threshold = std::numeric_limits<ScalarType>::max(), }
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * vertSource = NULL,
std::vector<VertexPointer> *InInterval=NULL)
{
bool isLeaf;
std::vector<VertDist> frontier;
VertexPointer curr,farthest=0,pw1;
ScalarType unreached = std::numeric_limits<ScalarType>::max();
VertexPointer pw; /*
This is the low level version of the geodesic computation framework.
Starting from the seeds, it assign a distance value to each vertex. The distance of a vertex is its
approximated geodesic distance to the closest seeds.
This is function is not meant to be called (although is not prevented). Instead, it is invoked by
wrapping function.
*/
static VertexPointer Visit(
MeshType & m,
std::vector<VertDist> & seedVec, // the set of seed to start from
bool farthestOnBorder = false,
ScalarType distance_threshold = std::numeric_limits<ScalarType>::max(), // cut off distance (do no compute anything farther than this value)
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * vertSource = NULL, // if present we put in this attribute the closest source for each vertex
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * vertParent = NULL, // if present we put in this attribute the parent in the path that goes from the vertex to the closest source
std::vector<VertexPointer> *InInterval=NULL)
{
std::vector<VertDist> frontier;
VertexPointer farthest=0,pw,pw1;
//Requirements //Requirements
assert(HasPerVertexVFAdjacency(m) && HasPerFaceVFAdjacency(m)); assert(HasPerVertexVFAdjacency(m) && HasPerFaceVFAdjacency(m));
assert(!seedVec.empty()); assert(!seedVec.empty());
TempDataType TD(m.vert,unreached); TempDataType TD(m.vert, std::numeric_limits<ScalarType>::max());
typename std::vector <VertDist >::iterator ifr; typename std::vector <VertDist >::iterator ifr;
for(ifr = seedVec.begin(); ifr != seedVec.end(); ++ifr){ for(ifr = seedVec.begin(); ifr != seedVec.end(); ++ifr){
TD[(*ifr).v].d = 0.0; (*ifr).d = 0.0;
(*ifr).d = 0.0; TD[(*ifr).v].d = 0.0;
TD[(*ifr).v].source = (*ifr).v; TD[(*ifr).v].source = (*ifr).v;
frontier.push_back(VertDist((*ifr).v,0.0)); TD[(*ifr).v].parent = (*ifr).v;
} frontier.push_back(VertDist((*ifr).v,0.0));
// initialize Heap }
make_heap(frontier.begin(),frontier.end(),pred()); // initialize Heap
make_heap(frontier.begin(),frontier.end(),pred());
ScalarType curr_d,d_curr = 0.0,d_heap; ScalarType curr_d,d_curr = 0.0,d_heap;
VertexPointer curr_s = NULL; ScalarType max_distance=0.0;
ScalarType max_distance=0.0;
typename std::vector<VertDist >:: iterator iv;
while(!frontier.empty() && max_distance < distance_threshold) while(!frontier.empty() && max_distance < distance_threshold)
{ {
pop_heap(frontier.begin(),frontier.end(),pred()); pop_heap(frontier.begin(),frontier.end(),pred());
curr = (frontier.back()).v; VertexPointer curr = (frontier.back()).v;
if (InInterval!=NULL) if (InInterval!=NULL) InInterval->push_back(curr);
InInterval->push_back(curr);
curr_s = TD[curr].source; if(vertSource!=NULL) (*vertSource)[curr] = TD[curr].source;
if(vertSource!=NULL) if(vertParent!=NULL) (*vertParent)[curr] = TD[curr].parent;
(*vertSource)[curr] = curr_s;
d_heap = (frontier.back()).d;
frontier.pop_back();
assert(TD[curr].d <= d_heap); d_heap = (frontier.back()).d;
assert(curr_s != NULL); frontier.pop_back();
if(TD[curr].d < d_heap )// a vertex whose distance has been improved after it was inserted in the queue
continue;
assert(TD[curr].d == d_heap);
d_curr = TD[curr].d; assert(TD[curr].d <= d_heap);
if(TD[curr].d < d_heap )// a vertex whose distance has been improved after it was inserted in the queue
continue;
assert(TD[curr].d == d_heap);
isLeaf = (!farthestOnBorder || curr->IsB()); d_curr = TD[curr].d;
face::VFIterator<FaceType> x;int k; bool isLeaf = (!farthestOnBorder || curr->IsB());
for( x.f = curr->VFp(), x.z = curr->VFi(); x.f!=0; ++x ) face::VFIterator<FaceType> x;int k;
for(k=0;k<2;++k)
{
if(k==0) {
pw = x.f->V1(x.z);
pw1=x.f->V2(x.z);
}
else {
pw = x.f->V2(x.z);
pw1=x.f->V1(x.z);
}
const ScalarType & d_pw1 = TD[pw1].d; for( x.f = curr->VFp(), x.z = curr->VFi(); x.f!=0; ++x )
{ for(k=0;k<2;++k)
const ScalarType inter = DistanceFunctor()(curr,pw1);//(curr->P() - pw1->P()).Norm(); {
const ScalarType tol = (inter + d_curr + d_pw1)*.0001f; if(k==0) {
pw = x.f->V1(x.z);
pw1=x.f->V2(x.z);
}
else {
pw = x.f->V2(x.z);
pw1=x.f->V1(x.z);
}
if ( (TD[pw1].source != TD[curr].source)||// not the same source const ScalarType & d_pw1 = TD[pw1].d;
(inter + d_curr < d_pw1 +tol ) || {
(inter + d_pw1 < d_curr +tol ) || const ScalarType inter = DistanceFunctor()(curr,pw1);//(curr->P() - pw1->P()).Norm();
(d_curr + d_pw1 < inter +tol ) // triangular inequality const ScalarType tol = (inter + d_curr + d_pw1)*.0001f;
)
curr_d = d_curr + DistanceFunctor()(pw,curr);//(pw->P()-curr->P()).Norm();
else
curr_d = Distance(pw,pw1,curr,d_pw1,d_curr);
}
if(TD[(pw)].d > curr_d){ if ( (TD[pw1].source != TD[curr].source)||// not the same source
TD[(pw)].d = curr_d; (inter + d_curr < d_pw1 +tol ) ||
TD[pw].source = curr_s; (inter + d_pw1 < d_curr +tol ) ||
frontier.push_back(VertDist(pw,curr_d)); (d_curr + d_pw1 < inter +tol ) // triangular inequality
push_heap(frontier.begin(),frontier.end(),pred()); )
} curr_d = d_curr + DistanceFunctor()(pw,curr);//(pw->P()-curr->P()).Norm();
if(isLeaf){ else
if(d_curr > max_distance){ curr_d = Distance(pw,pw1,curr,d_pw1,d_curr);
max_distance = d_curr; }
farthest = curr;
}
}
}
}// end while
// Copy found distance onto the Quality (\todo parametric!) if(TD[pw].d > curr_d){
if (InInterval==NULL) TD[pw].d = curr_d;
{ TD[pw].source = TD[curr].source;
for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) if(!(*vi).IsD()) TD[pw].parent = curr;
(*vi).Q() = TD[&(*vi)].d; frontier.push_back(VertDist(pw,curr_d));
} push_heap(frontier.begin(),frontier.end(),pred());
else }
{ if(isLeaf){
assert(InInterval->size()>0); if(d_curr > max_distance){
for(size_t i=0;i<InInterval->size();i++) max_distance = d_curr;
(*InInterval)[i]->Q() = TD[(*InInterval)[i]].d; farthest = curr;
} }
}
}
}// end while
return farthest; // Copy found distance onto the Quality (\todo parametric!)
} if (InInterval==NULL)
{
for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi) if(!(*vi).IsD())
(*vi).Q() = TD[&(*vi)].d;
}
else
{
assert(InInterval->size()>0);
for(size_t i=0;i<InInterval->size();i++)
(*InInterval)[i]->Q() = TD[(*InInterval)[i]].d;
}
return farthest;
}
public: public:
/* /*
Given a mesh and a vector of pointers to seed vertices, this function assigns the approximated geodesic Given a mesh and a vector of pointers to seed vertices, this function assigns the approximated geodesic
distance from the closest source to all the mesh vertices within the distance from the closest source to all the mesh vertices within the
specified interval and returns the found vertices writing on their Quality field the distance. specified interval and returns the found vertices writing on their Quality field the distance.
Optionally for each vertex it can store, in a passed attribute, its corresponding seed vertex. Optionally for each vertex it can store, in a passed attribute, its corresponding seed vertex.
To allocate such an attribute: To allocate such an attribute:
typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources; typename MeshType::template PerVertexAttributeHandle<VertexPointer> sources;
sources = tri::Allocator<CMeshO>::AddPerVertexAttribute<VertexPointer> (m,"sources"); sources = tri::Allocator<CMeshO>::AddPerVertexAttribute<VertexPointer> (m,"sources");
*/ */
static bool FarthestVertex( MeshType & m, static bool FarthestVertex( MeshType & m,
std::vector<VertexPointer> & seedVec, std::vector<VertexPointer> & seedVec,
VertexPointer & farthest_vert, VertexPointer & farthest_vert,
ScalarType distance_thr = std::numeric_limits<ScalarType>::max(), ScalarType distance_thr = std::numeric_limits<ScalarType>::max(),
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sourceSeed = NULL, typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sourceSeed = NULL,
std::vector<VertexPointer> *InInterval=NULL) typename MeshType::template PerVertexAttributeHandle<VertexPointer> * parentSeed = NULL,
{ std::vector<VertexPointer> *InInterval=NULL)
typename std::vector<VertexPointer>::iterator fi; {
std::vector<VertDist> vdSeedVec; typename std::vector<VertexPointer>::iterator fi;
if(seedVec.empty()) return false; std::vector<VertDist> vdSeedVec;
for( fi = seedVec.begin(); fi != seedVec.end() ; ++fi) if(seedVec.empty()) return false;
{ for( fi = seedVec.begin(); fi != seedVec.end() ; ++fi)
vdSeedVec.push_back(VertDist(*fi,0.0)); vdSeedVec.push_back(VertDist(*fi,0.0));
/* if (InInterval==NULL)continue; farthest_vert = Visit(m, vdSeedVec, false, distance_thr, sourceSeed, parentSeed, InInterval);
InInterval.push_back();*/ return true;
} }
farthest_vert = Visit(m, vdSeedVec, false, distance_thr, sourceSeed, InInterval); /*
return true; Given a mesh and a pointers to a vertex-source (source), assigns the approximated geodesic
} distance from the vertex-source to all the mesh vertices and returns the pointer to the farthest
/* Note: it updates the field Q() of the vertices
Given a mesh and a pointers to a vertex-source (source), assigns the approximated geodesic */
distance from the vertex-source to all the mesh vertices and returns the pointer to the farthest static bool FarthestVertex( MeshType & m, VertexPointer seed, ScalarType distance_thr = std::numeric_limits<ScalarType>::max())
Note: update the field Q() of the vertices {
*/ std::vector<VertexPointer> seedVec(1,seed);
static bool FarthestVertex( MeshType & m, VertexPointer v0;
VertexPointer seed, return FarthestVertex(m,seedVec,v0,distance_thr);
ScalarType distance_thr = std::numeric_limits<ScalarType>::max()) }
{
std::vector<VertexPointer> seedVec;
seedVec.push_back( seed );
VertexPointer v0;
return FarthestVertex(m,seedVec,v0,distance_thr);
}
/*
Same as FarthestPoint but the returned pointer is to a border vertex
Note: update the field Q() of the vertices
*/
static void FarthestBVertex(MeshType & m,
std::vector<VertexPointer> & seedVec,
VertexPointer & farthest,
ScalarType & distance,
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
){
typename std::vector<VertexPointer>::iterator fi; /*
std::vector<VertDist>fr; Same as FarthestPoint but the returned pointer is to a border vertex
Note: update the field Q() of the vertices
*/
static void FarthestBVertex(MeshType & m,
std::vector<VertexPointer> & seedVec,
VertexPointer & farthest,
ScalarType & distance,
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
)
{
std::vector<VertDist>fr;
for(typename std::vector<VertexPointer>::iterator fi = seedVec.begin(); fi != seedVec.end() ; ++fi)
fr.push_back(VertDist(*fi,0));
farthest = Visit(m,fr,distance,true,sources);
}
/*
Same as FarthestPoint but the returned pointer is to a border vertex
Note: update the field Q() of the vertices
*/
static void FarthestBVertex( MeshType & m,
VertexPointer seed,
VertexPointer & farthest,
ScalarType & distance,
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL)
{
std::vector<VertexPointer> fro(1,seed);
VertexPointer v0;
FarthestBVertex(m,fro,v0,distance,sources);
farthest = v0;
}
for( fi = seedVec.begin(); fi != seedVec.end() ; ++fi) /*
fr.push_back(VertDist(*fi,-1)); Assigns to each vertex of the mesh its distance to the closest vertex on the border
farthest = Visit(m,fr,distance,true,sources); Note: update the field Q() of the vertices
} Note: it needs the border bit set.
/* */
Same as FarthestPoint but the returned pointer is to a border vertex static bool DistanceFromBorder( MeshType & m, typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
Note: update the field Q() of the vertices ){
*/ std::vector<VertexPointer> fro;
static void FarthestBVertex( MeshType & m, VertexIterator vi;
VertexPointer seed, VertexPointer farthest;
VertexPointer & farthest, for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
ScalarType & distance, if( (*vi).IsB())
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL){ fro.push_back(&(*vi));
std::vector<VertexPointer> fro; if(fro.empty()) return false;
fro.push_back( seed );
VertexPointer v0;
FarthestBVertex(m,fro,v0,distance,sources);
farthest = v0;
}
/* tri::UpdateQuality<MeshType>::VertexConstant(m,0);
Assigns to each vertex of the mesh its distance to the closest vertex on the border
Note: update the field Q() of the vertices
*/
static bool DistanceFromBorder( MeshType & m,
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
){
std::vector<VertexPointer> fro;
VertexIterator vi;
VertexPointer farthest;
for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if( (*vi).IsB())
fro.push_back(&(*vi));
if(fro.empty()) return false;
tri::UpdateQuality<MeshType>::VertexConstant(m,0); return FarthestVertex(m,fro,farthest,std::numeric_limits<ScalarType>::max(),sources);
}
return FarthestVertex(m,fro,farthest,std::numeric_limits<ScalarType>::max(),sources); };// end class
} }// end namespace tri
}// end namespace vcg
};
};// end namespace tri
};// end namespace vcg
#endif #endif