diff --git a/vcg/complex/algorithms/create/plymc/plymc.h b/vcg/complex/algorithms/create/plymc/plymc.h
index ed43bfe7..2f3315d3 100644
--- a/vcg/complex/algorithms/create/plymc/plymc.h
+++ b/vcg/complex/algorithms/create/plymc/plymc.h
@@ -35,23 +35,13 @@
 #include <float.h>
 #include <math.h>
 
-#include <locale>
-#include <iostream>
-
-#include <list>
-#include <vcg/space/index/grid_static_ptr.h>
 #include <vcg/complex/complex.h>
 
-#include <vcg/complex/algorithms/update/position.h>
-#include <vcg/complex/algorithms/update/normal.h>
-#include <vcg/complex/algorithms/update/quality.h>
-#include <vcg/complex/algorithms/update/topology.h>
 #include <vcg/math/histogram.h>
-#include <vcg/complex/algorithms/clean.h>
 #include <vcg/complex/algorithms/geodesic.h>
 #include <wrap/io_trimesh/import.h>
 #include <wrap/io_trimesh/export_ply.h>
-#include <wrap/ply/plystuff.h>
+//#include <wrap/ply/plystuff.h>
 
 #include <vcg/complex/algorithms/create/marching_cubes.h>
 #include <vcg/complex/algorithms/create/mc_trivial_walker.h>
@@ -61,7 +51,6 @@
 #include <vcg/complex/algorithms/local_optimization/tri_edge_collapse.h>
 #include <vcg/complex/algorithms/local_optimization/tri_edge_collapse_quadric.h>
 
-//#include <vcg/simplex/edge/base.h>
 #include <stdarg.h>
 #include "volume.h"
 #include "tri_edge_collapse_mc.h"
@@ -73,6 +62,15 @@ template<class MeshType>
 void MCSimplify( MeshType &m, float perc, bool preserveBB=true, vcg::CallBackPos *cb=0);
 
 
+/** Surface Reconstruction
+ * 
+ *  To allow the managment of a very large set of meshes to be merged,
+ *  it is templated on a MeshProvider class that is able to provide the meshes to be merged. 
+ *  IT is the surface reconstrction algorithm that have been used for a long time inside the ISTI-Visual Computer Lab.
+ *  It is mostly a variant of the Curless et al. e.g. a volumetric approach with some original weighting schemes,"
+ *  a different expansion rule, and another approach to hole filling through volume dilation/relaxations.
+ */
+
 template < class SMesh, class MeshProvider>
 class PlyMC
 {
@@ -175,6 +173,7 @@ public:
   MeshProvider MP;
   Parameter p;
   Volume<Voxelf> VV;
+  char errorMessage[1024];
 
 /// PLYMC Methods
 
@@ -192,21 +191,36 @@ public:
     {
       if(!(loadmask & tri::io::Mask::IOM_VERTNORMAL))
       {
-        printf("Error, pointset MUST have normals");
-        return false;
+        if(m.FN()==0)
+        {
+          sprintf(errorMessage,"%sError: mesh has not per vertex normals\n",errorMessage);
+          return false;
+        }
+        else
+        {          
+          tri::Clean<SMesh>::RemoveUnreferencedVertex(m);
+          tri::Allocator<SMesh>::CompactEveryVector(m);
+          tri::UpdateNormal<SMesh>::PerVertexNormalizedPerFaceNormalized(m);          
+        }
       }
+      tri::UpdateNormal<SMesh>::NormalizePerVertex(m);    
+      int badNormalCnt=0;
       for(SVertexIterator vi=m.vert.begin(); vi!=m.vert.end();++vi)
         if(math::Abs(SquaredNorm((*vi).N())-1.0)>0.0001)
         {
-          printf("Error: mesh has not per vertex normalized normals\n");
+          badNormalCnt++;
+          tri::Allocator<SMesh>::DeleteVertex(m,*vi);
+        }
+      tri::Allocator<SMesh>::CompactEveryVector(m);      
+       if(badNormalCnt > m.VN()/10)
+        {
+          sprintf(errorMessage,"%sError: mesh has null normals\n",errorMessage);
           return false;
         }
-
+      
       if(!(loadmask & tri::io::Mask::IOM_VERTQUALITY))
         tri::UpdateQuality<SMesh>::VertexConstant(m,0);
       tri::UpdateNormal<SMesh>::PerVertexMatrix(m,Tr);
-      //if(!(loadmask & tri::io::Mask::IOM_VERTCOLOR))
-      //  saveMask &= ~tri::io::Mask::IOM_VERTCOLOR;
     }
     else // processing for triangle meshes
     {
@@ -223,7 +237,6 @@ public:
           tri::UpdateTopology<SMesh>::VertexFace(m);
           tri::UpdateFlags<SMesh>::FaceBorderFromVF(m);
           tri::Geodesic<SMesh>::DistanceFromBorder(m);
-          //          tri::UpdateQuality<SMesh>::VertexGeodesicFromBorder(m);
         }
       }
 
@@ -325,8 +338,9 @@ public:
     return true;
 }
 
-void Process(vcg::CallBackPos *cb=0)
+bool Process(vcg::CallBackPos *cb=0)
 {
+  sprintf(errorMessage,"");
   printf("bbox scanning...\n"); fflush(stdout);
   Matrix44f Id; Id.SetIdentity();
   MP.InitBBox();
@@ -344,7 +358,6 @@ void Process(vcg::CallBackPos *cb=0)
 
   voxdim = fullb.max - fullb.min;
 
-  int TotAdd=0,TotMC=0,TotSav=0;
   // if kcell==0 the number of cells is computed starting from required voxel size;
   __int64 cells;
   if(p.NCell>0) cells = (__int64)(p.NCell)*(__int64)(1000);
@@ -364,6 +377,7 @@ void Process(vcg::CallBackPos *cb=0)
   }
 
 
+  int TotAdd=0,TotMC=0,TotSav=0; // partial timings counter
 
   for(p.IPos[0]=p.IPosS[0];p.IPos[0]<=p.IPosE[0];++p.IPos[0])
     for(p.IPos[1]=p.IPosS[1];p.IPos[1]<=p.IPosE[1];++p.IPos[1])
@@ -405,8 +419,8 @@ void Process(vcg::CallBackPos *cb=0)
                 res = InitMesh(*sm,MP.MeshName(i).c_str(),MP.Tr(i));
                 if(!res)
                 {
-                  printf("Failed Init of mesh %s",MP.MeshName(i).c_str());
-                  return;
+                  sprintf(errorMessage,"%sFailed Init of mesh %s\n",errorMessage,MP.MeshName(i).c_str());
+                  return false ;
                 }
               }
               res |= AddMeshToVolumeM(*sm, MP.MeshName(i),MP.W(i));
@@ -452,26 +466,20 @@ void Process(vcg::CallBackPos *cb=0)
             VV.SlicedPPM("final","__",p.SliceNum);
             VV.SlicedPPMQ("final","__",p.SliceNum);
           }
-          //MCMesh me;
-          //
           MCMesh me;
           if(res)
           {
-            typedef vcg::tri::TrivialWalker<MCMesh, Volume <Voxelf> >	Walker;
+            typedef vcg::tri::TrivialWalker<MCMesh, Volume <Voxelf> >	  Walker;
             typedef vcg::tri::MarchingCubes<MCMesh, Walker>             MarchingCubes;
-            //typedef vcg::tri::ExtendedMarchingCubes<MCMesh, Walker> ExtendedMarchingCubes;
 
             Walker walker;
             MarchingCubes	mc(me, walker);
-            Box3i currentSubBox=VV.SubPartSafe;
-            Point3i currentSubBoxRes=VV.ssz;
             /**********************/
             if(cb) cb(50,"Step 2: Marching Cube...");
             else printf("Step 2: Marching Cube...\n");
             /**********************/
-            walker.Init(VV,currentSubBox);
+            walker.SetExtractionBox(VV.SubPartSafe);
             walker.BuildMesh(me,VV,mc,0);
-            //            walker.BuildMesh(me,VV,mc,currentSubBox,currentSubBoxRes);
 
             typename MCMesh::VertexIterator vi;
             Box3f bbb; bbb.Import(VV.SubPart);
@@ -481,8 +489,7 @@ void Process(vcg::CallBackPos *cb=0)
                 vcg::tri::Allocator< MCMesh >::DeleteVertex(me,*vi);
               VV.DeInterize((*vi).P());
             }
-            typename MCMesh::FaceIterator fi;
-            for (fi = me.face.begin(); fi != me.face.end(); ++fi)
+            for (typename MCMesh::FaceIterator fi = me.face.begin(); fi != me.face.end(); ++fi)
             {
               if((*fi).V(0)->IsD() || (*fi).V(1)->IsD() || (*fi).V(2)->IsD() )
                 vcg::tri::Allocator< MCMesh >::DeleteFace(me,*fi);
@@ -526,6 +533,7 @@ void Process(vcg::CallBackPos *cb=0)
         {
           printf("----------- skipping SubBlock %2i %2i %2i ----------\n",p.IPos[0],p.IPos[1],p.IPos[2]);
         }
+  return true;
 }
 
 
diff --git a/vcg/complex/algorithms/create/plymc/simplemeshprovider.h b/vcg/complex/algorithms/create/plymc/simplemeshprovider.h
new file mode 100644
index 00000000..1ce3be53
--- /dev/null
+++ b/vcg/complex/algorithms/create/plymc/simplemeshprovider.h
@@ -0,0 +1,233 @@
+#ifndef SIMPLEMESHPROVIDER_H
+#define SIMPLEMESHPROVIDER_H
+#include "../../meshlab/alnParser.h"
+#include <list>
+#include <vector>
+#include <vcg/space/box3.h>
+#include <wrap/ply/plystuff.h>
+#include <wrap/io_trimesh/import.h>
+
+/*
+ * A mesh provider class has the simpler role of passing the set of meshes to be merged to the surface reconstrcution algorithm. 
+ * The only reason for this abstraction is that, plymc can work in a out-of-core way and the loading of the needed range maps can be optimized with a high level caching system. 
+ */
+template<class TriMeshType>
+class MinimalMeshProvider
+{
+private:
+ 
+  std::vector< std::string > nameVec;  
+  std::vector< TriMeshType * > meshVec;
+  std::vector<vcg::Matrix44f> trVec;
+  std::vector<float> weightVec;		        // weight tot be applied to each mesh.
+  vcg::Box3f fullBBox;
+
+public:
+  bool Find(const std::string &name,  TriMeshType * &sm) 
+  {
+    for(int i=0;i<nameVec.size();++i) 
+      if(nameVec[i]==name) {
+        sm=meshVec[i]; 
+        return true;
+      } 
+    sm=0; return false;
+  }
+  
+};
+
+/**
+ * Cache based Loading of  meshes to avoid reloading an processing of the same mesh multiple times. 
+ * 
+ */
+
+namespace vcg {
+
+template<class TriMeshType>
+class MeshCache
+{
+  class Pair
+  {
+  public:
+    Pair(){used=0;}
+    TriMeshType *M;
+    std::string Name;
+    int used; // 'data' dell'ultimo accesso. si butta fuori quello lru
+  };
+  
+  std::list<Pair> MV;
+  
+public:
+  void clear();
+  
+  MeshCache() {MeshCacheSize=6;}
+  ~MeshCache() {
+    typename std::list<Pair>::iterator mi;
+    for(mi=MV.begin();mi!=MV.end();++mi)
+      delete (*mi).M;
+  }
+  
+
+  /**
+   * @brief Find load a mesh form the cache if it is in or from the disk otherwise
+   * @param name what mesh to find
+   * @param sm the pointer loaded mesh
+   * @return true if the mesh was already in cache
+   * 
+   */  
+  bool Find(const std::string &name,  TriMeshType * &sm)
+  {
+    typename std::list<Pair>::iterator mi;
+    typename std::list<Pair>::iterator oldest; // quello che e' piu' tempo che non viene acceduto.
+    int last;
+    
+    last = std::numeric_limits<int>::max();
+    oldest = MV.begin();
+    
+    for(mi=MV.begin();mi!=MV.end();++mi)
+    {
+      if((*mi).used<last)
+      {
+        last=(*mi).used;
+        oldest=mi;
+      }
+      if((*mi).Name==name) {
+        sm=(*mi).M;
+        (*mi).used++;
+        return true;
+      }
+    }
+    
+    // we have not found the requested mesh
+    // either allocate a new mesh or give back a previous mesh.
+    
+    if(MV.size()>MeshCacheSize)	{
+      sm=(*oldest).M;
+      (*oldest).used=0;
+      (*oldest).Name=name;
+    }	else	{
+      MV.push_back(Pair());
+      MV.back().Name=name;
+      MV.back().M=new TriMeshType();
+      sm=MV.back().M;
+    }
+    return false;
+  }
+  
+  
+  size_t MeshCacheSize;
+  size_t size() const {return MV.size();}
+};
+
+template<class TriMeshType>
+class SimpleMeshProvider
+{
+private:
+  std::vector< std::string > meshnames;
+  std::vector<vcg::Matrix44f> TrV;
+  std::vector<float> WV;		        // weight tot be applied to each mesh.
+  std::vector<vcg::Box3f> BBV;	    // bbox of the transformed meshes..
+  vcg::Box3f fullBBox;
+  MeshCache<TriMeshType> MC;
+  
+public:
+  
+  int size() {return meshnames.size();}
+  
+  int getCacheSize() {return MC.MeshCacheSize;}
+  int setCacheSize(size_t newsize)
+  {
+    if(newsize == MC.MeshCacheSize)
+      return MC.MeshCacheSize;
+    if(newsize <= 0)
+      return MC.MeshCacheSize;
+    
+    MC.MeshCacheSize = newsize;
+    return newsize;
+  }
+  
+  bool openALN (const char* alnName)
+  {
+    vector<RangeMap> rmaps;
+    ALNParser::ParseALN(rmaps, alnName);
+    
+    for(size_t i=0; i<rmaps.size(); i++)
+      AddSingleMesh(rmaps[i].filename.c_str(), rmaps[i].trasformation, rmaps[i].quality);
+    
+    return true;
+  }
+  
+  bool AddSingleMesh(const char* meshName, const Matrix44f &tr= Matrix44f::Identity(), float meshWeight=1)
+  {
+    assert(WV.size()==meshnames.size() && TrV.size() == WV.size());
+    TrV.push_back(tr);
+    meshnames.push_back(meshName);
+    WV.push_back(meshWeight);
+    BBV.push_back(Box3f());
+    return true;
+  }
+    
+  vcg::Box3f bb(int i) {return BBV[i];}
+  vcg::Box3f fullBB(){ return fullBBox;}
+  vcg::Matrix44f Tr(int i) const  {return TrV[i];}
+  std::string MeshName(int i) const {return meshnames[i];}
+  float W(int i) const {return WV[i];}
+  
+  void Clear()
+  {
+    meshnames.clear();
+    TrV.clear();
+    WV.clear();
+    BBV.clear();
+    fullBBox.SetNull();
+    MC.clear();
+  }
+  
+  bool Find(int i, TriMeshType * &sm)
+  {
+    return MC.Find(meshnames[i],sm);
+  }
+  
+  bool InitBBox()
+  {
+    fullBBox.SetNull();    
+    for(int i=0;i<int(meshnames.size());++i)
+    {
+      bool ret;      
+      printf("bbox scanning %4i/%i [%16s]      \r",i+1,(int)meshnames.size(), meshnames[i].c_str());
+      if(tri::io::Importer<TriMeshType>::FileExtension(meshnames[i],"PLY") || tri::io::Importer<TriMeshType>::FileExtension(meshnames[i],"ply"))
+      {
+        ret=ply::ScanBBox(meshnames[i].c_str(),BBV[i],TrV[i],true,0);
+      }
+      else
+      {
+        printf("Trying to import a non-ply file %s\n",meshnames[i].c_str());fflush(stdout);
+        TriMeshType m;
+        ret = (tri::io::Importer<TriMeshType>::Open(m,meshnames[i].c_str()) == tri::io::Importer<TriMeshType>::E_NOERROR);
+        tri::UpdatePosition<TriMeshType>::Matrix(m,TrV[i]);
+        tri::UpdateBounding<TriMeshType>::Box(m);
+        BBV[i].Import(m.bbox);
+      }
+      if( ! ret)
+      {
+        printf("\n\nwarning:\n file '%s' not found\n",meshnames[i].c_str());fflush(stdout);
+        continue;
+      }
+      fullBBox.Add(BBV[i]);
+    }
+    return true;
+  }
+  
+};
+
+class SVertex;
+class SFace;
+class SUsedTypes: public vcg::UsedTypes < vcg::Use<SVertex>::AsVertexType,
+    vcg::Use<SFace  >::AsFaceType >{};
+
+class SVertex     : public Vertex< SUsedTypes, vertex::Coord3f, vertex::Normal3f,vertex::VFAdj, vertex::BitFlags, vertex::Color4b, vertex::Qualityf>{};
+class SFace       : public Face< SUsedTypes, face::VertexRef, face::Normal3f,face::Qualityf, face::VFAdj, face::BitFlags> {};
+class SMesh       : public tri::TriMesh< std::vector< SVertex>, std::vector< SFace > > {};
+
+}
+
+#endif // SIMPLEMESHPROVIDER_H
diff --git a/vcg/complex/algorithms/create/plymc/volume.h b/vcg/complex/algorithms/create/plymc/volume.h
index d543c32e..0a9ba7fc 100644
--- a/vcg/complex/algorithms/create/plymc/volume.h
+++ b/vcg/complex/algorithms/create/plymc/volume.h
@@ -24,26 +24,10 @@
 #ifndef __VOLUME_H__
 #define __VOLUME_H__
 
-#ifdef __MINGW32__
-#define _int64 __int64
-#endif
-
 #include "voxel.h"
-#include "svoxel.h"
-#include <vector>
 #include <vcg/space/index/grid_static_ptr.h>
 
-//#define BLOCKSIDE() 8
-
-// Stato di un voxel
-
-// B() dice se ci sono dati in uno stadio usabile.
-// Cnt() dice quanti ce ne sono stati sommati (per la normalizzazione)
-
-// b==false cnt==0 totalmente non inzializzato (Zero)
-// b==false cnt >0 da normalizzare
-// b==true  cnt==0 gia' normalizzato
-// b==true  cnt >0 Errore!!!
+namespace vcg {
 
 // forward definition
 template < class VOL >
@@ -67,7 +51,7 @@ const char *SFormat( const char * f, ... )
 template<class VOX_TYPE, class SCALAR_TYPE=float>
 class Volume {
 public:
-    typedef SCALAR_TYPE scalar;
+  typedef SCALAR_TYPE scalar;
   typedef Point3<scalar> Point3x;
   typedef Box3<scalar> Box3x;
 
@@ -172,7 +156,7 @@ bool Verbose; // se true stampa un sacco di info in piu su logfp;
         for(size_t i=0;i<rv.size();++i)
             rv[i].resize(0,VOX_TYPE::Zero());
         SetDim(bb);
-    };
+    }
 
 private:
 
@@ -324,16 +308,16 @@ public:
     }
 
     /*
-    Data una posizione x,y,z restituisce true se tale posizione appartiene a  un blocco gia' allocato
-    In ogni caso mette in rpos la posizione del subbloc e in lpos la posizione all'interno del sottoblocco
-    */
+     * Compute the offset <lpos> inside the subblock <rpos> of voxel (x,y,z).    
+     * return true if the subblock is allocated.
+     */
     bool Pos(const int &_x,const int &_y,const int &_z, int & rpos,int &lpos) const
     {
         int x=_x-SubPartSafe.min[0];		int y=_y-SubPartSafe.min[1];		int z=_z-SubPartSafe.min[2];
 
         assert(_x>=SubPartSafe.min[0] && _x<SubPartSafe.max[0] &&
-                 _y>=SubPartSafe.min[1] && _y<SubPartSafe.max[1] &&
-                     _z>=SubPartSafe.min[2] && _z<SubPartSafe.max[2]);
+               _y>=SubPartSafe.min[1] && _y<SubPartSafe.max[1] &&
+               _z>=SubPartSafe.min[2] && _z<SubPartSafe.max[2]);
 
     //	assert(x>=0 && x<sz[0] && y>=0 && y<sz[1] && z>=0 && z<sz[2]);
 
@@ -1374,6 +1358,6 @@ class VolumeIterator
     }
 
 
-
 };
+}
 #endif
diff --git a/vcg/complex/algorithms/create/plymc/voxel.h b/vcg/complex/algorithms/create/plymc/voxel.h
index 3dc18642..92f18dfa 100644
--- a/vcg/complex/algorithms/create/plymc/voxel.h
+++ b/vcg/complex/algorithms/create/plymc/voxel.h
@@ -8,7 +8,7 @@
 *                                                                    \      *
 * All rights reserved.                                                      *
 *                                                                           *
-* This program is free software; you can redistribute it and/or modify      *   
+* This program is free software; you can redistribute it and/or modify      *
 * it under the terms of the GNU General Public License as published by      *
 * the Free Software Foundation; either version 2 of the License, or         *
 * (at your option) any later version.                                       *
@@ -22,105 +22,120 @@
 ****************************************************************************/
 #ifndef __VOXEL_H__
 #define __VOXEL_H__
+
+namespace vcg {
+
+// Stato di un voxel
+
+// B() dice se ci sono dati in uno stadio usabile.
+// Cnt() dice quanti ce ne sono stati sommati (per la normalizzazione)
+
+// b==false cnt==0 totalmente non inzializzato (Zero)
+// b==false cnt >0 da normalizzare
+// b==true  cnt==0 gia' normalizzato
+// b==true  cnt >0 Errore!!!
+
+
+/**
+ *
+ */
+
 template<class SCALAR_TYPE=float>
 class Voxel
 {
-    public:
-    typedef SCALAR_TYPE scalar;
+public:
+  typedef SCALAR_TYPE scalar;
 
-    Voxel(SCALAR_TYPE vv, bool bb, Point3<scalar> nn, short _cnt) {v=vv;b=bb;n=nn;cnt=_cnt;}
-    Voxel(SCALAR_TYPE vv, Point3<scalar> nn, scalar qq) {v=vv;b=true;n=nn;cnt=0;q=qq;}
+  Voxel(SCALAR_TYPE vv, bool bb, Point3<scalar> nn, short _cnt) {v=vv;b=bb;n=nn;cnt=_cnt;}
+  Voxel(SCALAR_TYPE vv, Point3<scalar> nn, scalar qq) {v=vv;b=true;n=nn;cnt=0;q=qq;}
 
-    const scalar &N(const int i) const { return n[i]; }
+  const Point3<scalar> &N() const  	{  return n; 		}
 
-    const Point3<scalar> &N() const  	{  return n; 		}
+  void SetN(const Point3<scalar> &nn) 	{ n=nn;		}
+  const scalar &V() const 		{ 	 return v; 		}
 
-    void SetN(const Point3<scalar> &nn) 	{ n=nn;		}
-    const scalar &V() const 		{ 	 return v; 		}
+  void SetV(const scalar &vv) 		{		 v=vv;		}
 
-    void SetV(const scalar &vv) 		{		 v=vv;		}
+  const scalar &Q() const		{		 return q;		}
 
-    const scalar &Q() const		{		 return q;		}
-
-    void SetQ(const scalar &qq) 		{		 q=qq;		}
+  void SetQ(const scalar &qq) 		{		 q=qq;		}
 
 
-    bool B() const {return b;};
-    void SetB(bool val) {b=val;}
-    int Cnt() const {return cnt;}
-    void SetCnt(int val) {cnt=val;}
-    inline void Blend( Voxel const & vx, scalar w)
+  bool B() const {return b;};
+  void SetB(bool val) {b=val;}
+  int Cnt() const {return cnt;}
+  void SetCnt(int val) {cnt=val;}
+  inline void Blend( Voxel const & vx, scalar w)
+  {
+    float w1=1.0-w;
+    v=v*w1+vx.v*w;
+    q=q*w1+vx.q*w;
+    n=n*w1+vx.n*w;
+    //return *this;
+  }
+
+  inline Voxel & operator += ( Voxel const & vx)
+  {
+    assert(!b);
+    if(cnt==0)
     {
-        float w1=1.0-w;
-        v=v*w1+vx.v*w;
-        q=q*w1+vx.q*w;
-        n=n*w1+vx.n*w;
-        //return *this;
+      v=vx.v;
+      q=vx.q;
+      n=vx.n;
+      cnt=1;
+      b=false;
     }
-
-    inline Voxel & operator += ( Voxel const & vx)
+    else
     {
-     if(cnt==0)
-         {
-             assert(!b);
-            v=vx.v;
-            q=vx.q;
-            n=vx.n;
-            cnt=1;
-            b=false;
-         }
-     else
-        {
-         assert(!b);
-         v+=vx.v;
-         q+=vx.q;
-         n+=vx.n;
-         ++cnt;
-         }
-     return *this;
+      v+=vx.v;
+      q+=vx.q;
+      n+=vx.n;
+      ++cnt;
     }
+    return *this;
+  }
 
-    inline bool Normalize(int thr)
+  inline bool Normalize(int thr)
+  {
+    assert(cnt>0);
+    assert(!B());
+    if(cnt<thr)
     {
-     assert(cnt>0);
-     assert(!B());
-     if(cnt<thr)
-     {
-         (*this) = Zero();
-         return false;
-     }
-     v/=cnt;
-     q/=cnt;
-     n/=cnt;
-     cnt=0;
-     b=true;
+      (*this) = Zero();
+      return false;
+    }
+    v/=cnt;
+    q/=cnt;
+    n/=cnt;
+    cnt=0;
+    b=true;
     return true;
-    }
+  }
 
-    static const Voxel &Zero() {
-        static Voxel tt(0,false,Point3f(0,0,0),0);
-        return tt;
-    }
-    void Merge(const Voxel &VOX)
-        {
-            v=(v*q+VOX.Q()*VOX.v)/(q+VOX.Q());
-            n=(n*q+VOX.n*VOX.Q())/(q+VOX.Q());
-            q=q+VOX.Q();
-        }
+  static const Voxel &Zero() {
+    static Voxel tt(0,false,Point3f(0,0,0),0);
+    return tt;
+  }
+  void Merge(const Voxel &VOX)
+  {
+    v=(v*q+VOX.Q()*VOX.v)/(q+VOX.Q());
+    n=(n*q+VOX.n*VOX.Q())/(q+VOX.Q());
+    q=q+VOX.Q();
+  }
 
-    void Set(const Voxel &VOX)
-        {
-            v=VOX.v;
-            n=VOX.n;
-            q=VOX.q;
-        }
+  void Set(const Voxel &VOX)
+  {
+    v=VOX.v;
+    n=VOX.n;
+    q=VOX.q;
+  }
 
-    protected:
-    bool b;
-    short cnt;
-    scalar v;
-    scalar q;
-    Point3<SCALAR_TYPE> n;
+protected:
+  bool b;
+  short cnt;
+  scalar v;
+  scalar q;
+  Point3<SCALAR_TYPE> n;
 
 };
 
@@ -129,67 +144,67 @@ class Voxelfc :public Voxel<float>
 {
 public:
 
-    Voxelfc(float vv, bool bb, Point3f nn, short _cnt) :Voxel<float>(vv,bb,nn,_cnt){}
-    Voxelfc(float vv, Point3f nn, scalar qq) :Voxel<float>(vv,nn,qq) {}
-    Voxelfc(float vv, Point3f nn, scalar qq,Color4b cc) :Voxel<float>(vv,nn,qq)
-    {
-        c[0]=cc[0];
-        c[1]=cc[1];
-        c[2]=cc[2];
-    }
+  Voxelfc(float vv, bool bb, Point3f nn, short _cnt) :Voxel<float>(vv,bb,nn,_cnt){}
+  Voxelfc(float vv, Point3f nn, scalar qq) :Voxel<float>(vv,nn,qq) {}
+  Voxelfc(float vv, Point3f nn, scalar qq,Color4b cc) :Voxel<float>(vv,nn,qq)
+  {
+    c[0]=cc[0];
+    c[1]=cc[1];
+    c[2]=cc[2];
+  }
 
-    inline bool Normalize(int thr)
-    {
-        if(cnt>=thr)  c/=cnt;
-        return Voxel<float>::Normalize(thr);
-    }
+  inline bool Normalize(int thr)
+  {
+    if(cnt>=thr)  c/=cnt;
+    return Voxel<float>::Normalize(thr);
+  }
 
-    static const Voxelfc &Zero() {
-        static Voxelfc tt(0,false,Point3f(0,0,0),0);
-        return tt;
-    }
+  static const Voxelfc &Zero() {
+    static Voxelfc tt(0,false,Point3f(0,0,0),0);
+    return tt;
+  }
 
-    void Merge(const Voxelfc &VOX)
-    {
-        c=( c*q + VOX.C()*VOX.Q() )/(q+VOX.Q());
-        Voxel<float>::Merge(VOX);
-    }
+  void Merge(const Voxelfc &VOX)
+  {
+    c=( c*q + VOX.C()*VOX.Q() )/(q+VOX.Q());
+    Voxel<float>::Merge(VOX);
+  }
 
-    void Set(const Voxelfc &VOX)
-    {
-        Voxel<float>::Set(VOX);
-        c=VOX.c;
-    }
+  void Set(const Voxelfc &VOX)
+  {
+    Voxel<float>::Set(VOX);
+    c=VOX.c;
+  }
 
-    const float &C(const int i) const { return c[i]; }
-    const Point3f &C() const  	{  return c; 		}
-    void SetC(const Point3f &cc) 	{ c=cc;		}
-    Color4b C4b() const
-        {
-            static Color4b cc;
-            cc=Color4b(c[0],c[1],c[2],255);
-            return cc;
-        }
-    inline void Blend( Voxelfc const & vx, scalar w)
-    {
-        float w1=1.0-w;
-        v=v*w1+vx.v*w;
-        q=q*w1+vx.q*w;
-        n=n*w1+vx.n*w;
-        c=c*w1+vx.c*w;
-        //return *this;
-    }
+  const float &C(const int i) const { return c[i]; }
+  const Point3f &C() const  	{  return c; 		}
+  void SetC(const Point3f &cc) 	{ c=cc;		}
+  Color4b C4b() const
+  {
+    static Color4b cc;
+    cc=Color4b(c[0],c[1],c[2],255);
+    return cc;
+  }
+  inline void Blend( Voxelfc const & vx, scalar w)
+  {
+    float w1=1.0-w;
+    v=v*w1+vx.v*w;
+    q=q*w1+vx.q*w;
+    n=n*w1+vx.n*w;
+    c=c*w1+vx.c*w;
+    //return *this;
+  }
 
-    inline Voxelfc & operator += ( Voxelfc const & vx)
-    {
-            Voxel<float>::operator +=(vx);
-            if(cnt==1)	c =vx.c;
-                         else c+=vx.c;
-            return *this;
-    }
+  inline Voxelfc & operator += ( Voxelfc const & vx)
+  {
+    Voxel<float>::operator +=(vx);
+    if(cnt==1)	c =vx.c;
+    else c+=vx.c;
+    return *this;
+  }
 
 private:
- Point3f c;
+  Point3f c;
 };
-
+}
 #endif