23 #include <vcg/simplex/face/pos.h>
24 #include <vcg/simplex/face/topology.h>
25 #include <vcg/complex/algorithms/update/quality.h>
28 #ifndef __VCGLIB_GEODESIC
29 #define __VCGLIB_GEODESIC
34 template <
class MeshType>
35 struct EuclideanDistance{
36 typedef typename MeshType::VertexType VertexType;
37 typedef typename MeshType::ScalarType ScalarType;
38 typedef typename MeshType::FacePointer FacePointer;
42 ScalarType operator()(
const VertexType * v0,
const VertexType * v1)
const
43 {
return vcg::Distance(v0->cP(),v1->cP());}
45 ScalarType operator()(
const FacePointer f0,
const FacePointer f1)
const
46 {
return vcg::Distance(Barycenter(*f0),Barycenter(*f1));}
50 template <
class MeshType>
51 class IsotropicDistance{
54 typename MeshType::template PerVertexAttributeHandle<float> wH;
56 typedef typename MeshType::VertexType VertexType;
57 typedef typename MeshType::ScalarType ScalarType;
58 typedef typename MeshType::FacePointer FacePointer;
66 IsotropicDistance(MeshType &m,
float variance)
69 wH = tri::Allocator<MeshType>:: template GetPerVertexAttribute<float> (m,
"distW");
70 float qmin = 1.0f/variance;
71 float qmax = variance;
72 float qrange = qmax-qmin;
73 std::pair<float,float> minmax = Stat<MeshType>::ComputePerVertexQualityMinMax(m);
74 float range = minmax.second-minmax.first;
75 for(
size_t i=0;i<m.vert.size();++i)
76 wH[i]=qmin+((m.vert[i].Q()-minmax.first)/range)*qrange;
81 ScalarType operator()( VertexType * v0, VertexType * v1)
83 float scale = (wH[v0]+wH[v1])/2.0f;
84 return (1.0f/scale)*vcg::Distance(v0->cP(),v1->cP());
89 template <
class MeshType>
90 struct BasicCrossFunctor
92 BasicCrossFunctor(MeshType &m) { tri::RequirePerVertexCurvatureDir(m); }
93 typedef typename MeshType::VertexType VertexType;
95 typename MeshType::CoordType D1(VertexType &v) {
return v.PD1(); }
96 typename MeshType::CoordType D2(VertexType &v) {
return v.PD2(); }
107 template <
class MeshType>
109 typedef typename MeshType::VertexType VertexType;
110 typedef typename MeshType::ScalarType ScalarType;
111 typedef typename MeshType::CoordType CoordType;
112 typedef typename MeshType::VertexIterator VertexIterator;
114 typename MeshType::template PerVertexAttributeHandle<CoordType> wxH,wyH;
121 for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
128 ScalarType operator()( VertexType * v0, VertexType * v1)
130 CoordType dd = v0->cP()-v1->cP();
131 float x = (fabs(dd * wxH[v0])+fabs(dd *wxH[v1]))/2.0f;
132 float y = (fabs(dd * wyH[v0])+fabs(dd *wyH[v1]))/2.0f;
134 return sqrt(x*x+y*y);
153 template <
class MeshType>
158 typedef typename MeshType::VertexType VertexType;
159 typedef typename MeshType::VertexIterator VertexIterator;
160 typedef typename MeshType::VertexPointer VertexPointer;
161 typedef typename MeshType::FacePointer FacePointer;
162 typedef typename MeshType::FaceType FaceType;
163 typedef typename MeshType::CoordType CoordType;
164 typedef typename MeshType::ScalarType ScalarType;
171 VertDist(VertexPointer _v, ScalarType _d):v(_v),d(_d){}
179 DIJKDist(VertexPointer _v):v(_v), q(_v->Q()){}
183 bool operator < (
const DIJKDist &o)
const
193 FaceDist(FacePointer _f):f(_f){}
195 bool operator < (
const FaceDist &o)
const
197 if( f->Q() != o.f->Q())
198 return f->Q() > o.f->Q();
206 TempData(
const ScalarType & _d):d(_d),source(0),parent(0){}
209 VertexPointer source;
210 VertexPointer parent;
213 typedef SimpleTempData<std::vector<VertexType>, TempData > TempDataType;
216 struct pred:
public std::binary_function<VertDist,VertDist,bool>{
218 bool operator()(
const VertDist& v0,
const VertDist& v1)
const
219 {
return (v0.d > v1.d);}
243 template <
class DistanceFunctor>
244 static ScalarType Distance(DistanceFunctor &distFunc,
245 const VertexPointer &pw,
246 const VertexPointer &pw1,
247 const VertexPointer &curr,
248 const ScalarType &d_pw1,
249 const ScalarType &d_curr)
253 ScalarType ew_c = distFunc(pw,curr);
254 ScalarType ew_w1 = distFunc(pw,pw1);
255 ScalarType ec_w1 = distFunc(pw1,curr);
256 CoordType w_c = (pw->cP()-curr->cP()).Normalize() * ew_c;
257 CoordType w_w1 = (pw->cP() - pw1->cP()).Normalize() * ew_w1;
258 CoordType w1_c = (pw1->cP() - curr->cP()).Normalize() * ec_w1;
260 ScalarType alpha,alpha_, beta,beta_,theta,h,delta,s,a,b;
262 alpha = acos((w_c.dot(w1_c))/(ew_c*ec_w1));
263 s = (d_curr + d_pw1+ec_w1)/2;
266 alpha_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_pw1)/d_curr)));
268 if ( alpha+alpha_ > M_PI){
269 curr_d = d_curr + ew_c;
272 beta_ = 2*acos ( std::min<ScalarType>(1.0,sqrt( (b- a* d_curr)/d_pw1)));
273 beta = acos((w_w1).dot(-w1_c)/(ew_w1*ec_w1));
275 if ( beta+beta_ > M_PI)
276 curr_d = d_pw1 + ew_w1;
279 theta = ScalarType(M_PI)-alpha-alpha_;
280 delta = cos(theta)* ew_c;
281 h = sin(theta)* ew_c;
282 curr_d = sqrt( pow(h,2)+ pow(d_curr + delta,2));
299 template <
class DistanceFunctor>
300 static VertexPointer Visit(
302 std::vector<VertDist> & seedVec,
303 DistanceFunctor &distFunc,
304 ScalarType distance_threshold = std::numeric_limits<ScalarType>::max(),
305 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * vertSource = NULL,
306 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * vertParent = NULL,
307 std::vector<VertexPointer> *InInterval=NULL)
309 VertexPointer farthest=0;
312 tri::RequireVFAdjacency(m);
313 tri::RequirePerVertexQuality(m);
315 assert(!seedVec.empty());
317 TempDataType TD(m.vert, std::numeric_limits<ScalarType>::max());
320 std::vector<VertDist> frontierHeap;
321 typename std::vector <VertDist >::iterator ifr;
322 for(ifr = seedVec.begin(); ifr != seedVec.end(); ++ifr){
323 TD[(*ifr).v].d = (*ifr).d;
324 TD[(*ifr).v].source = (*ifr).v;
325 TD[(*ifr).v].parent = (*ifr).v;
326 frontierHeap.push_back(*ifr);
328 make_heap(frontierHeap.begin(),frontierHeap.end(),pred());
330 ScalarType curr_d,d_curr = 0.0,d_heap;
331 ScalarType max_distance=0.0;
333 while(!frontierHeap.empty() && max_distance < distance_threshold)
335 pop_heap(frontierHeap.begin(),frontierHeap.end(),pred());
336 VertexPointer curr = (frontierHeap.back()).v;
337 if (InInterval!=NULL) InInterval->push_back(curr);
339 if(vertSource!=NULL) (*vertSource)[curr] = TD[curr].source;
340 if(vertParent!=NULL) (*vertParent)[curr] = TD[curr].parent;
342 d_heap = (frontierHeap.back()).d;
343 frontierHeap.pop_back();
345 assert(TD[curr].d <= d_heap);
346 if(TD[curr].d < d_heap )
348 assert(TD[curr].d == d_heap);
352 for(face::VFIterator<FaceType> vfi(curr) ; vfi.f!=0; ++vfi )
356 VertexPointer pw,pw1;
358 pw = vfi.f->V1(vfi.z);
359 pw1= vfi.f->V2(vfi.z);
362 pw = vfi.f->V2(vfi.z);
363 pw1= vfi.f->V1(vfi.z);
366 const ScalarType & d_pw1 = TD[pw1].d;
368 const ScalarType inter = distFunc(curr,pw1);
369 const ScalarType tol = (inter + d_curr + d_pw1)*.0001f;
371 if ( (TD[pw1].source != TD[curr].source)||
372 (inter + d_curr < d_pw1 +tol ) ||
373 (inter + d_pw1 < d_curr +tol ) ||
374 (d_curr + d_pw1 < inter +tol )
376 curr_d = d_curr + distFunc(pw,curr);
378 curr_d = Distance(distFunc,pw,pw1,curr,d_pw1,d_curr);
381 if(TD[pw].d > curr_d){
383 TD[pw].source = TD[curr].source;
384 TD[pw].parent = curr;
385 frontierHeap.push_back(VertDist(pw,curr_d));
386 push_heap(frontierHeap.begin(),frontierHeap.end(),pred());
389 if(d_curr > max_distance){
390 max_distance = d_curr;
400 if (InInterval==NULL)
402 for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
if(!(*vi).IsD())
403 (*vi).Q() = TD[&(*vi)].d;
407 assert(InInterval->size()>0);
408 for(
size_t i=0;i<InInterval->size();i++)
409 (*InInterval)[i]->Q() = TD[(*InInterval)[i]].d;
449 const std::vector<VertexPointer> & seedVec)
451 EuclideanDistance<MeshType> dd;
455 template <
class DistanceFunctor>
456 static bool Compute( MeshType & m,
457 const std::vector<VertexPointer> & seedVec,
458 DistanceFunctor &distFunc,
459 ScalarType maxDistanceThr = std::numeric_limits<ScalarType>::max(),
460 std::vector<VertexPointer> *withinDistanceVec=NULL,
461 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sourceSeed = NULL,
462 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * parentSeed = NULL
465 if(seedVec.empty())
return false;
466 std::vector<VertDist> vdSeedVec;
467 typename std::vector<VertexPointer>::const_iterator fi;
468 for( fi = seedVec.begin(); fi != seedVec.end() ; ++fi)
469 vdSeedVec.push_back(VertDist(*fi,0.0));
470 Visit(m, vdSeedVec, distFunc, maxDistanceThr, sourceSeed, parentSeed, withinDistanceVec);
481 static bool DistanceFromBorder( MeshType & m,
typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL)
483 std::vector<VertexPointer> fro;
484 for(VertexIterator vi = m.vert.begin(); vi != m.vert.end(); ++vi)
486 fro.push_back(&(*vi));
487 if(fro.empty())
return false;
488 EuclideanDistance<MeshType> dd;
490 return Compute(m,fro,dd,std::numeric_limits<ScalarType>::max(),0,sources);
494 static bool ConvertPerVertexSeedToPerFaceSeed(MeshType &m,
const std::vector<VertexPointer> &vertexSeedVec,
495 std::vector<FacePointer> &faceSeedVec)
497 tri::RequireVFAdjacency(m);
498 tri::RequirePerFaceMark(m);
502 for(
size_t i=0;i<vertexSeedVec.size();++i)
504 for(face::VFIterator<FaceType> vfi(vertexSeedVec[i]);!vfi.End();++vfi)
506 if(tri::IsMarked(m,vfi.F()))
return false;
507 faceSeedVec.push_back(vfi.F());
508 tri::Mark(m,vfi.F());
514 static inline std::string sourcesAttributeName(
void) {
return "sources"; }
515 static inline std::string parentsAttributeName(
void) {
return "parent"; }
517 template <
class DistanceFunctor>
518 static void PerFaceDijkstraCompute(MeshType &m,
const std::vector<FacePointer> &seedVec,
519 DistanceFunctor &distFunc,
520 ScalarType maxDistanceThr = std::numeric_limits<ScalarType>::max(),
521 std::vector<FacePointer> *InInterval=NULL,
522 FacePointer FaceTarget=NULL,
523 bool avoid_selected=
false)
525 tri::RequireFFAdjacency(m);
526 tri::RequirePerFaceMark(m);
527 tri::RequirePerFaceQuality(m);
529 typename MeshType::template PerFaceAttributeHandle<FacePointer> sourceHandle
530 = tri::Allocator<MeshType>::template GetPerFaceAttribute<FacePointer>(m, sourcesAttributeName());
532 typename MeshType::template PerFaceAttributeHandle<FacePointer> parentHandle
533 = tri::Allocator<MeshType>::template GetPerFaceAttribute<FacePointer>(m, parentsAttributeName());
535 std::vector<FaceDist> Heap;
537 for(
size_t i=0;i<seedVec.size();++i)
539 tri::Mark(m,seedVec[i]);
541 sourceHandle[seedVec[i]]=seedVec[i];
542 parentHandle[seedVec[i]]=seedVec[i];
543 Heap.push_back(FaceDist(seedVec[i]));
544 if (InInterval!=NULL) InInterval->push_back(seedVec[i]);
547 std::make_heap(Heap.begin(),Heap.end());
550 pop_heap(Heap.begin(),Heap.end());
551 FacePointer curr = (Heap.back()).f;
552 if ((FaceTarget!=NULL)&&(curr==FaceTarget))
return;
559 FacePointer nextF = curr->FFp(i);
560 ScalarType nextDist = curr->Q() + distFunc(curr,nextF);
561 if( (nextDist < maxDistanceThr) &&
562 (!tri::IsMarked(m,nextF) || nextDist < nextF->Q()) )
564 nextF->Q() = nextDist;
565 if ((avoid_selected)&&(nextF->IsS()))
continue;
567 Heap.push_back(FaceDist(nextF));
568 push_heap(Heap.begin(),Heap.end());
569 if (InInterval!=NULL) InInterval->push_back(nextF);
570 sourceHandle[nextF] = sourceHandle[curr];
571 parentHandle[nextF] = curr;
582 template <
class DistanceFunctor>
583 static void PerVertexDijkstraCompute(MeshType &m,
const std::vector<VertexPointer> &seedVec,
584 DistanceFunctor &distFunc,
585 ScalarType maxDistanceThr = std::numeric_limits<ScalarType>::max(),
586 std::vector<VertexPointer> *InInterval=NULL,
587 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sourceHandle= NULL,
588 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * parentHandle=NULL,
589 bool avoid_selected=
false,
590 VertexPointer target=NULL)
592 tri::RequireVFAdjacency(m);
593 tri::RequirePerVertexMark(m);
594 tri::RequirePerVertexQuality(m);
602 std::vector<DIJKDist> Heap;
606 for(
size_t i=0;i<seedVec.size();++i)
608 assert(!tri::IsMarked(m,seedVec[i]));
609 tri::Mark(m,seedVec[i]);
611 if (sourceHandle!=NULL)
612 (*sourceHandle)[seedVec[i]]=seedVec[i];
613 if (parentHandle!=NULL)
614 (*parentHandle)[seedVec[i]]=seedVec[i];
615 Heap.push_back(DIJKDist(seedVec[i]));
616 if (InInterval!=NULL) InInterval->push_back(seedVec[i]);
619 std::make_heap(Heap.begin(),Heap.end());
622 pop_heap(Heap.begin(),Heap.end());
623 VertexPointer curr = (Heap.back()).v;
624 if ((target!=NULL)&&(target==curr))
return;
626 std::vector<VertexPointer> vertVec;
627 face::VVStarVF<FaceType>(curr,vertVec);
628 for(
size_t i=0;i<vertVec.size();++i)
630 VertexPointer nextV = vertVec[i];
631 if ((avoid_selected)&&(nextV->IsS()))
continue;
632 ScalarType nextDist = curr->Q() + distFunc(curr,nextV);
633 if( (nextDist < maxDistanceThr) &&
634 (!tri::IsMarked(m,nextV) || nextDist < nextV->Q()) )
636 nextV->Q() = nextDist;
638 Heap.push_back(DIJKDist(nextV));
639 push_heap(Heap.begin(),Heap.end());
640 if (InInterval!=NULL) InInterval->push_back(nextV);
641 if (sourceHandle!=NULL)
642 (*sourceHandle)[nextV] = (*sourceHandle)[curr];
643 if (parentHandle!=NULL)
644 (*parentHandle)[nextV] = curr;