vcglib: smooth.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.19  2008/05/08 23:50:44  cignoni
00027 renamed vertex quality smoothing
00028 added face normal smoothing FF (and added a VF to the previous face normal smoothing)
00029 
00030 Revision 1.18  2008/05/02 09:43:25  cignoni
00031 Added color smoothing, scale dependent laplacian changed a SD_old into SD fujumori, improved comments.
00032 
00033 Revision 1.17  2008/04/18 17:48:29  cignoni
00034 added facenormal smoothing
00035 
00036 Revision 1.16  2008/02/07 10:24:51  cignoni
00037 added a missing IsD() check
00038 
00039 Revision 1.15  2007/11/05 23:47:20  cignoni
00040 added selection to the pasodoble smoothing
00041 
00042 Revision 1.14  2007/03/27 09:40:47  cignoni
00043 Changed use of selected to visited flags. Improved variable namings and comments
00044 
00045 Revision 1.13  2006/11/07 15:13:56  zifnab1974
00046 Necessary changes for compilation with gcc 3.4.6. Especially the hash function is a problem
00047 
00048 Revision 1.12  2006/11/07 11:28:02  cignoni
00049 Added Quality weighted laplacian smoothing
00050 
00051 Revision 1.11  2006/10/19 07:33:03  cignoni
00052 Corrected Laplacian, Added selection to HCSmooth
00053 
00054 Revision 1.10  2006/09/25 09:41:41  cignoni
00055 Added new version of pasodoble smoothing
00056 
00057 Revision 1.9  2006/02/06 10:45:47  cignoni
00058 Added missing typenames
00059 
00060 Revision 1.7  2006/01/24 13:23:22  pietroni
00061 used template types instead of point3f and float inside function calls
00062 
00063 Revision 1.6  2005/12/06 17:55:16  pietroni
00064 1 bug corrected
00065 
00066 Revision 1.5  2005/12/02 16:24:56  pietroni
00067 corrected 1 bug in Cross Prod Gradient
00068 
00069 Revision 1.4  2005/11/23 16:24:44  pietroni
00070 corrected CrossProdGradient( )
00071 
00072 Revision 1.3  2005/07/11 13:12:05  cignoni
00073 small gcc-related compiling issues (typenames,ending cr, initialization order)
00074 
00075 Revision 1.2  2005/03/16 16:14:12  spinelli
00076 aggiunta funzione PasoDobleSmooth e relative:
00077 
00078 - FitMesh
00079 - FaceErrorGrad
00080 - CrossProdGradient
00081 - TriAreaGradient
00082 - NormalSmooth
00083 
00084 e le classi:
00085 
00086 - PDVertInfo
00087 - PDFaceInfo
00088 
00089 necessarie per utilizzare SimpleTempData
00090 
00091 Revision 1.1  2004/12/11 14:53:19  ganovelli
00092 first partial porting: compiled gcc,intel and msvc
00093 
00094 
00095 ****************************************************************************/
00096 
00097 
00098 #ifndef __VCGLIB__SMOOTH
00099 #define __VCGLIB__SMOOTH
00100 
00101 #include <wrap/callback.h>
00102 #include <vcg/space/point3.h>
00103 #include <vcg/space/line3.h>
00104 #include <vcg/container/simple_temporary_data.h>
00105 #include <vcg/complex/trimesh/update/normal.h>
00106 
00107 namespace vcg
00108 {
00109 namespace tri
00110 {
00112 
00114 
00115 
00116 template <class SmoothMeshType>
00117 class Smooth
00118 {
00119 
00120 public:
00121                         typedef SmoothMeshType MeshType;
00122                         typedef typename MeshType::VertexType     VertexType;
00123                         typedef typename MeshType::VertexType::CoordType     CoordType;
00124                         typedef typename MeshType::VertexPointer  VertexPointer;
00125                         typedef typename MeshType::VertexIterator VertexIterator;
00126                         typedef typename MeshType::ScalarType                   ScalarType;
00127                         typedef typename MeshType::FaceType       FaceType;
00128                         typedef typename MeshType::FacePointer    FacePointer;
00129                         typedef typename MeshType::FaceIterator   FaceIterator;
00130                         typedef typename MeshType::FaceContainer  FaceContainer;
00131       typedef typename vcg::Box3<ScalarType>  Box3Type;
00132         typedef typename vcg::face::VFIterator<FaceType> VFLocalIterator;
00133 
00134 class ScaleLaplacianInfo
00135 {
00136 public:
00137         CoordType PntSum;
00138         ScalarType LenSum;
00139 };
00140 
00141 // This is precisely what curvature flow does.
00142 // Curvature flow smoothes the surface by moving along the surface
00143 // normal n with a speed equal to the mean curvature
00144 void VertexCoordLaplacianCurvatureFlow(MeshType &m, int step, ScalarType delta)
00145 {
00146 
00147 }
00148 
00149 // Another Laplacian smoothing variant,
00150 // here we sum the baricenter of the faces incidents on each vertex weighting them with the angle
00151 
00152 static void VertexCoordLaplacianAngleWeighted(MeshType &m, int step, ScalarType delta)
00153 {
00154         ScaleLaplacianInfo lpz;
00155         lpz.PntSum=CoordType(0,0,0);
00156         lpz.LenSum=0;
00157         SimpleTempData<typename MeshType::VertContainer, ScaleLaplacianInfo > TD(m.vert,lpz);
00158         FaceIterator fi;
00159         for(int i=0;i<step;++i)
00160         {
00161                 VertexIterator vi;
00162                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00163                          TD[*vi]=lpz;
00164                 ScalarType a[3];
00165                 for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
00166                 {
00167                         CoordType  mp=((*fi).V(0)->P() + (*fi).V(1)->P() + (*fi).V(2)->P())/3.0;
00168                         CoordType  e0=((*fi).V(0)->P() - (*fi).V(1)->P()).Normalize();
00169                         CoordType  e1=((*fi).V(1)->P() - (*fi).V(2)->P()).Normalize();
00170                         CoordType  e2=((*fi).V(2)->P() - (*fi).V(0)->P()).Normalize();
00171 
00172                         a[0]=AngleN(-e0,e2);
00173                         a[1]=AngleN(-e1,e0);
00174                         a[2]=AngleN(-e2,e1);
00175                         //assert(fabs(M_PI -a[0] -a[1] -a[2])<0.0000001);
00176 
00177                         for(int j=0;j<3;++j){
00178                                                 CoordType dir= (mp-(*fi).V(j)->P()).Normalize();
00179                                                 TD[(*fi).V(j)].PntSum+=dir*a[j];
00180                                                 TD[(*fi).V(j)].LenSum+=a[j]; // well, it should be named angleSum
00181                         }
00182                 }
00183                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00184                                 if(!(*vi).IsD() && TD[*vi].LenSum>0 )
00185                                  (*vi).P() = (*vi).P() +  (TD[*vi].PntSum/TD[*vi].LenSum ) * delta;
00186 
00187         }
00188 };
00189 
00190 // Scale dependent laplacian smoothing [Fujiwara 95]
00191 // as described in
00192 // Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow
00193 // Mathieu Desbrun, Mark Meyer, Peter Schroeder, Alan H. Barr
00194 // SIGGRAPH 99
00195 // REQUIREMENTS: Border Flags.
00196 //
00197 // Note the delta parameter is in a absolute unit
00198 // to get stability it should be a small percentage of the shortest edge.
00199 
00200 static void VertexCoordScaleDependentLaplacian_Fujiwara(MeshType &m, int step, ScalarType delta)
00201 {
00202         SimpleTempData<typename MeshType::VertContainer, ScaleLaplacianInfo > TD(m.vert);
00203         ScaleLaplacianInfo lpz;
00204         lpz.PntSum=CoordType(0,0,0);
00205         lpz.LenSum=0;
00206         FaceIterator fi;
00207         for(int i=0;i<step;++i)
00208         {
00209                         VertexIterator vi;
00210                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00211                                  TD[*vi]=lpz;
00212 
00213                         for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
00214                                 for(int j=0;j<3;++j)
00215                                         if(!(*fi).IsB(j)) {
00216                                                 CoordType edge= (*fi).V1(j)->P() -(*fi).V(j)->P();
00217                                                 ScalarType len=Norm(edge);
00218                                                 edge/=len;
00219                                                 TD[(*fi).V(j)].PntSum+=edge;
00220                                                 TD[(*fi).V1(j)].PntSum-=edge;
00221                                                 TD[(*fi).V(j)].LenSum+=len;
00222                                                 TD[(*fi).V1(j)].LenSum+=len;
00223                                         }
00224 
00225                         for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
00226                                 for(int j=0;j<3;++j)
00227                                         // se l'edge j e' di bordo si riazzera tutto e si riparte
00228                                         if((*fi).IsB(j)) {
00229                                                         TD[(*fi).V(j)].PntSum=CoordType(0,0,0);
00230                                                         TD[(*fi).V1(j)].PntSum=CoordType(0,0,0);
00231                                                         TD[(*fi).V(j)].LenSum=0;
00232                                                         TD[(*fi).V1(j)].LenSum=0;
00233                                         }
00234 
00235 
00236                         for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
00237                                         for(int j=0;j<3;++j)
00238                                                 if((*fi).IsB(j))
00239                                                 {
00240                                                         CoordType edge= (*fi).V1(j)->P() -(*fi).V(j)->P();
00241                                                         ScalarType len=Norm(edge);
00242                                                         edge/=len;
00243                                                         TD[(*fi).V(j)].PntSum+=edge;
00244                                                         TD[(*fi).V1(j)].PntSum-=edge;
00245                                                         TD[(*fi).V(j)].LenSum+=len;
00246                                                         TD[(*fi).V1(j)].LenSum+=len;
00247                                                 }
00248   // The fundamental part:
00249         // We move the new point of a quantity
00250         //
00251         //  L(M) = 1/Sum(edgelen) * Sum(Normalized edges)
00252         //
00253 
00254                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00255                                 if(!(*vi).IsD() && TD[*vi].LenSum>0 )
00256                                  (*vi).P() = (*vi).P() + (TD[*vi].PntSum/TD[*vi].LenSum)*delta;
00257         }
00258 };
00259 
00260 
00261 class LaplacianInfo
00262 {
00263 public:
00264         LaplacianInfo(const CoordType &_p, const int _n):sum(_p),cnt(_n) {}
00265         LaplacianInfo() {}
00266         CoordType sum;
00267         ScalarType cnt;
00268 };
00269 
00270 // Classical Laplacian Smoothing. Each vertex can be moved onto the average of the adjacent vertices.
00271 // Can smooth only the selected vertices and weight the smoothing according to the quality
00272 // In the latter case 0 means that the vertex is not moved and 1 means that the vertex is moved onto the computed position.
00273 //
00274 // From the Taubin definition "A signal proc approach to fair surface design"
00275 // We define the discrete Laplacian of a discrete surface signal by weighted averages over the neighborhoods
00276 //          \delta xi  = \Sum wij (xj - xi) ;
00277 // where xj are the adjacent vertices of xi and wij is usually 1/n_adj
00278 //
00279 // This function simply accumulate over a TempData all the positions of the ajacent vertices
00280 //
00281 static void AccumulateLaplacianInfo(MeshType &m, SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > &TD)
00282 {
00283                         FaceIterator fi;
00284                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00285                         {
00286                                 if(!(*fi).IsD())
00287                                         for(int j=0;j<3;++j)
00288                                                 if(!(*fi).IsB(j))
00289                                                 {
00290                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->P();
00291                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->P();
00292                                                         ++TD[(*fi).V(j)].cnt;
00293                                                         ++TD[(*fi).V1(j)].cnt;
00294                                                 }
00295                         }
00296                         // si azzaera i dati per i vertici di bordo
00297                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00298                         {
00299                                 if(!(*fi).IsD())
00300                                         for(int j=0;j<3;++j)
00301                                                 if((*fi).IsB(j))
00302                                                 {
00303                                                         TD[(*fi).V0(j)].sum=(*fi).P0(j);
00304                                                         TD[(*fi).V1(j)].sum=(*fi).P1(j);
00305                                                         TD[(*fi).V0(j)].cnt=1;
00306                                                         TD[(*fi).V1(j)].cnt=1;
00307                                                 }
00308                         }
00309 
00310                         // se l'edge j e' di bordo si deve mediare solo con gli adiacenti
00311                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00312                         {
00313                                 if(!(*fi).IsD())
00314                                         for(int j=0;j<3;++j)
00315                                                 if((*fi).IsB(j))
00316                                                 {
00317                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->P();
00318                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->P();
00319                                                         ++TD[(*fi).V(j)].cnt;
00320                                                         ++TD[(*fi).V1(j)].cnt;
00321                                                 }
00322                         }
00323 }
00324 
00325 static void VertexCoordLaplacian(MeshType &m, int step, bool SmoothSelected=false, vcg::CallBackPos * cb=0)
00326 {
00327   VertexIterator vi;
00328         LaplacianInfo lpz(CoordType(0,0,0),0);
00329         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00330         for(int i=0;i<step;++i)
00331                 {
00332                         if(cb)cb(100*i/step, "Classic Laplacian Smoothing");
00333                         TD.Init(lpz);
00334                         AccumulateLaplacianInfo(m,TD);
00335                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00336                                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00337                                 {
00338                                                         if(!SmoothSelected || (*vi).IsS())
00339                                                                         (*vi).P() = ( (*vi).P() + TD[*vi].sum)/(TD[*vi].cnt+1);
00340                                 }
00341                 }
00342 }
00343 
00344 // Same of above but moves only the vertices that do not change FaceOrientation more that the given threshold
00345 static void VertexCoordPlanarLaplacian(MeshType &m, int step, float AngleThrRad = math::ToRad(1.0), bool SmoothSelected=false, vcg::CallBackPos * cb=0)
00346 {
00347   VertexIterator vi;
00348         FaceIterator fi;
00349         LaplacianInfo lpz(CoordType(0,0,0),0);
00350         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00351         for(int i=0;i<step;++i)
00352                 {
00353                 if(cb)cb(100*i/step, "Planar Laplacian Smoothing");
00354                 TD.Init(lpz);
00355                 AccumulateLaplacianInfo(m,TD);
00356                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00357                         if(!(*vi).IsD() && TD[*vi].cnt>0 )
00358                         {
00359                                 if(!SmoothSelected || (*vi).IsS())
00360                                         TD[*vi].sum = ( (*vi).P() + TD[*vi].sum)/(TD[*vi].cnt+1);
00361                         }
00362 
00363                 for(fi=m.face.begin();fi!=m.face.end();++fi){
00364                                 if(!(*fi).IsD()){
00365                                         for (int j = 0; j < 3; ++j) {
00366                                                 if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, (*fi).P1(j), (*fi).P2(j) ),
00367                                                                                         NormalizedNormal(   (*fi).P0(j)     , (*fi).P1(j), (*fi).P2(j) ) ) > AngleThrRad )
00368                                                         TD[(*fi).V0(j)].sum = (*fi).P0(j);
00369                                         }
00370                                 }
00371                         }
00372                         for(fi=m.face.begin();fi!=m.face.end();++fi){
00373                                 if(!(*fi).IsD()){
00374                                         for (int j = 0; j < 3; ++j) {
00375                                                 if(Angle( NormalizedNormal(TD[(*fi).V0(j)].sum, TD[(*fi).V1(j)].sum, (*fi).P2(j) ),
00376                                                                                         NormalizedNormal(   (*fi).P0(j)     ,    (*fi).P1(j),      (*fi).P2(j) ) ) > AngleThrRad )
00377                                                 {
00378                                                         TD[(*fi).V0(j)].sum = (*fi).P0(j);
00379                                                         TD[(*fi).V1(j)].sum = (*fi).P1(j);
00380                                                 }
00381                                         }
00382                                 }
00383                         }
00384 
00385                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00386                                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00387                                                 (*vi).P()= TD[*vi].sum;
00388 
00389 
00390 
00391                 }// end step
00392 
00393 
00394 }
00395 
00396 static void VertexCoordLaplacianBlend(MeshType &m, int step, float alpha, bool SmoothSelected=false)
00397 {
00398         VertexIterator vi;
00399         LaplacianInfo lpz(CoordType(0,0,0),0);
00400         assert (alpha<= 1.0);
00401         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert);
00402 
00403         for(int i=0;i<step;++i)
00404                 {
00405                 TD.Init(lpz);
00406                 AccumulateLaplacianInfo(m,TD);
00407                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00408                         if(!(*vi).IsD() && TD[*vi].cnt>0 )
00409                         {
00410                                 if(!SmoothSelected || (*vi).IsS())
00411                                 {
00412                                         CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P();
00413                                         (*vi).P() = (*vi).P() + Delta*alpha;
00414                                 }
00415                         }
00416                 }
00417 }
00418 
00419 /* a couple of notes about the lambda mu values
00420 We assume that 0 < lambda , and mu  is a negative scale factor such that mu < - lambda.
00421 Holds mu+lambda < 0 (e.g in absolute value mu is greater)
00422 
00423 let kpb be the pass-band frequency, taubin says that:
00424                         kpb = 1/lambda + 1/mu >0
00425 
00426 Values of kpb from 0.01 to 0.1 produce good results according to the original paper.
00427 
00428 kpb * mu - mu/lambda = 1
00429 mu = 1/(kpb-1/lambda )
00430 
00431 So if
00432 * lambda == 0.5, kpb==0.1  -> mu = 1/(0.1 - 2)  = -0.526
00433 * lambda == 0.5, kpb==0.01 -> mu = 1/(0.01 - 2) = -0.502
00434 */
00435 
00436 
00437 static void VertexCoordTaubin(MeshType &m, int step, float lambda, float mu, bool SmoothSelected=false, vcg::CallBackPos * cb=0)
00438 {
00439         LaplacianInfo lpz(CoordType(0,0,0),0);
00440         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00441         VertexIterator vi;
00442         for(int i=0;i<step;++i)
00443                 {
00444                   if(cb) cb(100*i/step, "Taubin Smoothing");
00445                         TD.Init(lpz);
00446                         AccumulateLaplacianInfo(m,TD);
00447                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00448                                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00449                                 {
00450                                         if(!SmoothSelected || (*vi).IsS())
00451                                         {
00452                                                 CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P();
00453                                                 (*vi).P() = (*vi).P() + Delta*lambda ;
00454                                         }
00455                                 }
00456                         TD.Init(lpz);
00457                         AccumulateLaplacianInfo(m,TD);
00458                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00459                                         if(!(*vi).IsD() && TD[*vi].cnt>0 )
00460                                         {
00461                                                 if(!SmoothSelected || (*vi).IsS())
00462                                                 {
00463                                                         CoordType Delta = TD[*vi].sum/TD[*vi].cnt - (*vi).P();
00464                                                         (*vi).P() = (*vi).P() + Delta*mu ;
00465                                                 }
00466                                         }
00467                         } // end for step
00468 }
00469 
00470 
00471 static void VertexCoordLaplacianQuality(MeshType &m, int step, bool SmoothSelected=false)
00472 {
00473   LaplacianInfo lpz;
00474         lpz.sum=CoordType(0,0,0);
00475         lpz.cnt=1;
00476         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00477         for(int i=0;i<step;++i)
00478                 {
00479                         for(VertexIterator vi=m.vert.begin();vi!=m.vert.end();++vi)
00480                                         if(!(*vi).IsD() && TD[*vi].cnt>0 )
00481                                                 if(!SmoothSelected || (*vi).IsS())
00482                                                 {
00483                                                         float q=(*vi).Q();
00484                                                         (*vi).P()=(*vi).P()*q + (TD[*vi].sum/TD[*vi].cnt)*(1.0-q);
00485                                                 }
00486                 } // end for
00487 };
00488 
00489 /*
00490   Improved Laplacian Smoothing of Noisy Surface Meshes
00491   J. Vollmer, R. Mencl, and H. M�ller
00492   EUROGRAPHICS Volume 18 (1999), Number 3
00493 */
00494 
00495 class HCSmoothInfo
00496 {
00497 public:
00498         CoordType dif;
00499         CoordType sum;
00500         int cnt;
00501 };
00502 
00503 static void VertexCoordLaplacianHC(MeshType &m, int step, bool SmoothSelected=false )
00504 {
00505         ScalarType beta=0.5;
00506   HCSmoothInfo lpz;
00507         lpz.sum=CoordType(0,0,0);
00508         lpz.dif=CoordType(0,0,0);
00509         lpz.cnt=0;
00510                 for(int i=0;i<step;++i)
00511                 {
00512                         SimpleTempData<typename MeshType::VertContainer,HCSmoothInfo > TD(m.vert,lpz);
00513                         // First Loop compute the laplacian
00514                         FaceIterator fi;
00515                         for(fi=m.face.begin();fi!=m.face.end();++fi)if(!(*fi).IsD())
00516                                 {
00517                                         for(int j=0;j<3;++j)
00518                                         {
00519                                                 TD[(*fi).V(j)].sum+=(*fi).V1(j)->P();
00520                                                 TD[(*fi).V1(j)].sum+=(*fi).V(j)->P();
00521                                                 ++TD[(*fi).V(j)].cnt;
00522                                                 ++TD[(*fi).V1(j)].cnt;
00523                                                 // se l'edge j e' di bordo si deve sommare due volte
00524                                                 if((*fi).IsB(j))
00525                                                 {
00526                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->P();
00527                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->P();
00528                                                         ++TD[(*fi).V(j)].cnt;
00529                                                         ++TD[(*fi).V1(j)].cnt;
00530                                                 }
00531                                         }
00532                                 }
00533                         VertexIterator vi;
00534                         for(vi=m.vert.begin();vi!=m.vert.end();++vi) if(!(*vi).IsD())
00535                                  TD[*vi].sum/=(float)TD[*vi].cnt;
00536 
00537                         // Second Loop compute average difference
00538                         for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
00539                                 {
00540                                         for(int j=0;j<3;++j)
00541                                         {
00542                                                 TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P();
00543                                                 TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P();
00544                                                 // se l'edge j e' di bordo si deve sommare due volte
00545                                                 if((*fi).IsB(j))
00546                                                 {
00547                                                         TD[(*fi).V(j)].dif +=TD[(*fi).V1(j)].sum-(*fi).V1(j)->P();
00548                                                         TD[(*fi).V1(j)].dif+=TD[(*fi).V(j)].sum-(*fi).V(j)->P();
00549                                                 }
00550                                         }
00551                                 }
00552 
00553                         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00554                                 {
00555                                  TD[*vi].dif/=(float)TD[*vi].cnt;
00556                                  if(!SmoothSelected || (*vi).IsS())
00557                                                 (*vi).P()= TD[*vi].sum - (TD[*vi].sum - (*vi).P())*beta  + (TD[*vi].dif)*(1.f-beta);
00558                                 }
00559                 } // end for step
00560 };
00561 
00562 // Laplacian smooth of the quality.
00563 
00564 
00565 class ColorSmoothInfo
00566 {
00567 public:
00568         unsigned int r;
00569         unsigned int g;
00570         unsigned int b;
00571         unsigned int a;
00572         int cnt;
00573 };
00574 
00575 static void VertexColorLaplacian(MeshType &m, int step, bool SmoothSelected=false, vcg::CallBackPos * cb=0)
00576 {
00577         ColorSmoothInfo csi;
00578         csi.r=0; csi.g=0; csi.b=0; csi.cnt=0;
00579         SimpleTempData<typename MeshType::VertContainer, ColorSmoothInfo> TD(m.vert,csi);
00580 
00581         for(int i=0;i<step;++i)
00582         {
00583                 if(cb) cb(100*i/step, "Vertex Color Laplacian Smoothing");
00584                 VertexIterator vi;
00585                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00586                         TD[*vi]=csi;
00587 
00588                 FaceIterator fi;
00589                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00590                         if(!(*fi).IsD())
00591                                 for(int j=0;j<3;++j)
00592                                         if(!(*fi).IsB(j))
00593                                         {
00594                                                 TD[(*fi).V(j)].r+=(*fi).V1(j)->C()[0];
00595                                                 TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1];
00596                                                 TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2];
00597                                                 TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3];
00598 
00599                                                 TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0];
00600                                                 TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1];
00601                                                 TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2];
00602                                                 TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3];
00603 
00604                                                 ++TD[(*fi).V(j)].cnt;
00605                                                 ++TD[(*fi).V1(j)].cnt;
00606                                         }
00607 
00608                 // si azzaera i dati per i vertici di bordo
00609                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00610                         if(!(*fi).IsD())
00611                                 for(int j=0;j<3;++j)
00612                                         if((*fi).IsB(j))
00613                                         {
00614                                                 TD[(*fi).V(j)]=csi;
00615                                                 TD[(*fi).V1(j)]=csi;
00616                                         }
00617 
00618                 // se l'edge j e' di bordo si deve mediare solo con gli adiacenti
00619                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00620                         if(!(*fi).IsD())
00621                                 for(int j=0;j<3;++j)
00622                                         if((*fi).IsB(j))
00623                                         {
00624                                                 TD[(*fi).V(j)].r+=(*fi).V1(j)->C()[0];
00625                                                 TD[(*fi).V(j)].g+=(*fi).V1(j)->C()[1];
00626                                                 TD[(*fi).V(j)].b+=(*fi).V1(j)->C()[2];
00627                                                 TD[(*fi).V(j)].a+=(*fi).V1(j)->C()[3];
00628 
00629                                                 TD[(*fi).V1(j)].r+=(*fi).V(j)->C()[0];
00630                                                 TD[(*fi).V1(j)].g+=(*fi).V(j)->C()[1];
00631                                                 TD[(*fi).V1(j)].b+=(*fi).V(j)->C()[2];
00632                                                 TD[(*fi).V1(j)].a+=(*fi).V(j)->C()[3];
00633 
00634                                                 ++TD[(*fi).V(j)].cnt;
00635                                                 ++TD[(*fi).V1(j)].cnt;
00636                                         }
00637 
00638                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00639                         if(!(*vi).IsD() && TD[*vi].cnt>0 )
00640                                 if(!SmoothSelected || (*vi).IsS())
00641                                 {
00642                                         (*vi).C()[0] = (unsigned int) ceil((double) (TD[*vi].r / TD[*vi].cnt));
00643                                         (*vi).C()[1] = (unsigned int) ceil((double) (TD[*vi].g / TD[*vi].cnt));
00644                                         (*vi).C()[2] = (unsigned int) ceil((double) (TD[*vi].b / TD[*vi].cnt));
00645                                         (*vi).C()[3] = (unsigned int) ceil((double) (TD[*vi].a / TD[*vi].cnt));
00646                                 }
00647         } // end for step
00648 };
00649 
00650 static void FaceColorLaplacian(MeshType &m, int step, bool SmoothSelected=false, vcg::CallBackPos * cb=0)
00651 {
00652         ColorSmoothInfo csi;
00653         csi.r=0; csi.g=0; csi.b=0; csi.cnt=0;
00654         SimpleTempData<typename MeshType::FaceContainer, ColorSmoothInfo> TD(m.face,csi);
00655 
00656         for(int i=0;i<step;++i)
00657         {
00658                 if(cb) cb(100*i/step, "Face Color Laplacian Smoothing");
00659                 FaceIterator fi;
00660                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00661                         TD[*fi]=csi;
00662 
00663                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00664                 {
00665                         if(!(*fi).IsD())
00666                                 for(int j=0;j<3;++j)
00667                                         if(!(*fi).IsB(j))
00668                                         {
00669                                                 TD[*fi].r+=(*fi).FFp(j)->C()[0];
00670                                                 TD[*fi].g+=(*fi).FFp(j)->C()[1];
00671                                                 TD[*fi].b+=(*fi).FFp(j)->C()[2];
00672                                                 TD[*fi].a+=(*fi).FFp(j)->C()[3];
00673                                                 ++TD[*fi].cnt;
00674                                         }
00675                 }
00676                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00677                         if(!(*fi).IsD() && TD[*fi].cnt>0 )
00678                                 if(!SmoothSelected || (*fi).IsS())
00679                                 {
00680                                         (*fi).C()[0] = (unsigned int) ceil((float) (TD[*fi].r / TD[*fi].cnt));
00681                                         (*fi).C()[1] = (unsigned int) ceil((float) (TD[*fi].g / TD[*fi].cnt));
00682                                         (*fi).C()[2] = (unsigned int) ceil((float) (TD[*fi].b / TD[*fi].cnt));
00683                                         (*fi).C()[3] = (unsigned int) ceil((float) (TD[*fi].a / TD[*fi].cnt));
00684                                 }
00685         } // end for step
00686 };
00687 
00688 // Laplacian smooth of the quality.
00689 
00690 class QualitySmoothInfo
00691 {
00692 public:
00693         ScalarType sum;
00694         int cnt;
00695 };
00696 
00697 
00698 static void VertexQualityLaplacian(MeshType &m, int step=1, bool SmoothSelected=false)
00699 {
00700   QualitySmoothInfo lpz;
00701         lpz.sum=0;
00702         lpz.cnt=0;
00703         SimpleTempData<typename MeshType::VertContainer,QualitySmoothInfo> TD(m.vert,lpz);
00704         //TD.Start(lpz);
00705         for(int i=0;i<step;++i)
00706         {
00707                 VertexIterator vi;
00708                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00709                          TD[*vi]=lpz;
00710 
00711                 FaceIterator fi;
00712                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00713                         if(!(*fi).IsD())
00714                                 for(int j=0;j<3;++j)
00715                                         if(!(*fi).IsB(j))
00716                                                 {
00717                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q();
00718                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q();
00719                                                         ++TD[(*fi).V(j)].cnt;
00720                                                         ++TD[(*fi).V1(j)].cnt;
00721                                         }
00722 
00723                         // si azzaera i dati per i vertici di bordo
00724                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00725                                 if(!(*fi).IsD())
00726                                         for(int j=0;j<3;++j)
00727                                                 if((*fi).IsB(j))
00728                                                         {
00729                                                                 TD[(*fi).V(j)]=lpz;
00730                                                                 TD[(*fi).V1(j)]=lpz;
00731                                                         }
00732 
00733                         // se l'edge j e' di bordo si deve mediare solo con gli adiacenti
00734                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00735                                 if(!(*fi).IsD())
00736                                         for(int j=0;j<3;++j)
00737                                                 if((*fi).IsB(j))
00738                                                         {
00739                                                                 TD[(*fi).V(j)].sum+=(*fi).V1(j)->Q();
00740                                                                 TD[(*fi).V1(j)].sum+=(*fi).V(j)->Q();
00741                                                                 ++TD[(*fi).V(j)].cnt;
00742                                                                 ++TD[(*fi).V1(j)].cnt;
00743                                                 }
00744 
00745         //VertexIterator vi;
00746         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00747                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00748                         if(!SmoothSelected || (*vi).IsS())
00749                                         (*vi).Q()=TD[*vi].sum/TD[*vi].cnt;
00750         }
00751 
00752         //TD.Stop();
00753 };
00754 
00755 static void VertexNormalLaplacian(MeshType &m, int step,bool SmoothSelected=false)
00756 {
00757     LaplacianInfo lpz;
00758       lpz.sum=CoordType(0,0,0);
00759       lpz.cnt=0;
00760     SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00761     for(int i=0;i<step;++i)
00762         {
00763                 VertexIterator vi;
00764                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00765                          TD[*vi]=lpz;
00766 
00767                 FaceIterator fi;
00768                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00769                         if(!(*fi).IsD())
00770                                 for(int j=0;j<3;++j)
00771                                         if(!(*fi).IsB(j))
00772                                                 {
00773                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->N();
00774                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->N();
00775                                                         ++TD[(*fi).V(j)].cnt;
00776                                                         ++TD[(*fi).V1(j)].cnt;
00777                                         }
00778 
00779                         // si azzaera i dati per i vertici di bordo
00780                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00781                                 if(!(*fi).IsD())
00782                                         for(int j=0;j<3;++j)
00783                                                 if((*fi).IsB(j))
00784                                                         {
00785                                                                 TD[(*fi).V(j)]=lpz;
00786                                                                 TD[(*fi).V1(j)]=lpz;
00787                                                         }
00788 
00789                         // se l'edge j e' di bordo si deve mediare solo con gli adiacenti
00790                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00791                                 if(!(*fi).IsD())
00792                                         for(int j=0;j<3;++j)
00793                                                 if((*fi).IsB(j))
00794                                                         {
00795                                                                 TD[(*fi).V(j)].sum+=(*fi).V1(j)->N();
00796                                                                 TD[(*fi).V1(j)].sum+=(*fi).V(j)->N();
00797                                                                 ++TD[(*fi).V(j)].cnt;
00798                                                                 ++TD[(*fi).V1(j)].cnt;
00799                                                 }
00800 
00801         //VertexIterator vi;
00802         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00803                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00804                         if(!SmoothSelected || (*vi).IsS())
00805                                         (*vi).N()=TD[*vi].sum/TD[*vi].cnt;
00806         }
00807 
00808 };
00809 
00810 // Smooth solo lungo la direzione di vista
00811         // alpha e' compreso fra 0(no smoot) e 1 (tutto smoot)
00812   // Nota che se smootare il bordo puo far fare bandierine.
00813 static void VertexCoordViewDepth(MeshType &m,
00814                                  const CoordType & viewpoint,
00815                                  const ScalarType alpha,
00816                                  int step, bool SmoothBorder=false )
00817 {
00818         LaplacianInfo lpz;
00819         lpz.sum=CoordType(0,0,0);
00820         lpz.cnt=0;
00821         SimpleTempData<typename MeshType::VertContainer,LaplacianInfo > TD(m.vert,lpz);
00822         for(int i=0;i<step;++i)
00823         {
00824                 VertexIterator vi;
00825                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00826                          TD[*vi]=lpz;
00827 
00828                 FaceIterator fi;
00829                 for(fi=m.face.begin();fi!=m.face.end();++fi)
00830                         if(!(*fi).IsD())
00831                                 for(int j=0;j<3;++j)
00832                                         if(!(*fi).IsB(j))
00833                                                 {
00834                                                         TD[(*fi).V(j)].sum+=(*fi).V1(j)->cP();
00835                                                         TD[(*fi).V1(j)].sum+=(*fi).V(j)->cP();
00836                                                         ++TD[(*fi).V(j)].cnt;
00837                                                         ++TD[(*fi).V1(j)].cnt;
00838                                         }
00839 
00840                         // si azzaera i dati per i vertici di bordo
00841                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00842                                 if(!(*fi).IsD())
00843                                         for(int j=0;j<3;++j)
00844                                                 if((*fi).IsB(j))
00845                                                         {
00846                                                                 TD[(*fi).V(j)]=lpz;
00847                                                                 TD[(*fi).V1(j)]=lpz;
00848                                                         }
00849 
00850                         // se l'edge j e' di bordo si deve mediare solo con gli adiacenti
00851      if(SmoothBorder)
00852                         for(fi=m.face.begin();fi!=m.face.end();++fi)
00853                                 if(!(*fi).IsD())
00854                                         for(int j=0;j<3;++j)
00855                                                 if((*fi).IsB(j))
00856                                                         {
00857                                                                 TD[(*fi).V(j)].sum+=(*fi).V1(j)->cP();
00858                                                                 TD[(*fi).V1(j)].sum+=(*fi).V(j)->cP();
00859                                                                 ++TD[(*fi).V(j)].cnt;
00860                                                                 ++TD[(*fi).V1(j)].cnt;
00861                                                 }
00862 
00863         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
00864                 if(!(*vi).IsD() && TD[*vi].cnt>0 )
00865                         {
00866                                 CoordType np = TD[*vi].sum/TD[*vi].cnt;
00867                                 CoordType d = (*vi).cP() - viewpoint; d.Normalize();
00868                                 ScalarType s = d .dot ( np - (*vi).cP() );
00869                                 (*vi).P() += d * (s*alpha);
00870                         }
00871         }
00872 }
00873 
00874 
00875 
00876 /****************************************************************************************************************/
00877 /****************************************************************************************************************/
00878 // Paso Double Smoothing
00879 //  The proposed
00880 // approach is a two step method where  in the first step the face normals
00881 // are adjusted and then, in a second phase, the vertex positions are updated.
00882 /****************************************************************************************************************/
00883 /****************************************************************************************************************/
00884 // Classi di info
00885 
00886 class PDVertInfo
00887 {
00888 public:
00889         CoordType np;
00890 };
00891 
00892 
00893 class PDFaceInfo
00894 {
00895 public:
00896         CoordType m;
00897 };
00898 /***************************************************************************/
00899 // Paso Doble Step 1 compute the smoothed normals
00900 /***************************************************************************/
00901 // Requirements:
00902 // VF Topology
00903 // Normalized Face Normals
00904 //
00905 // This is the Normal Smoothing approach of Shen and Berner
00906 // Fuzzy Vector Median-Based Surface Smoothing TVCG 2004
00907 
00908 
00909 void FaceNormalFuzzyVectorSB(MeshType &m,
00910                                   SimpleTempData<typename MeshType::FaceContainer,PDFaceInfo > &TD,
00911                                   ScalarType sigma)
00912 {
00913         int i;
00914 
00915         FaceIterator fi;
00916 
00917         for(fi=m.face.begin();fi!=m.face.end();++fi)
00918         {
00919                 CoordType bc=(*fi).Barycenter();
00920     // 1) Clear all the visited flag of faces that are vertex-adjacent to fi
00921                 for(i=0;i<3;++i)
00922                 {
00923         vcg::face::VFIterator<FaceType> ep(&*fi,i);
00924                   while (!ep.End())
00925                         {
00926                                 ep.f->ClearV();
00927                                 ++ep;
00928                         }
00929                 }
00930 
00931     // 1) Effectively average the normals weighting them with
00932     (*fi).SetV();
00933                 CoordType mm=CoordType(0,0,0);
00934                 for(i=0;i<3;++i)
00935                 {
00936                   vcg::face::VFIterator<FaceType> ep(&*fi,i);
00937                   while (!ep.End())
00938                         {
00939                                 if(! (*ep.f).IsV() )
00940                                 {
00941                                         if(sigma>0)
00942                                         {
00943                                                 ScalarType dd=SquaredDistance(ep.f->Barycenter(),bc);
00944                                                 ScalarType ang=AngleN(ep.f->N(),(*fi).N());
00945                                                 mm+=ep.f->N()*exp((-sigma)*ang*ang/dd);
00946                                         }
00947                                         else mm+=ep.f->N();
00948                                         (*ep.f).SetV();
00949                                 }
00950                                 ++ep;
00951                         }
00952                 }
00953                 mm.Normalize();
00954                 TD[*fi].m=mm;
00955         }
00956 }
00957 
00958 // Replace the normal of the face with the average of normals of the vertex adijacent faces.
00959 // Normals are weighted with face area.
00960 // It assumes that:
00961 // Normals are normalized:
00962 // VF adjacency is present.
00963 
00964 static void FaceNormalLaplacianVF(MeshType &m)
00965 {
00966         SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face);
00967 
00968         PDFaceInfo lpzf;
00969         lpzf.m=CoordType(0,0,0);
00970 
00971         assert(tri::HasPerVertexVFAdjacency(m) && tri::HasPerFaceVFAdjacency(m) );
00972         TDF.Start(lpzf);
00973         int i;
00974 
00975   FaceIterator fi;
00976 
00977         tri::UpdateNormals<MeshType>::AreaNormalizeFace(m);
00978 
00979         for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
00980         {
00981                                 CoordType bc=Barycenter<FaceType>(*fi);
00982                                 // 1) Clear all the visited flag of faces that are vertex-adjacent to fi
00983                                 for(i=0;i<3;++i)
00984                                 {
00985                                         VFLocalIterator ep(&*fi,i);
00986                                         for (;!ep.End();++ep)
00987                                                 ep.f->ClearV();
00988                                 }
00989 
00990                                 // 2) Effectively average the normals
00991                                 CoordType normalSum=(*fi).N();
00992 
00993                                 for(i=0;i<3;++i)
00994                                 {
00995                                         VFLocalIterator ep(&*fi,i);
00996                                         for (;!ep.End();++ep)
00997                                         {
00998                                                 if(! (*ep.f).IsV() )
00999                                                                 {
01000                                                                                  normalSum += ep.f->N();
01001                                                                                  (*ep.f).SetV();
01002                                                                 }
01003                                         }
01004                                 }
01005                                 normalSum.Normalize();
01006                                 TDF[*fi].m=normalSum;
01007         }
01008         for(fi=m.face.begin();fi!=m.face.end();++fi)
01009                 (*fi).N()=TDF[*fi].m;
01010 
01011         tri::UpdateNormals<MeshType>::NormalizeFace(m);
01012 
01013         TDF.Stop();
01014 }
01015 
01016 // Replace the normal of the face with the average of normals of the face adijacent faces.
01017 // Normals are weighted with face area.
01018 // It assumes that:
01019 // Normals are normalized:
01020 // FF adjacency is present.
01021 
01022 
01023 static void FaceNormalLaplacianFF(MeshType &m, int step=1, bool SmoothSelected=false )
01024 {
01025         PDFaceInfo lpzf;
01026         lpzf.m=CoordType(0,0,0);
01027         SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face,lpzf);
01028         assert(tri::HasFFAdjacency(m));
01029 
01030   FaceIterator fi;
01031         tri::UpdateNormals<MeshType>::AreaNormalizeFace(m);
01032   for(int i=0;i<step;++i)
01033                 {
01034                                 for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
01035                                 {
01036                                                         CoordType normalSum=(*fi).N();
01037 
01038                                                         for(i=0;i<3;++i)
01039                                                                                 normalSum+=(*fi).FFp(i)->N();
01040 
01041                                                         TDF[*fi].m=normalSum;
01042                                 }
01043                                 for(fi=m.face.begin();fi!=m.face.end();++fi)
01044                                                 if(!SmoothSelected || (*fi).IsS())
01045                                                                 (*fi).N()=TDF[*fi].m;
01046 
01047                                 tri::UpdateNormals<MeshType>::NormalizeFace(m);
01048                 }
01049 }
01050 
01051 
01052 /***************************************************************************/
01053 // Paso Doble Step 1 compute the smoothed normals
01054 /***************************************************************************/
01055 // Requirements:
01056 // VF Topology
01057 // Normalized Face Normals
01058 //
01059 // This is the Normal Smoothing approach bsased on a angle thresholded weighting
01060 // sigma is in the 0 .. 1 range, it represent the cosine of a threshold angle.
01061 // sigma == 0 All the normals are averaged
01062 // sigma == 1 Nothing is averaged.
01063 // Only within the specified range are averaged toghether. The averagin is weighted with the
01064 
01065 
01066 static void FaceNormalAngleThreshold(MeshType &m,
01067                                   SimpleTempData<typename MeshType::FaceContainer,PDFaceInfo> &TD,
01068                                   ScalarType sigma)
01069 {
01070         int i;
01071 
01072 
01073   FaceIterator fi;
01074 
01075         for(fi=m.face.begin();fi!=m.face.end();++fi) if(!(*fi).IsD())
01076         {
01077                 CoordType bc=Barycenter<FaceType>(*fi);
01078     // 1) Clear all the visited flag of faces that are vertex-adjacent to fi
01079                 for(i=0;i<3;++i)
01080                 {
01081         VFLocalIterator ep(&*fi,i);
01082                   for (;!ep.End();++ep)
01083                                 ep.f->ClearV();
01084                 }
01085 
01086     // 1) Effectively average the normals weighting them with the squared difference of the angle similarity
01087                 // sigma is the cosine of a threshold angle. sigma \in 0..1
01088                 // sigma == 0 All the normals are averaged
01089                 // sigma == 1 Nothing is averaged.
01090                 // The averaging is weighted with the difference between normals. more similar the normal more important they are.
01091 
01092     CoordType normalSum=CoordType(0,0,0);
01093                 for(i=0;i<3;++i)
01094                 {
01095         VFLocalIterator ep(&*fi,i);
01096                   for (;!ep.End();++ep)
01097                         {
01098                                 if(! (*ep.f).IsV() )
01099                                 {
01100           ScalarType cosang=ep.f->N().dot((*fi).N());
01101                                         // Note that if two faces form an angle larger than 90 deg, their contribution should be very very small.
01102                                         // Without this clamping
01103                                         cosang = math::Clamp(cosang,0.0001f,1.f);
01104           if(cosang >= sigma)
01105           {
01106             ScalarType w = cosang-sigma;
01107                                                 normalSum += ep.f->N()*(w*w);   // similar normals have a cosang very close to 1 so cosang - sigma is maximized
01108           }
01109                                         (*ep.f).SetV();
01110                                 }
01111                         }
01112                 }
01113                 normalSum.Normalize();
01114                 TD[*fi].m=normalSum;
01115         }
01116 
01117   for(fi=m.face.begin();fi!=m.face.end();++fi)
01118     (*fi).N()=TD[*fi].m;
01119 }
01120 
01121 /****************************************************************************************************************/
01122 // Restituisce il gradiente dell'area del triangolo nel punto p.
01123 // Nota che dovrebbe essere sempre un vettore che giace nel piano del triangolo e perpendicolare al lato opposto al vertice p.
01124 // Ottimizzato con Maple e poi pesantemente a mano.
01125 
01126 CoordType TriAreaGradient(CoordType &p,CoordType &p0,CoordType &p1)
01127 {
01128         CoordType dd = p1-p0;
01129         CoordType d0 = p-p0;
01130         CoordType d1 = p-p1;
01131         CoordType grad;
01132 
01133         ScalarType t16 =  d0[1]* d1[2] - d0[2]* d1[1];
01134         ScalarType t5  = -d0[2]* d1[0] + d0[0]* d1[2];
01135         ScalarType t4  = -d0[0]* d1[1] + d0[1]* d1[0];
01136 
01137         ScalarType delta= sqrtf(t4*t4 + t5*t5 +t16*t16);
01138 
01139         grad[0]= (t5  * (-dd[2]) + t4 * ( dd[1]))/delta;
01140         grad[1]= (t16 * (-dd[2]) + t4 * (-dd[0]))/delta;
01141         grad[2]= (t16 * ( dd[1]) + t5 * ( dd[0]))/delta;
01142 
01143         return grad;
01144 }
01145 
01146 template <class ScalarType>
01147 CoordType CrossProdGradient(CoordType &p, CoordType &p0, CoordType &p1, CoordType &m)
01148 {
01149         CoordType grad;
01150         CoordType p00=p0-p;
01151         CoordType p01=p1-p;
01152         grad[0] = (-p00[2] + p01[2])*m[1] + (-p01[1] + p00[1])*m[2];
01153         grad[1] = (-p01[2] + p00[2])*m[0] + (-p00[0] + p01[0])*m[2];
01154         grad[2] = (-p00[1] + p01[1])*m[0] + (-p01[0] + p00[0])*m[1];
01155 
01156         return grad;
01157 }
01158 
01159 /*
01160 Deve Calcolare il gradiente di
01161 E(p) = A(p,p0,p1) (n - m)^2 =
01162 A(...) (2-2nm)   =
01163 (p0-p)^(p1-p)
01164 2A - 2A * ------------- m  =
01165 2A
01166 
01167 2A  -  2 (p0-p)^(p1-p) * m
01168 */
01169 
01170 CoordType FaceErrorGrad(CoordType &p,CoordType &p0,CoordType &p1, CoordType &m)
01171 {
01172         return     TriAreaGradient(p,p0,p1) *2.0f
01173                 - CrossProdGradient(p,p0,p1,m) *2.0f ;
01174 }
01175 /***************************************************************************/
01176 // Paso Doble Step 2 Fitta la mesh a un dato insieme di normali
01177 /***************************************************************************/
01178 
01179 
01180 void FitMesh(MeshType &m,
01181                          SimpleTempData<typename MeshType::VertContainer, PDVertInfo> &TDV,
01182                          SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TDF,
01183                          float lambda)
01184 {
01185         //vcg::face::Pos<FaceType> ep;
01186         vcg::face::VFIterator<FaceType> ep;
01187         VertexIterator vi;
01188         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
01189         {
01190                 CoordType ErrGrad=CoordType(0,0,0);
01191 
01192                 ep.f=(*vi).VFp();
01193                 ep.z=(*vi).VFi();
01194                 while (!ep.End())
01195                 {
01196                         ErrGrad+=FaceErrorGrad(ep.f->V(ep.z)->P(),ep.f->V1(ep.z)->P(),ep.f->V2(ep.z)->P(),TDF[ep.f].m);
01197                         ++ep;
01198                 }
01199                 TDV[*vi].np=(*vi).P()-ErrGrad*(ScalarType)lambda;
01200         }
01201 
01202         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
01203                 (*vi).P()=TDV[*vi].np;
01204 
01205 }
01206 /****************************************************************************************************************/
01207 
01208 
01209 
01210 static void FastFitMesh(MeshType &m,
01211                          SimpleTempData<typename MeshType::VertContainer, PDVertInfo> &TDV,
01212                          SimpleTempData<typename MeshType::FaceContainer, PDFaceInfo> &TDF,
01213                          bool OnlySelected=false)
01214 {
01215         //vcg::face::Pos<FaceType> ep;
01216         vcg::face::VFIterator<FaceType> ep;
01217         VertexIterator vi;
01218 
01219         for(vi=m.vert.begin();vi!=m.vert.end();++vi)
01220         {
01221    CoordType Sum(0,0,0);
01222    ScalarType cnt=0;
01223    VFLocalIterator ep(&*vi);
01224          for (;!ep.End();++ep)
01225                 {
01226       CoordType bc=Barycenter<FaceType>(*ep.F());
01227       Sum += ep.F()->N()*(ep.F()->N().dot(bc - (*vi).P()));
01228       ++cnt;
01229                 }
01230                 TDV[*vi].np=(*vi).P()+ Sum*(1.0/cnt);
01231         }
01232 
01233         if(OnlySelected)
01234         {
01235                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
01236                                 if((*vi).IsS()) (*vi).P()=TDV[*vi].np;
01237         }
01238         else
01239         {
01240                 for(vi=m.vert.begin();vi!=m.vert.end();++vi)
01241                 (*vi).P()=TDV[*vi].np;
01242         }
01243 }
01244 
01245 
01246 
01247 static void VertexCoordPasoDoble(MeshType &m, int step, typename MeshType::ScalarType Sigma=0, int FitStep=10, typename MeshType::ScalarType FitLambda=0.05)
01248 {
01249         SimpleTempData< typename MeshType::VertContainer, PDVertInfo> TDV(m.vert);
01250         SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face);
01251         PDVertInfo lpzv;
01252         lpzv.np=CoordType(0,0,0);
01253         PDFaceInfo lpzf;
01254         lpzf.m=CoordType(0,0,0);
01255 
01256         assert(m.HasVFTopology());
01257         m.HasVFTopology();
01258         TDV.Start(lpzv);
01259         TDF.Start(lpzf);
01260         for(int j=0;j<step;++j)
01261         {
01262 
01263                 vcg::tri::UpdateNormals<MeshType>::PerFace(m);
01264                 FaceNormalAngleThreshold(m,TDF,Sigma);
01265                 for(int k=0;k<FitStep;k++)
01266                         FitMesh(m,TDV,TDF,FitLambda);
01267         }
01268 
01269         TDF.Stop();
01270         TDV.Stop();
01271 
01272 }
01273 
01274 // The sigma parameter affect the normal smoothing step
01275 
01276 static void VertexCoordPasoDobleFast(MeshType &m, int NormalSmoothStep, typename MeshType::ScalarType Sigma=0, int FitStep=50, bool SmoothSelected =false)
01277 {
01278         PDVertInfo lpzv;
01279         lpzv.np=CoordType(0,0,0);
01280         PDFaceInfo lpzf;
01281         lpzf.m=CoordType(0,0,0);
01282 
01283         assert(HasPerVertexVFAdjacency(m) && HasPerFaceVFAdjacency(m));
01284         SimpleTempData< typename MeshType::VertContainer, PDVertInfo> TDV(m.vert,lpzv);
01285         SimpleTempData< typename MeshType::FaceContainer, PDFaceInfo> TDF(m.face,lpzf);
01286 
01287   for(int j=0;j<NormalSmoothStep;++j)
01288            FaceNormalAngleThreshold(m,TDF,Sigma);
01289 
01290   for(int j=0;j<FitStep;++j)
01291           FastFitMesh(m,TDV,TDF,SmoothSelected);
01292 }
01293 
01294 }; //end Smooth class
01295 
01296 }               // End namespace tri
01297 }               // End namespace vcg
01298 
01299 #endif //  VCG_SMOOTH