vcglib: geodesic.h Source File

00001 /****************************************************************************
00002 * VCGLib                                                            o o     *
00003 * Visual and Computer Graphics Library                            o     o   *
00004 *                                                                _   O  _   *
00005 * Copyright(C) 2004                                                \/)\/    *
00006 * Visual Computing Lab                                            /\/|      *
00007 * ISTI - Italian National Research Council                           |      *
00008 *                                                                    \      *
00009 * All rights reserved.                                                      *
00010 *                                                                           *
00011 * This program is free software; you can redistribute it and/or modify      *   
00012 * it under the terms of the GNU General Public License as published by      *
00013 * the Free Software Foundation; either version 2 of the License, or         *
00014 * (at your option) any later version.                                       *
00015 *                                                                           *
00016 * This program is distributed in the hope that it will be useful,           *
00017 * but WITHOUT ANY WARRANTY; without even the implied warranty of            *
00018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
00019 * GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
00020 * for more details.                                                         *
00021 *                                                                           *
00022 ****************************************************************************/
00023 /*#**************************************************************************
00024   History
00025         $Log: not supported by cvs2svn $
00026         Revision 1.6  2008/01/12 19:07:05  ganovelli
00027         Recompiled from previous out of date version. Still to revise but working
00028         
00029         Revision 1.5  2005/12/13 17:17:19  ganovelli
00030         first importing from old version. NOT optimized! It works with VertexFace Adjacency even over non manifolds
00031         
00032 
00033  *#**************************************************************************/
00034 
00035 
00036 
00037 #include <assert.h>
00038 #include <vcg/math/base.h>
00039 #include <vcg/container/simple_temporary_data.h>
00040 #include <vcg/simplex/face/pos.h>
00041 #include <vcg/simplex/face/topology.h>
00042 #include <vcg/complex/trimesh/update/quality.h>
00043 #include <deque>
00044 #include <vector>
00045 #include <list>
00046 #include <functional>
00047 
00048 /*
00049 class for computing approximated geodesic distances on a mesh.
00050 
00051 basic example: farthest vertex from a specified one
00052         MyMesh m;
00053         MyMesh::VertexPointer seed,far;
00054         MyMesh::ScalarType dist;
00055 
00056   vcg::Geo<MyMesh> g;
00057         g.FarthestVertex(m,seed,far,d);
00058 
00059 */
00060 namespace vcg{
00061 namespace tri{
00062 
00063 template <class MeshType>
00064 struct EuclideanDistance{
00065         typedef typename MeshType::VertexType VertexType;
00066         typedef typename MeshType::ScalarType  ScalarType;
00067 
00068         EuclideanDistance(){}
00069         ScalarType operator()(const VertexType * v0, const VertexType * v1) const
00070         {return vcg::Distance(v0->cP(),v1->cP());}
00071 };
00072 
00073 template <class MeshType, class DistanceFunctor = EuclideanDistance<MeshType> >
00074 class Geo{
00075 
00076         public:
00077 
00078         typedef typename MeshType::VertexType VertexType;
00079         typedef typename MeshType::VertexIterator VertexIterator;
00080         typedef typename MeshType::VertexPointer VertexPointer;
00081         typedef typename MeshType::FaceType  FaceType;
00082         typedef typename MeshType::CoordType  CoordType;
00083         typedef typename MeshType::ScalarType  ScalarType;
00084         
00085 
00086 
00087         /* Auxiliary class for keeping the heap of vertices to visit and their estimated distance
00088         */
00089         struct VertDist{
00090                 VertDist(){}
00091                 VertDist(VertexPointer _v, ScalarType _d):v(_v),d(_d){}
00092                 VertexPointer v;
00093                 ScalarType d;
00094         };
00095 
00096 
00097         /* Temporary data to associate to all the vertices: estimated distance and boolean flag
00098         */
00099         struct TempData{
00100                 TempData(){}
00101                 TempData(const ScalarType & d_){d=d_;source = NULL;}
00102                 ScalarType d;
00103                 VertexPointer source;//closest source
00104 
00105                 };
00106 
00107         typedef SimpleTempData<std::vector<VertexType>, TempData >  TempDataType;
00108         //TempDataType  * TD;
00109         
00110         
00111         struct pred: public std::binary_function<VertDist,VertDist,bool>{
00112                 pred(){};
00113                         bool operator()(const VertDist& v0, const VertDist& v1) const
00114                                 {return (v0.d > v1.d);}
00115                 };
00116         struct pred_addr: public std::binary_function<VertDist,VertDist,bool>{
00117                 pred_addr(){};
00118                         bool operator()(const VertDist& v0, const VertDist& v1) const
00119                                 {return (v0.v > v1.v);}
00120                 };
00121 
00122         //************** calcolo della distanza di pw in base alle distanze note di pw1 e curr
00123         //************** sapendo che (curr,pw,pw1) e'una faccia della mesh
00124         //************** (vedi figura in file distance.gif)
00125         static ScalarType Distance(const VertexPointer &pw,
00126                                                                                    const VertexPointer &pw1,
00127                                                                                    const VertexPointer &curr,
00128                                                                                    const ScalarType &d_pw1,
00129                                                                                    const ScalarType &d_curr)
00130         {
00131                 ScalarType curr_d=0;
00132 
00133                 ScalarType ew_c  = DistanceFunctor()(pw,curr);
00134                 ScalarType ew_w1 = DistanceFunctor()(pw,pw1);
00135                 ScalarType ec_w1 = DistanceFunctor()(pw1,curr);
00136                 CoordType w_c =  (pw->cP()-curr->cP()).Normalize() * ew_c;
00137                 CoordType w_w1 = (pw->cP() - pw1->cP()).Normalize() * ew_w1;
00138                 CoordType w1_c =  (pw1->cP() - curr->cP()).Normalize() * ec_w1;
00139 
00140                 ScalarType      alpha,alpha_, beta,beta_,theta,h,delta,s,a,b;
00141 
00142                 alpha = acos((w_c.dot(w1_c))/(ew_c*ec_w1));
00143                 s = (d_curr + d_pw1+ec_w1)/2;
00144                 a = s/ec_w1;
00145                 b = a*s;
00146                 alpha_ = 2*acos ( math::Min<ScalarType>(1.0,sqrt(  (b- a* d_pw1)/d_curr)));
00147 
00148                 if ( alpha+alpha_ > M_PI){
00149                         curr_d = d_curr + ew_c;         
00150                         }else
00151                         {
00152                                 beta_ = 2*acos ( math::Min<ScalarType>(1.0,sqrt(  (b- a* d_curr)/d_pw1)));
00153                                 beta  = acos((w_w1).dot(-w1_c)/(ew_w1*ec_w1));
00154 
00155                                 if ( beta+beta_ > M_PI)
00156                                         curr_d = d_pw1  + ew_w1;
00157                                 else 
00158                                         {
00159                                         theta   = ScalarType(M_PI)-alpha-alpha_;
00160                                         delta   = cos(theta)* ew_c;
00161                                         h               = sin(theta)* ew_c;
00162                                         curr_d = sqrt( pow(h,2)+ pow(d_curr + delta,2));
00163                                         }
00164                         }
00165                 return (curr_d);
00166         }
00167 
00168         /*
00169           starting from the seeds, it assign a distance value to each vertex. The distance of a vertex is its 
00170           approximated geodesic distance to the closest seeds.
00171           This is function is not meant to be called (although is not prevented). Instead, it is invoked by 
00172           wrapping function.
00173         */
00174         static  VertexPointer Visit( 
00175                  MeshType & m,
00176                  std::vector<VertDist> & _inputfrontier,
00177                  ScalarType & max_distance,
00178                  bool farthestOnBorder = false,
00179                  typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
00180                  )
00181 {
00182         bool isLeaf;
00183         std::vector<VertDist> frontier;
00184         VertexIterator ii;
00185         std::list<VertexPointer> children;
00186   VertexPointer curr,farthest=0,pw1;
00187         typename std::list<VertexPointer>::iterator is;
00188         std::deque<VertexPointer> leaves;
00189         std::vector<VertDist> _frontier;
00190         ScalarType unreached  = std::numeric_limits<ScalarType>::max();
00191 
00192         std::vector <std::pair<VertexPointer,ScalarType> > expansion;
00193         typename std::vector <VertDist >::iterator ifr;
00194   face::VFIterator<FaceType> x;
00195         VertexPointer pw;
00196 
00197         //Requirements
00198         assert(m.HasVFTopology());
00199         assert(!_inputfrontier.empty());
00200 
00201         TempDataType  * TD;
00202         TD = new TempDataType(m.vert,unreached);
00203 
00204         for(ifr = _inputfrontier.begin(); ifr != _inputfrontier.end(); ++ifr){
00205                 (*TD)[(*ifr).v].d = 0.0;        
00206                 (*ifr).d = 0.0;
00207                 (*TD)[(*ifr).v].source  = (*ifr).v;
00208                 frontier.push_back(VertDist((*ifr).v,0.0));
00209                 }
00210                 // initialize Heap
00211                 make_heap(frontier.begin(),frontier.end(),pred());      
00212 
00213                 ScalarType curr_d,d_curr = 0.0,d_heap;
00214                 VertexPointer curr_s = NULL;
00215                 max_distance=0.0;
00216                 typename std::vector<VertDist >:: iterator iv;
00217 
00218                 while(!frontier.empty())
00219                         {  
00220                                 pop_heap(frontier.begin(),frontier.end(),pred());
00221                                 curr = (frontier.back()).v;
00222                                 curr_s = (*TD)[curr].source;
00223                                 if(sources!=NULL) 
00224                                         (*sources)[curr] = curr_s;
00225                                 d_heap = (frontier.back()).d;
00226                                 frontier.pop_back();
00227 
00228                                 assert((*TD)[curr].d <= d_heap); 
00229                                 assert(curr_s != NULL);
00230                                 if((*TD)[curr].d < d_heap )// a vertex whose distance has been improved after it was inserted in the queue
00231                                         continue; 
00232                                 assert((*TD)[curr].d == d_heap); 
00233 
00234                                 d_curr =  (*TD)[curr].d;
00235 
00236                                 isLeaf = (!farthestOnBorder || curr->IsB());
00237 
00238                                 face::VFIterator<FaceType> x;int k;
00239                         
00240                                 for( x.f = curr->VFp(), x.z = curr->VFi(); x.f!=0; ++x )
00241                                         for(k=0;k<2;++k)
00242                                         {
00243                                         if(k==0) {
00244                                                 pw = x.f->V1(x.z);
00245                                                 pw1=x.f->V2(x.z);
00246                                                 }
00247                                         else {
00248                                                 pw = x.f->V2(x.z);
00249                                                 pw1=x.f->V1(x.z);
00250                                                 }
00251 
00252                                                 const ScalarType & d_pw1  =  (*TD)[pw1].d;
00253                                                 {
00254 
00255                                                         const ScalarType inter  = DistanceFunctor()(curr,pw1);//(curr->P() - pw1->P()).Norm();
00256                                                         const ScalarType tol = (inter + d_curr + d_pw1)*.0001f;
00257 
00258                                                         if (    ((*TD)[pw1].source != (*TD)[curr].source)||// not the same source
00259                                                                         (inter + d_curr < d_pw1  +tol   ) ||
00260                                                                         (inter + d_pw1  < d_curr +tol  ) ||
00261                                                                         (d_curr + d_pw1  < inter +tol  )   // triangular inequality
00262                                                                          )  
00263                                                                 curr_d = d_curr + DistanceFunctor()(pw,curr);//(pw->P()-curr->P()).Norm();
00264                                                         else
00265                                                                 curr_d = Distance(pw,pw1,curr,d_pw1,d_curr);
00266                                                         
00267                                                 }
00268                                  
00269                                         if((*TD)[(pw)].d > curr_d){
00270                                                         (*TD)[(pw)].d = curr_d; 
00271                                                         (*TD)[pw].source = curr_s;
00272                                                         frontier.push_back(VertDist(pw,curr_d));
00273                                                         push_heap(frontier.begin(),frontier.end(),pred());
00274                                                 } 
00275                                         if(isLeaf){
00276                                                 if(d_curr > max_distance){
00277                                                         max_distance = d_curr;
00278                                                         farthest = curr;
00279                                                         }
00280                                                 }
00281                                         }
00282         }// end while
00283 
00284         // scrivi le distanze sul campo qualita' (nn: farlo parametrico)
00285         VertexIterator vi;
00286         for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
00287                 (*vi).Q() =  (*TD)[&(*vi)].d; 
00288 
00289         delete TD;
00290   assert(farthest);
00291         return farthest;
00292 
00293  }
00294         
00295 
00296 public:
00297         /*
00298         Given a mesh and  a vector of pointers to vertices (sources), assigns the approximated geodesic
00299         distance from the cloasest source to all the mesh vertices and returns the pointer to the farthest.
00300         Note: update the field Q() of the vertices
00301         */
00302  static bool FarthestVertex( MeshType & m,
00303                                                                         std::vector<VertexPointer> & fro, 
00304                                                                         VertexPointer & farthest,                
00305                                                                         ScalarType & distance,
00306                                                                         typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL){                                  
00307 
00308                                                                         typename std::vector<VertexPointer>::iterator fi; 
00309                                                                         std::vector<VertDist>fr;
00310                                                                         if(fro.empty()) return false;
00311                                                                                 
00312                                                                         for( fi  = fro.begin(); fi != fro.end() ; ++fi)
00313                                                                                 fr.push_back(VertDist(*fi,0.0));
00314                                                                         farthest = Visit(m,fr,distance,false,sources); 
00315                                                                         return true;
00316                                           }
00317         /*
00318         Given a mesh and  a  pointers to a vertex-source (source), assigns the approximated geodesic
00319         distance from the vertex-source to all the mesh vertices and returns the pointer to the farthest
00320         Note: update the field Q() of the vertices
00321         */
00322  static void FarthestVertex( MeshType & m, 
00323                                                                         VertexPointer seed,
00324                                                                         VertexPointer & farthest,                
00325                                                                         ScalarType & distance){ 
00326         std::vector<VertexPointer>  fro;
00327         fro.push_back( seed );
00328         VertexPointer v0;
00329         FarthestVertex(m,fro,v0,distance);
00330         farthest = v0;
00331 }
00332 
00333 /* 
00334         Same as FarthestPoint but the returned pointer is to a border vertex
00335         Note: update the field Q() of the vertices
00336 */
00337  static void FarthestBVertex(MeshType & m,
00338                                                                                 std::vector<VertexPointer> & fro,  
00339                                                                                 VertexPointer & farthest,            
00340                                                                                 ScalarType & distance,
00341                                                                                 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
00342                                                                                 ){
00343 
00344         typename std::vector<VertexPointer>::iterator fi; 
00345         std::vector<VertDist>fr;
00346 
00347         for( fi  = fro.begin(); fi != fro.end() ; ++fi)
00348                 fr.push_back(VertDist(*fi,-1));
00349         farthest =  Visit(m,fr,distance,true,sources); 
00350 }
00351 /* 
00352         Same as FarthestPoint but the returned pointer is to a border vertex
00353         Note: update the field Q() of the vertices
00354 */
00355  static void FarthestBVertex(   MeshType & m, 
00356                                                                 VertexPointer seed,
00357                                                                 VertexPointer & farthest,                
00358                                                                 ScalarType & distance,
00359                                                                 typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL){  
00360         std::vector<VertexPointer>  fro;
00361         fro.push_back( seed );
00362         VertexPointer v0;
00363         FarthestBVertex(m,fro,v0,distance,sources);
00364         farthest = v0;
00365  }
00366 
00367 /* 
00368         Assigns to each vertex of the mesh its distance to the closest vertex on the border
00369         Note: update the field Q() of the vertices
00370 */
00371  static bool DistanceFromBorder(        MeshType & m,
00372                                                                         ScalarType & distance,
00373                                                                         typename MeshType::template PerVertexAttributeHandle<VertexPointer> * sources = NULL
00374                                         ){
00375         std::vector<VertexPointer> fro;
00376         VertexIterator vi;
00377         VertexPointer farthest;
00378         for(vi = m.vert.begin(); vi != m.vert.end(); ++vi)
00379                 if( (*vi).IsB())
00380                         fro.push_back(&(*vi));
00381         if(fro.empty()) return false;
00382         
00383         tri::UpdateQuality<MeshType>::VertexConstant(m,0);
00384                 
00385         return FarthestVertex(m,fro,farthest,distance,sources);
00386 }
00387 
00388  };
00389 };// end namespace tri
00390 };// end namespace vcg