vcglib: distance.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 
00026 $Log: not supported by cvs2svn $
00027 Revision 1.10  2006/01/22 10:06:23  cignoni
00028 Corrected use of Area with the unambiguous DoubleArea
00029 
00030 Revision 1.9  2005/09/28 19:35:06  m_di_benedetto
00031 Added class PointDistanceFunctor.
00032 
00033 Revision 1.8  2005/09/14 12:58:44  pietroni
00034 changed min calls to Min<ScalarType> of math.h of vcglib
00035 
00036 Revision 1.7  2005/09/14 09:58:32  pietroni
00037 removed vcg::math::Min<ScalarType> definition generate warnings
00038 
00039 Revision 1.6  2005/09/14 09:03:54  pietroni
00040 added definition of vcg::math::Min<ScalarType> function
00041 
00042 Revision 1.5  2005/02/02 16:44:34  pietroni
00043 1 warning corrected added casting in const ScalarType EPSILON = ScalarType( 0.000001);
00044 
00045 Revision 1.4  2005/01/28 12:00:33  cignoni
00046 small gcc compiling issues for namespaces
00047 
00048 Revision 1.3  2005/01/24 15:35:25  cignoni
00049 Removed a 'using namespace'
00050 
00051 Revision 1.2  2005/01/21 17:11:03  pietroni
00052 changed Dist Function to PointDistance... the function is on vcg::face::PointDistance this file will contain all distance functions between a face and othe entities
00053 
00054 Revision 1.1  2004/05/12 18:50:25  ganovelli
00055 created
00056 
00057 
00058 ****************************************************************************/
00059 
00060 #ifndef __VCGLIB_FACE_DISTANCE
00061 #define __VCGLIB_FACE_DISTANCE
00062 
00063 #include <vcg/math/base.h>
00064 #include <vcg/space/point3.h>
00065 #include <vcg/space/segment3.h>
00066 
00067 
00068 namespace vcg {
00069         namespace face{
00070 /*
00071    Point face distance
00072    trova il punto <p> sulla faccia piu' vicino a <q>, con possibilit� di 
00073    rejection veloce su se la distanza trovata � maggiore di <rejdist>
00074 
00075  Commenti del 12/11/02
00076  Funziona solo se la faccia e di quelle di tipo E (con edge e piano per faccia gia' calcolati)
00077  algoritmo:
00078         1) si calcola la proiezione <p> di q sul piano della faccia
00079         2) se la distanza punto piano e' > rejdist ritorna
00080         3) si lavora sul piano migliore e si cerca di capire se il punto sta dentro il triangolo:
00081            a) prodotto vettore tra edge triangolo (v[i+1]-v[i]) e (p-v[i])
00082                  b) se il risultato e' negativo (gira in senso orario) allora il punto
00083                     sta fuori da quella parte e si fa la distanza punto segmento.
00084      c) se il risultato sempre positivo allora sta dentro il triangolo
00085         4) e si restituisce la distanza punto /piano gia` calcolata 
00086 
00087         Note sulla robustezza:
00088         il calcolo del prodotto vettore e` la cosa piu` delicata:
00089         possibili fallimenti quando a^b ~= 0
00090         1) doveva essere <= 0 e viene positivo (q era fuori o sulla linea dell'edge)
00091            allora capita che si faccia la distanza punto piano anziche` la distanza punto seg
00092   2) doveva essere > 0 e viene <=0 (q era dentro il triangolo)
00093 
00094 */
00095         template <class FaceType>
00096         bool PointDistance(     const FaceType &f, 
00097                                                         const vcg::Point3<typename FaceType::ScalarType> & q, 
00098                                                         typename FaceType::ScalarType & dist, 
00099                                                         vcg::Point3<typename FaceType::ScalarType> & p )
00100         {
00101                 typedef typename FaceType::ScalarType ScalarType;
00102                 
00103                 const ScalarType EPS = ScalarType( 0.000001);
00104 
00105                 //const ScalarType EPSILON = 0.00000001;
00106                 ScalarType b,b0,b1,b2;
00107                         // Calcolo distanza punto piano
00108                 ScalarType d = Distance( f.cPlane(), q );
00109                 if( d>dist || d<-dist )                 // Risultato peggiore: niente di fatto
00110                         return false;
00111 
00112                         // Calcolo del punto sul piano
00113                 // NOTA: aggiunto un '-d' in fondo Paolo C.
00114                 Point3<ScalarType> t = f.cPlane().Direction();
00115                 t[0] *= -d;
00116                 t[1] *= -d;
00117                 t[2] *= -d;
00118                 p = q; p += t;
00119                     
00120                 switch( f.Flags() & (FaceType::NORMX|FaceType::NORMY|FaceType::NORMZ) )
00121                 {
00122                 case FaceType::NORMX:
00123                         b0 = f.cEdge(1)[1]*(p[2] - f.cP(1)[2]) - f.cEdge(1)[2]*(p[1] - f.cP(1)[1]);
00124                         if(b0<=0)
00125                         {
00126                                 b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00127                                 if(dist>b0) { dist = b0; return true; }
00128                                 else return false;
00129                         }
00130                         b1 = f.cEdge(2)[1]*(p[2] - f.cP(2)[2]) - f.cEdge(2)[2]*(p[1] - f.cP(2)[1]);
00131                         if(b1<=0)
00132                         {
00133                                 b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00134                                 if(dist>b1) { dist = b1; return true; }
00135                                 else return false;
00136                         }
00137                         b2 = f.cEdge(0)[1]*(p[2] - f.cP(0)[2]) - f.cEdge(0)[2]*(p[1] - f.cP(0)[1]);
00138                         if(b2<=0)
00139                         {
00140                                 b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00141                                 if(dist>b2) { dist = b2; return true; }
00142                                 else return false;
00143                         }
00144                         // sono tutti e tre > 0 quindi dovrebbe essere dentro;
00145                         // per sicurezza se il piu' piccolo dei tre e' < epsilon (scalato rispetto all'area della faccia
00146                         // per renderlo dimension independent.) allora si usa ancora la distanza punto 
00147                         // segmento che e' piu robusta della punto piano, e si fa dalla parte a cui siamo piu' 
00148                         // vicini (come prodotto vettore)
00149                         // Nota: si potrebbe rendere un pochino piu' veloce sostituendo Area()
00150                         // con il prodotto vettore dei due edge in 2d lungo il piano migliore.
00151       if( (b=std::min(b0,std::min(b1,b2)) ) < EPS*DoubleArea(f))
00152       {
00153                                 ScalarType bt;
00154                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00155                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00156         else {          assert(b==b2);
00157                         bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00158                       }
00159                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPS*Area(),bt);
00160                                 if(dist>bt) { dist = bt; return true; }
00161                                 else return false;
00162                         }
00163                         break;
00164 
00165                   case FaceType::NORMY:
00166                         b0 = f.cEdge(1)[2]*(p[0] - f.cP(1)[0]) - f.cEdge(1)[0]*(p[2] - f.cP(1)[2]);
00167                         if(b0<=0)
00168                         {
00169                                 b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00170                                 if(dist>b0) { dist = b0; return true; }
00171                                 else return false;
00172                         }
00173                         b1 = f.cEdge(2)[2]*(p[0] - f.cP(2)[0]) - f.cEdge(2)[0]*(p[2] - f.cP(2)[2]);
00174                         if(b1<=0)
00175                         {
00176                                 b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00177                                 if(dist>b1) { dist = b1; return true; }
00178                                 else return false;
00179                         }
00180                         b2 = f.cEdge(0)[2]*(p[0] - f.cP(0)[0]) - f.cEdge(0)[0]*(p[2] - f.cP(0)[2]);
00181                         if(b2<=0)
00182                         {
00183                                 b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00184                                 if(dist>b2) { dist = b2; return true; }
00185                                 else return false;
00186                         }
00187                         if( (b=vcg::math::Min<ScalarType>(b0,vcg::math::Min<ScalarType>(b1,b2))) < EPS*DoubleArea(f))
00188       {
00189                                 ScalarType bt;
00190                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00191                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00192         else { assert(b==b2);
00193                         bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00194         }
00195                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
00196                                 if(dist>bt) { dist = bt; return true; }
00197                                 else return false;
00198                         }
00199                         break;
00200 
00201                   case FaceType::NORMZ:
00202                         b0 = f.cEdge(1)[0]*(p[1] - f.cP(1)[1]) - f.cEdge(1)[1]*(p[0] - f.cP(1)[0]);
00203                         if(b0<=0)
00204                         {
00205                                 b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00206                                 if(dist>b0) { dist = b0; return true; }
00207                                 else return false;
00208                         }
00209                         b1 = f.cEdge(2)[0]*(p[1] - f.cP(2)[1]) - f.cEdge(2)[1]*(p[0] - f.cP(2)[0]);
00210                         if(b1<=0)
00211                         {
00212                                 b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00213                                 if(dist>b1) { dist = b1; return true; }
00214                                 else return false;
00215                         }
00216                         b2 = f.cEdge(0)[0]*(p[1] - f.cP(0)[1]) - f.cEdge(0)[1]*(p[0] - f.cP(0)[0]);
00217                         if(b2<=0)
00218                         {
00219                                 b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00220                                 if(dist>b2) { dist = b2; return true; }
00221                                 else return false;
00222                         }
00223                         if( (b=vcg::math::Min<ScalarType>(b0,vcg::math::Min<ScalarType>(b1,b2))) < EPS*DoubleArea(f))
00224       {
00225                                 ScalarType bt;
00226                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00227                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00228         else { assert(b==b2);
00229                         bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00230         }
00231                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
00232                                 
00233                                 if(dist>bt) { dist = bt; return true; }
00234                                 else return false;
00235                         }
00236                         break;
00237 
00238                 }
00239 
00240                 dist = ScalarType(fabs(d));
00241                 //dist = Distance(p,q);
00242                 return true;
00243         }
00244 
00245         template <class S>
00246         class PointDistanceFunctor {
00247         public:
00248                 typedef S ScalarType;
00249                 typedef Point3<ScalarType> QueryType;
00250                 static inline const Point3<ScalarType> &  Pos(const QueryType & qt)  {return qt;}
00251 
00252                 template <class FACETYPE, class SCALARTYPE>
00253                 inline bool operator () (const FACETYPE & f, const Point3<SCALARTYPE> & p, SCALARTYPE & minDist, Point3<SCALARTYPE> & q) {
00254                         const Point3<typename FACETYPE::ScalarType> fp = Point3<typename FACETYPE::ScalarType>::Construct(p);
00255                         Point3<typename FACETYPE::ScalarType> fq;
00256                         typename FACETYPE::ScalarType md = (typename FACETYPE::ScalarType)(minDist);
00257                         const bool ret = PointDistance(f, fp, md, fq);
00258                         minDist = (SCALARTYPE)(md);
00259                         q = Point3<SCALARTYPE>::Construct(fq);
00260                         return (ret);
00261                 }
00262         };
00263 
00264 
00265 
00266         template <class S>
00267         class PointNormalDistanceFunctor {
00268         public:
00269                 typedef typename S::ScalarType ScalarType;
00270                 typedef S QueryType;
00271                 static inline const Point3<ScalarType> &  Pos(const QueryType & qt)  {return qt.P();}
00272 
00273 
00274                 static ScalarType & Alpha(){static ScalarType alpha = 1.0; return alpha;}
00275                 static ScalarType & Beta (){static ScalarType beta  = 1.0; return beta;}
00276                 static ScalarType & Gamma(){static ScalarType gamma = 1.0; return gamma;}
00277                 static ScalarType & InterPoint(){static ScalarType interpoint= 1.0; return interpoint;} 
00278 
00279 
00280                 template <class FACETYPE, class SCALARTYPE>
00281                 inline bool operator () (const FACETYPE &f, const typename FACETYPE::VertexType &p, 
00282                         SCALARTYPE & minDist,Point3<SCALARTYPE> & q)
00283                 {               
00284                         const Point3<typename FACETYPE::ScalarType> fp = Point3<typename FACETYPE::ScalarType>::Construct(p.cP());
00285                         const Point3<typename FACETYPE::ScalarType> fn = Point3<typename FACETYPE::ScalarType>::Construct(p.cN());
00286                         Point3<typename FACETYPE::ScalarType> fq;
00287                         typename FACETYPE::ScalarType md = (typename FACETYPE::ScalarType)(minDist);
00288                         const bool ret=PointDistance(f,fp,md,fq);
00289 
00290                         SCALARTYPE  dev=InterPoint()*(pow((ScalarType)(1-f.cN().dot(fn)),(ScalarType)Beta())/(Gamma()*md+0.1));
00291 
00292                         if (md+dev < minDist){
00293                                 minDist = (SCALARTYPE)(md+dev);
00294                                 q = Point3<SCALARTYPE>::Construct(fq);
00295                                 //q.N() = f.N();
00296                                 return (ret);
00297                         }
00298                         return false;
00299                 }
00300         };
00301                 
00305                 // UpdateNormals::PerFaceNormalized(m)
00306                 // UpdateFlags<>::FaceProjection(m);
00307                 
00308                 template <class FaceType>
00309                         bool PointDistanceBase(
00310                                                                                                         const FaceType &f,                                                                                                                                              
00311                                                                                                         const vcg::Point3<typename FaceType::ScalarType> & q, 
00312                                                                                                         typename FaceType::ScalarType & dist,                 
00313                                                                                                         vcg::Point3<typename FaceType::ScalarType> & p )      
00314                 {
00315                                 typedef typename FaceType::ScalarType ScalarType;
00316                                 // remember that the macro NDEBUG is defined when you want to optimize a lot. 
00317                                 #ifndef NDEBUG
00318                                 static int staticCnt=0; // small piece of code that sometime check that face normals are really normalized
00319                                 if((staticCnt++%100)==0) 
00320             assert((f.cN().SquaredNorm() ==0) || (f.cN().SquaredNorm() > 0.9999 && f.cN().SquaredNorm()<1.0001)); // if you get this assert you have forgot to make a UpdateNormals::PerFaceNormalized(m)
00321         #endif
00322 
00323         if(f.cN()==Point3<ScalarType>(0,0,0)) // to correctly manage the case of degenerate triangles we consider them as segments.
00324         {
00325             Box3<ScalarType> bb;
00326             f.GetBBox(bb);
00327             Segment3<ScalarType> degenTri(bb.min,bb.max);
00328             Point3<ScalarType> closest= ClosestPoint( degenTri, q );
00329             ScalarType d = Distance(closest, q);
00330             if( d>dist || d<-dist )                     // Risultato peggiore: niente di fatto
00331                       return false;
00332             dist=d;
00333             p=closest;
00334             return true;
00335         }
00336 
00337                                 Plane3<ScalarType> fPlane;
00338                                 fPlane.Init(f.cP(0),f.cN());
00339         const ScalarType EPS = ScalarType( 0.000001);
00340                                 ScalarType b,b0,b1,b2;
00341                                 // Calcolo distanza punto piano
00342                                 ScalarType d = Distance( fPlane, q );
00343                                 if( d>dist || d<-dist )                 // Risultato peggiore: niente di fatto
00344                                         return false;
00345                                 
00346                                 // Calcolo del punto sul piano
00347                                 // NOTA: aggiunto un '-d' in fondo Paolo C.
00348                                 Point3<ScalarType> t = fPlane.Direction();
00349                                 t[0] *= -d;
00350                                 t[1] *= -d;
00351                                 t[2] *= -d;
00352                                 p = q; p += t;
00353 
00354                                 Point3<ScalarType> fEdge[3];                            
00355                                 fEdge[0] = f.cP(1); fEdge[0] -= f.cP(0);
00356                                 fEdge[1] = f.cP(2); fEdge[1] -= f.cP(1);
00357                                 fEdge[2] = f.cP(0); fEdge[2] -= f.cP(2);
00358                                 
00359                                 /* 
00360                                 This piece of code is part of the EdgePlane initialization structure: note that the edges are scaled!. 
00361                                  
00362                                 if(nx>ny && nx>nz) { f.Flags() |= FaceType::NORMX; d = 1/f.Plane().Direction()[0]; }
00363                                 else if(ny>nz)     { f.Flags() |= FaceType::NORMY; d = 1/f.Plane().Direction()[1]; }
00364                                 else               { f.Flags() |= FaceType::NORMZ; d = 1/f.Plane().Direction()[2]; }
00365                                 f.Edge(0)*=d; f.Edge(1)*=d;f.Edge(2)*=d;
00366                                  
00367                                 So we must apply the same scaling according to the plane orientation, eg in the case of NORMX
00368                                 
00369                                 scaleFactor= 1/fPlane.Direction()[0];
00370                                 fEdge[0]*=d; fEdge[1]*=d;fEdge[2]*=d;
00371                                 */
00372                                 
00373                                 ScalarType scaleFactor;
00374 
00375                                 switch( f.Flags() & (FaceType::NORMX|FaceType::NORMY|FaceType::NORMZ) )
00376                                 {
00377                                         case FaceType::NORMX:
00378                                                 scaleFactor= 1/fPlane.Direction()[0];
00379                                                 fEdge[0]*=scaleFactor; fEdge[1]*=scaleFactor; fEdge[2]*=scaleFactor;
00380                                                 
00381                                                 b0 = fEdge[1][1]*(p[2] - f.cP(1)[2]) - fEdge[1][2]*(p[1] - f.cP(1)[1]);
00382                                                 if(b0<=0)
00383                                                 {
00384                                                         b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00385                                                         if(dist>b0) { dist = b0; return true; }
00386                                                         else return false;
00387                                                 }
00388                                                         b1 = fEdge[2][1]*(p[2] - f.cP(2)[2]) - fEdge[2][2]*(p[1] - f.cP(2)[1]);
00389                                                 if(b1<=0)
00390                                                 {
00391                                                         b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00392                                                         if(dist>b1) { dist = b1; return true; }
00393                                                         else return false;
00394                                                 }
00395                                                         b2 = fEdge[0][1]*(p[2] - f.cP(0)[2]) - fEdge[0][2]*(p[1] - f.cP(0)[1]);
00396                                                 if(b2<=0)
00397                                                 {
00398                                                         b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00399                                                         if(dist>b2) { dist = b2; return true; }
00400                                                         else return false;
00401                                                 }
00402                                                         // sono tutti e tre > 0 quindi dovrebbe essere dentro;
00403                                                         // per sicurezza se il piu' piccolo dei tre e' < epsilon (scalato rispetto all'area della faccia
00404                                                         // per renderlo dimension independent.) allora si usa ancora la distanza punto 
00405                                                         // segmento che e' piu robusta della punto piano, e si fa dalla parte a cui siamo piu' 
00406                                                         // vicini (come prodotto vettore)
00407                                                         // Nota: si potrebbe rendere un pochino piu' veloce sostituendo Area()
00408                                                         // con il prodotto vettore dei due edge in 2d lungo il piano migliore.
00409               if( (b=vcg::math::Min<ScalarType>(b0,vcg::math::Min<ScalarType>(b1,b2))) < EPS*DoubleArea(f))
00410                                                         {
00411                                                                 ScalarType bt;
00412                                                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00413                                                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00414                                                                 else if(b==b2)  bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00415                                                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
00416                                                                 if(dist>bt) { dist = bt; return true; }
00417                                                                 else return false;
00418                                                         }
00419                                                         break;
00420                                                 
00421                                         case FaceType::NORMY:
00422                                                 scaleFactor= 1/fPlane.Direction()[1];
00423                                                 fEdge[0]*=scaleFactor; fEdge[1]*=scaleFactor; fEdge[2]*=scaleFactor;
00424                                                 
00425                                                 b0 = fEdge[1][2]*(p[0] - f.cP(1)[0]) - fEdge[1][0]*(p[2] - f.cP(1)[2]);
00426                                                 if(b0<=0)
00427                                                 {
00428                                                         b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00429                                                         if(dist>b0) { dist = b0; return true; }
00430                                                         else return false;
00431                                                 }
00432                                                         b1 = fEdge[2][2]*(p[0] - f.cP(2)[0]) - fEdge[2][0]*(p[2] - f.cP(2)[2]);
00433                                                 if(b1<=0)
00434                                                 {
00435                                                         b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00436                                                         if(dist>b1) { dist = b1; return true; }
00437                                                         else return false;
00438                                                 }
00439                                                         b2 = fEdge[0][2]*(p[0] - f.cP(0)[0]) - fEdge[0][0]*(p[2] - f.cP(0)[2]);
00440                                                 if(b2<=0)
00441                                                 {
00442                                                         b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00443                                                         if(dist>b2) { dist = b2; return true; }
00444                                                         else return false;
00445                                                 }
00446             if( (b=vcg::math::Min<ScalarType>(b0,vcg::math::Min<ScalarType>(b1,b2))) < EPS*DoubleArea(f))
00447                                                         {
00448                                                                 ScalarType bt;
00449                                                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00450                                                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00451                                                                 else if(b==b2)  bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00452                                                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
00453                                                                 if(dist>bt) { dist = bt; return true; }
00454                                                                 else return false;
00455                                                         }
00456                                                         break;
00457                                                 
00458                                         case FaceType::NORMZ:
00459                                                 scaleFactor= 1/fPlane.Direction()[2];
00460                                                 fEdge[0]*=scaleFactor; fEdge[1]*=scaleFactor; fEdge[2]*=scaleFactor;
00461                                                 
00462                                                 b0 = fEdge[1][0]*(p[1] - f.cP(1)[1]) - fEdge[1][1]*(p[0] - f.cP(1)[0]);
00463                                                 if(b0<=0)
00464                                                 {
00465                                                         b0 = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00466                                                         if(dist>b0) { dist = b0; return true; }
00467                                                         else return false;
00468                                                 }
00469                                                         b1 = fEdge[2][0]*(p[1] - f.cP(2)[1]) - fEdge[2][1]*(p[0] - f.cP(2)[0]);
00470                                                 if(b1<=0)
00471                                                 {
00472                                                         b1 = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00473                                                         if(dist>b1) { dist = b1; return true; }
00474                                                         else return false;
00475                                                 }
00476                                                         b2 = fEdge[0][0]*(p[1] - f.cP(0)[1]) - fEdge[0][1]*(p[0] - f.cP(0)[0]);
00477                                                 if(b2<=0)
00478                                                 {
00479                                                         b2 = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00480                                                         if(dist>b2) { dist = b2; return true; }
00481                                                         else return false;
00482                                                 }
00483             if( (b=vcg::math::Min<ScalarType>(b0,vcg::math::Min<ScalarType>(b1,b2))) < EPS*DoubleArea(f))
00484                                                         {
00485                                                                 ScalarType bt;
00486                                                                 if(b==b0)           bt = PSDist(q,f.V(1)->cP(),f.V(2)->cP(),p);
00487                                                                 else if(b==b1)  bt = PSDist(q,f.V(2)->cP(),f.V(0)->cP(),p);
00488                                                                 else if(b==b2)  bt = PSDist(q,f.V(0)->cP(),f.V(1)->cP(),p);
00489                                                                 //printf("Warning area:%g %g %g %g thr:%g bt:%g\n",Area(), b0,b1,b2,EPSILON*Area(),bt);
00490                                                                 
00491                                                                 if(dist>bt) { dist = bt; return true; }
00492                                                                 else return false;
00493                                                         }
00494                                                         break;
00495                                                 default: assert(0); // if you get this assert it means that you forgot to set the required UpdateFlags<MeshType>::FaceProjection(m);
00496 
00497                                 }
00498                                 
00499                                 dist = ScalarType(fabs(d));
00500                                 //dist = Distance(p,q);
00501                                 return true;
00502                 }
00503                 
00504         template <class S>
00505         class PointDistanceBaseFunctor {
00506 public:
00507                         typedef S ScalarType;
00508                         typedef Point3<ScalarType> QueryType;
00509 
00510                   static inline const Point3<ScalarType> & Pos(const Point3<ScalarType> & qt)  {return qt;}
00511                         template <class FACETYPE, class SCALARTYPE>
00512                         inline bool operator () (const FACETYPE & f, const Point3<SCALARTYPE> & p, SCALARTYPE & minDist, Point3<SCALARTYPE> & q) {
00513                                 const Point3<typename FACETYPE::ScalarType> fp = Point3<typename FACETYPE::ScalarType>::Construct(p);
00514                                 Point3<typename FACETYPE::ScalarType> fq;
00515                                 typename FACETYPE::ScalarType md = (typename FACETYPE::ScalarType)(minDist);
00516                                 const bool ret = PointDistanceBase(f, fp, md, fq);
00517                                 minDist = (SCALARTYPE)(md);
00518                                 q = Point3<SCALARTYPE>::Construct(fq);
00519                                 return (ret);
00520                         }
00521                 };
00522                 
00523                 
00524                 
00525 }        // end namespace face
00526         
00527 }        // end namespace vcg
00528 
00529 
00530 #endif
00531