Refactoring and cleaning of the plymc surface reconstruction algorithm

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");
+        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());
+                  sprintf(errorMessage,"%sFailed Init of mesh %s\n",errorMessage,MP.MeshName(i).c_str());
-                  return;
+                  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::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 (typename MCMesh::FaceIterator fi = me.face.begin(); fi != me.face.end(); ++fi)
-            for (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;
 }
 
 

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

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
+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!!!
 
 // forward definition
 template < class VOL >
@@ -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,8 +308,8 @@ public:
     }
 
     /*
-    Data una posizione x,y,z restituisce true se tale posizione appartiene a  un blocco gia' allocato
+     * Compute the offset <lpos> inside the subblock <rpos> of voxel (x,y,z).    
-    In ogni caso mette in rpos la posizione del subbloc e in lpos la posizione all'interno del sottoblocco
+     * return true if the subblock is allocated.
      */
     bool Pos(const int &_x,const int &_y,const int &_z, int & rpos,int &lpos) const
     {
@@ -1374,6 +1358,6 @@ class VolumeIterator
     }
 
 
-
 };
+}
 #endif

vcg/complex/algorithms/create/plymc/voxel.h

@@ -22,6 +22,24 @@
 ****************************************************************************/
 #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
 {
@@ -31,8 +49,6 @@ class Voxel
   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; 		}
 
   void SetN(const Point3<scalar> &nn) 	{ n=nn;		}
@@ -60,9 +76,9 @@ class Voxel
 
   inline Voxel & operator += ( Voxel const & vx)
   {
+    assert(!b);
     if(cnt==0)
     {
-             assert(!b);
       v=vx.v;
       q=vx.q;
       n=vx.n;
@@ -71,7 +87,6 @@ class Voxel
     }
     else
     {
-         assert(!b);
       v+=vx.v;
       q+=vx.q;
       n+=vx.n;
@@ -191,5 +206,5 @@ public:
 private:
   Point3f c;
 };
-
+}
 #endif