From a6ba20c33875ac69643afd13cd3d66cd82b0c582 Mon Sep 17 00:00:00 2001
From: cignoni <paolo.cignoni@isti.cnr.it>
Date: Tue, 29 Dec 2015 07:22:13 +0000
Subject: [PATCH] First version of the Curve On Manifold managment class.

---
 vcg/complex/algorithms/curve_on_manifold.h | 817 +++++++++++++++++++++
 1 file changed, 817 insertions(+)
 create mode 100644 vcg/complex/algorithms/curve_on_manifold.h

diff --git a/vcg/complex/algorithms/curve_on_manifold.h b/vcg/complex/algorithms/curve_on_manifold.h
new file mode 100644
index 00000000..73e061b0
--- /dev/null
+++ b/vcg/complex/algorithms/curve_on_manifold.h
@@ -0,0 +1,817 @@
+/****************************************************************************
+* VCGLib                                                            o o     *
+* Visual and Computer Graphics Library                            o     o   *
+*                                                                _   O  _   *
+* Copyright(C) 2004                                                \/)\/    *
+* Visual Computing Lab                                            /\/|      *
+* ISTI - Italian National Research Council                           |      *
+*                                                                    \      *
+* All rights reserved.                                                      *
+*                                                                           *
+* 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.                                       *
+*                                                                           *
+* This program is distributed in the hope that it will be useful,           *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of            *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the             *
+* GNU General Public License (http://www.gnu.org/licenses/gpl.txt)          *
+* for more details.                                                         *
+*                                                                           *
+****************************************************************************/
+#ifndef __VCGLIB_CURVE_ON_SURF_H
+#define __VCGLIB_CURVE_ON_SURF_H
+
+#include<vcg/complex/complex.h>
+#include<vcg/simplex/face/topology.h>
+#include<vcg/complex/algorithms/update/topology.h>
+#include<vcg/complex/algorithms/update/color.h>
+#include<vcg/complex/algorithms/update/normal.h>
+#include<vcg/complex/algorithms/update/quality.h>
+#include<vcg/complex/algorithms/clean.h>
+#include<vcg/complex/algorithms/refine.h>
+#include<vcg/complex/algorithms/create/platonic.h>
+#include <vcg/space/index/grid_static_ptr.h>
+#include <vcg/space/index/kdtree/kdtree.h>
+#include <vcg/math/histogram.h>
+#include<vcg/space/distance3.h>
+#include<eigenlib/Eigen/Core>
+
+namespace vcg {
+namespace tri {
+/// \ingroup trimesh
+/// \brief A class for managing curves on a 2manifold.
+/**
+  This class is used to project/simplify/smooth polylines over a triangulated surface. 
+  
+*/
+
+template <class MeshType>
+class CoM
+{
+public:
+  typedef typename MeshType::ScalarType     ScalarType;
+  typedef typename MeshType::VertexType     VertexType;
+  typedef typename MeshType::VertexPointer  VertexPointer;
+  typedef typename MeshType::VertexIterator VertexIterator;
+  typedef typename MeshType::EdgeIterator   EdgeIterator;
+  typedef typename MeshType::EdgeType       EdgeType;
+  typedef typename MeshType::FaceType       FaceType;
+  typedef typename MeshType::FacePointer    FacePointer;
+  typedef typename MeshType::FaceIterator   FaceIterator;
+  typedef Box3  <ScalarType>               Box3Type;
+  typedef typename vcg::GridStaticPtr<FaceType, ScalarType> MeshGrid;  
+  typedef typename vcg::GridStaticPtr<EdgeType, ScalarType> EdgeGrid;
+  typedef typename face::Pos<FaceType> PosType;
+  
+  class Param 
+  {
+  public:
+    
+    ScalarType surfDistThr; // Distance between surface and curve; used in simplify and refine
+    ScalarType polyDistThr; // Distance between the 
+    ScalarType minRefEdgeLen;  // Minimal admitted Edge Lenght (used in refine: never make edge shorther than this value) 
+    ScalarType maxSimpEdgeLen; // Minimal admitted Edge Lenght (used in simplify: never make edges longer than this value) 
+    ScalarType maxSmoothDelta; // The maximum movement that is admitted during smoothing.
+    
+    Param(MeshType &m) { Default(m);}
+    
+    void Default(MeshType &m)
+    {
+      surfDistThr   = m.bbox.Diag()/10000.0;
+      polyDistThr   = m.bbox.Diag()/1000.0;
+      minRefEdgeLen    = m.bbox.Diag()/2000.0;
+      maxSimpEdgeLen    = m.bbox.Diag()/1000.0;
+      maxSmoothDelta = m.bbox.Diag()/100.0;
+    }
+    void Dump() const
+    {
+      printf("surfDistThr    = %6.3f\n",surfDistThr   );
+      printf("polyDistThr    = %6.3f\n",polyDistThr   );
+      printf("minEdgeLen     = %6.3f\n",minRefEdgeLen    );
+      printf("maxSmoothDelta = %6.3f\n",maxSmoothDelta);
+    }
+  };
+  
+  
+  
+  // The Data Members
+  
+  MeshType &base; 
+  MeshGrid uniformGrid;
+  Param par; 
+  CoM(MeshType &_m) :base(_m),par(_m){}
+ 
+  
+  //******** CUT TREE GENERATION 
+
+  int countNonVisitedEdges(VertexType * vp, EdgeType * &ep)
+  {
+    std::vector<EdgeType *> starVec;
+    edge::VEStarVE(&*vp,starVec);
+  
+    int cnt =0;
+    for(size_t i=0;i<starVec.size();++i)
+    {
+      if(!starVec[i]->IsV())
+      {
+        cnt++;
+        ep = starVec[i];
+      }
+    }
+  
+    return cnt;
+  }
+  
+  void Retract(MeshType &t)
+  {
+    tri::Clean<MeshType>::RemoveDuplicateVertex(t);
+    printf("Retracting a mesh of %i edges and %i vertices\n",t.en,t.vn);
+    tri::UpdateTopology<MeshType>::VertexEdge(t);
+  
+    std::stack<VertexType *> vertStack;
+  
+    for(VertexIterator vi=t.vert.begin();vi!=t.vert.end();++vi)
+    {
+      std::vector<EdgeType *> starVec;
+      edge::VEStarVE(&*vi,starVec);
+      if(starVec.size()==1)
+        vertStack.push(&*vi);
+    }
+  
+    tri::UpdateFlags<MeshType>::EdgeClearV(t);
+    tri::UpdateFlags<MeshType>::VertexClearV(t);
+  
+    while(!vertStack.empty())
+    {
+      VertexType *vp = vertStack.top();
+      vertStack.pop();
+      vp->C()=Color4b::Blue;
+      EdgeType *ep=0;
+  
+      int eCnt =  countNonVisitedEdges(vp,ep);
+      if(eCnt==1) // We have only one non visited edge over vp
+      {
+        assert(!ep->IsV());
+        ep->SetV();
+        VertexType *otherVertP;
+        if(ep->V(0)==vp) otherVertP = ep->V(1);
+        else otherVertP = ep->V(0);
+        vertStack.push(otherVertP);
+      }
+    }
+    
+    for(size_t i =0; i<t.edge.size();++i)
+      if(t.edge[i].IsV()) tri::Allocator<MeshType>::DeleteEdge(t,t.edge[i]);
+  
+    tri::Clean<MeshType>::RemoveUnreferencedVertex(t);
+    tri::Allocator<MeshType>::CompactEveryVector(t);
+  
+  }
+  
+  
+  // 
+  void BuildVisitTree(MeshType &dualMesh)
+  {
+    tri::UpdateTopology<MeshType>::FaceFace(base);
+    tri::UpdateFlags<MeshType>::FaceClearV(base);
+    
+    std::vector<face::Pos<FaceType> > visitStack; // the stack contain the pos on the 'starting' face. 
+    
+    
+    MeshType treeMesh; // this mesh will contain the set of the non traversed edge 
+    
+    base.face[0].SetV();
+    for(int i=0;i<3;++i)
+      visitStack.push_back(PosType(&(base.face[0]),i,base.face[0].V(i)));
+  
+    int cnt=1;
+    
+    while(!visitStack.empty())
+    {
+      std::swap(visitStack.back(),visitStack[rand()%visitStack.size()]);
+      PosType c = visitStack.back();
+      visitStack.pop_back();
+      assert(c.F()->IsV());
+      c.F()->C() = Color4b::ColorRamp(0,base.fn,cnt);
+      c.FlipF();
+      if(!c.F()->IsV())
+      {
+        tri::Allocator<MeshType>::AddEdge(treeMesh,Barycenter(*(c.FFlip())),Barycenter(*(c.F())));
+        ++cnt;
+        c.F()->SetV();
+        c.FlipE();c.FlipV();
+        visitStack.push_back(c);
+        c.FlipE();c.FlipV();
+        visitStack.push_back(c);
+      }
+      else
+      {
+        tri::Allocator<MeshType>::AddEdge(dualMesh,c.V()->P(),c.VFlip()->P());
+      }
+    }
+    assert(cnt==base.fn);
+    
+    Retract(dualMesh);
+}  
+  
+  // Given an edge of a mesh, supposedly intersecting the polyline, 
+  // we search on it the closest point to the polyline
+  static float MinDistOnEdge(VertexType *v0,VertexType *v1, EdgeGrid &edgeGrid, MeshType &poly, Point3f &closestPoint)
+  {
+    float minPolyDist = std::numeric_limits<ScalarType>::max();
+    const float sampleNum = 10;
+    const float maxDist = poly.bbox.Diag()/10.0;
+    for(float k = 0;k<sampleNum+1;++k)
+    {
+      float polyDist;
+      Point3f closestPPoly;
+      Point3f samplePnt = (v0->P()*k +v1->P()*(sampleNum-k))/sampleNum;          
+      
+      EdgeType *cep = vcg::tri::GetClosestEdgeBase(poly,edgeGrid,samplePnt,maxDist,polyDist,closestPPoly);        
+      
+      if(polyDist < minPolyDist)
+      {
+        minPolyDist = polyDist;
+        closestPoint = closestPPoly;
+      }
+    }
+    return minPolyDist;    
+  }
+  
+  class QualitySign
+  {
+  public:
+      EdgeGrid &edgeGrid;
+      MeshType &poly;
+      Param &par;
+      QualitySign(EdgeGrid &_e,MeshType &_poly, Param &_par):edgeGrid(_e),poly(_poly),par(_par) {};
+      bool operator()(face::Pos<FaceType> ep) const
+      {
+        VertexType *v0 = ep.V();
+        VertexType *v1 = ep.VFlip();
+        if(v0->Q() * v1->Q() < 0)
+        { 
+          //Point3f pp = CoS::QLerp(v0,v1);
+          Point3f closestP;
+          float minDist = MinDistOnEdge(v0,v1,edgeGrid,poly,closestP);
+          if(minDist < par.polyDistThr) return true;
+        }
+        return false;
+      }
+  };
+  
+  static Point3f QLerp(VertexType *v0, VertexType *v1)
+  {
+    
+    float qSum = fabs(v0->Q())+fabs(v1->Q());      
+    float w0 = (qSum - fabs(v0->Q()))/qSum;
+    float w1 = (qSum - fabs(v1->Q()))/qSum;      
+    return v0->P()*w0 + v1->P()*w1;      
+  }
+  
+  struct QualitySignSplit : public std::unary_function<face::Pos<FaceType> ,  Point3f>
+  {
+    EdgeGrid &edgeGrid;
+    MeshType &poly;
+    Param &par;
+    QualitySignSplit(EdgeGrid &_e,MeshType &_p, Param &_par):edgeGrid(_e),poly(_p),par(_par) {};
+      void operator()(VertexType &nv, face::Pos<FaceType> ep)
+      {
+        VertexType *v0 = ep.V();
+        VertexType *v1 = ep.VFlip();
+        Point3f closestP;
+        float minDist = MinDistOnEdge(v0,v1,edgeGrid,poly,closestP);
+        nv.P()=closestP;        
+//        nv.P() = CoS::QLerp(v0,v1);        
+      }
+      Color4b WedgeInterp(Color4b &c0, Color4b &c1)
+      {
+          Color4b cc;
+          cc.lerp(c0,c1,0.5f);
+          return Color4b::Red;
+      }
+      TexCoord2f WedgeInterp(TexCoord2f &t0, TexCoord2f &t1)
+      {
+          TexCoord2f tmp;
+          assert(t0.n()== t1.n());
+          tmp.n()=t0.n();
+          tmp.t()=(t0.t()+t1.t())/2.0;
+          return tmp;
+      }
+  };
+  
+  void DumpPlanes(MeshType &poly, std::vector<Plane3f> &planeVec)
+  {
+    MeshType full;
+    for(int i=0;i<planeVec.size();++i)
+    {
+      float radius=edge::Length(poly.edge[i])/2.0;
+      MeshType t;
+      OrientedDisk(t,16,edge::Center(poly.edge[i]),planeVec[i].Direction(),radius);
+      tri::Append<MeshType,MeshType>::Mesh(full,t);
+    }
+    tri::io::ExporterPLY<MeshType>::Save(full,"planes.ply");  
+  }
+
+  Plane3f ComputeEdgePlane(VertexType *v0, VertexType *v1)
+  {
+    Plane3f pl;
+    Point3f mid = (v0->P()+v1->P())/2.0;
+    Point3f delta = (v1->P()-v0->P()).Normalize();
+    Point3f n0 =  (v0->N() ^ delta).Normalize();
+    Point3f n1 =  (v1->N() ^ delta).Normalize();
+    Point3f n = (n0+n1)/2.0;
+    pl.Init(mid,n);
+    return pl;
+  }
+
+  void ComputePlaneField(MeshType &poly, EdgeGrid &edgeGrid)
+  {
+    // First Compute per-edge planes
+    std::vector<Plane3f> planeVec(poly.en);
+    for(int i=0;i<poly.en;++i)
+    {
+      planeVec[i] = ComputeEdgePlane(poly.edge[i].V(0), poly.edge[i].V(1));
+    }
+    
+    DumpPlanes(poly,planeVec);
+    edgeGrid.Set(poly.edge.begin(), poly.edge.end());    
+    const float maxDist= base.bbox.Diag()/10.0;
+    UpdateSelection<MeshType>::VertexClear(base);
+    
+    for(VertexIterator vi=base.vert.begin();vi!=base.vert.end();++vi)
+    {
+      Point3<ScalarType> p = vi->P();
+      
+      float minDist=maxDist;
+      Point3f closestP;
+      EdgeType *cep = vcg::tri::GetClosestEdgeBase(poly,edgeGrid,p,maxDist,minDist,closestP);        
+      if(cep)
+      {
+        int ind = tri::Index(poly,cep);        
+        vi->Q() = SignedDistancePlanePoint(planeVec[ind],p);
+        Point3f proj = planeVec[ind].Projection(p);
+        
+//        if(Distance(proj,closestP)>edge::Length(*cep))
+//        {
+//          vi->Q() =1;
+//          vi->SetS();
+//        }
+      }
+      else {
+        vi->Q() =1;
+        vi->SetS();
+      }
+    } 
+  }
+
+ 
+  void CutAlongPolyLineUsingField(MeshType &poly,EdgeGrid &edgeGrid)
+{
+  QualitySign qsPred(edgeGrid,poly,par);
+  QualitySignSplit qsSplit(edgeGrid,poly,par);
+  tri::UpdateTopology<MeshType>::FaceFace(base);
+  tri::RefineE(base,qsSplit,qsPred);     
+ } 
+
+
+  
+
+ 
+ 
+
+  
+  /*
+   * Make an edge mesh 1-manifold by splitting all the
+   * vertexes that have more than two incident edges
+   * 
+   * It performs the split in three steps. 
+   * First it collects and counts the vertices to be splitten. 
+   * Then it adds the vertices to the mesh and lastly it updates the poly with the newly added vertices. 
+   * 
+   * singSplitFlag allow to ubersplit each singularity in a number of vertex of the same order of its degreee. 
+   * 
+   */
+  
+  void DecomposeNonManifoldTree(MeshType &poly, bool singSplitFlag = true)
+  {
+    std::vector<int> degreeVec(poly.vn);
+    tri::UpdateTopology<MeshType>::VertexEdge(poly);
+    int neededVert=0;
+    int delta;
+    if(singSplitFlag) delta = 1;
+        else delta =2;
+      
+    for(VertexIterator vi=poly.vert.begin(); vi!=poly.vert.end();++vi)
+    {
+      std::vector<EdgeType *> starVec;
+      edge::VEStarVE(&*vi,starVec);
+      degreeVec[tri::Index(poly, *vi)] = starVec.size();
+      if(starVec.size()>2)
+        neededVert += starVec.size()-delta;
+    }
+  
+    VertexIterator firstVi = tri::Allocator<MeshType>::AddVertices(poly,neededVert);
+  
+    for(size_t i=0;i<degreeVec.size();++i)
+    {
+      if(degreeVec[i]>2)
+      {
+        std::vector<EdgeType *> edgeStarVec;
+        edge::VEStarVE(&(poly.vert[i]),edgeStarVec);
+        assert(edgeStarVec.size() == degreeVec[i]);
+        for(size_t j=delta;j<edgeStarVec.size();++j)
+        {
+          EdgeType *ep = edgeStarVec[j];
+          int ind; // index of the vertex to be changed
+          if(tri::Index(poly,ep->V(0)) == i) ind = 0;
+              else ind = 1;
+  
+          ep->V(ind) = &*firstVi;
+          ep->V(ind)->P() = poly.vert[i].P();
+          ep->V(ind)->N() = poly.vert[i].N();
+          ++firstVi;
+        }
+      }
+    }
+  }
+  
+  /*
+   * Given a manifold edgemesh it returns the boundary/terminal endpoints.
+   */
+  static void FindTerminalPoints(MeshType &poly, std::vector<VertexType *> &vec)
+  {
+    tri::UpdateTopology<MeshType>::VertexEdge(poly);
+    for(VertexIterator vi=poly.vert.begin(); vi!=poly.vert.end();++vi)
+    {
+      if(edge::VEDegree<EdgeType>(&*vi)==1)
+         vec.push_back(&*vi);
+    }
+  }
+  
+  
+  // This function will decompose the input edge mesh into a set of 
+  // connected components.
+  // the vector will contain, for each connected component, a vector with all the edge indexes. 
+  void BuildConnectedComponentVectors(MeshType &poly, std::vector< std::vector< int> > &ccVec)
+  {
+    tri::UpdateFlags<MeshType>::EdgeClearV(poly);
+    tri::UpdateTopology<MeshType>::EdgeEdge(poly);
+  
+    int visitedEdgeNum=0 ;
+    int ccCnt=0;
+    
+    EdgeIterator eIt = poly.edge.begin();  
+    while(visitedEdgeNum < poly.en)
+    {
+      ccVec.resize(ccVec.size()+1);
+      while(eIt->IsV()) ++eIt;
+      EdgeType *startE = &*eIt;
+      
+      EdgeType *curEp = &*eIt;
+      int     curEi = 0;
+//      printf("Starting Visit of connected Component %i from edge %i\n",ccCnt,tri::Index(m,*eIt));
+      while( (curEp->EEp(curEi) != startE) &&
+             (curEp->EEp(curEi) != curEp) )
+      {
+        EdgeType *nextEp = curEp->EEp(curEi);
+        int nextEi =  (curEp->EEi(curEi)+1)%2;
+        curEp = nextEp;
+        curEi = nextEi;
+      }   
+  
+      curEp->SetV();
+      curEi = (curEi +1)%2; // Flip the visit direction!
+      visitedEdgeNum++;
+      ccVec[ccCnt].push_back(tri::Index(poly,curEp));        
+      while(! curEp->EEp(curEi)->IsV())
+      {
+        EdgeType *nextEp = curEp->EEp(curEi);
+        int nextEi =  (curEp->EEi(curEi)+1)%2;
+        curEp->SetV();
+        curEp = nextEp;
+        curEi = nextEi;
+        curEp->V(curEi)->C() = Color4b::Scatter(30,ccCnt%30);
+        if(!curEp->IsV()) {
+          ccVec[ccCnt].push_back(tri::Index(poly,curEp));      
+          visitedEdgeNum++;      
+        }          
+      }
+      ccCnt++;    
+    }
+    printf("en %i - VisitedEdgeNum = %i\n",poly.en, visitedEdgeNum);
+  }
+  
+  void ExtractSubMesh(MeshType &poly, std::vector<int> &ind, MeshType &subPoly)
+  {
+    subPoly.Clear();
+    std::vector<int> remap(poly.vert.size(),-1);
+    for(size_t i=0;i<ind.size();++i)
+    {
+      int v0 = tri::Index(poly,poly.edge[ind[i]].V(0));
+      int v1 = tri::Index(poly,poly.edge[ind[i]].V(1));
+      if(remap[v0]==-1) 
+      {
+        remap[v0]=subPoly.vn;
+        tri::Allocator<MeshType>::AddVertex(subPoly,poly.vert[v0].P(),poly.vert[v0].N());      
+      }
+      if(remap[v1]==-1) 
+      {
+        remap[v1]=subPoly.vn;
+        tri::Allocator<MeshType>::AddVertex(subPoly,poly.vert[v1].P(),poly.vert[v1].N());      
+      }
+      tri::Allocator<MeshType>::AddEdge(subPoly, &subPoly.vert[remap[v0]], &subPoly.vert[remap[v1]]);   
+    }  
+  }
+  
+  
+  void Reorient(MeshType &poly, std::vector< std::vector< int> > &ccVec)
+  {
+    UpdateTopology<MeshType>::EdgeEdge(poly);
+    for(size_t i=0;i<ccVec.size();++i)
+    {
+      std::vector<bool> toFlipVec(ccVec[i].size(),false);
+      
+      for(int j=0;j<ccVec[i].size();++j)
+      {
+        EdgeType *cur  = & poly.edge[ccVec[i][j]];
+        EdgeType *prev;
+        if(j==0) prev = cur;
+            else prev = & poly.edge[ccVec[i][j-1]];
+        
+        if(cur->EEp(0) != prev)
+        {
+          toFlipVec[j] = true;
+          assert(cur->EEp(1) == prev);
+        }      
+      }
+      for(int j=0;j<ccVec[i].size();++j)
+        if(toFlipVec[j]) 
+          std::swap(poly.edge[ccVec[i][j]].V(0), poly.edge[ccVec[i][j]].V(1));    
+    }
+  }
+  
+  
+  /*
+   * Given a mesh and vector of vertex pointers it splits the mesh with the given points.
+   * To avoid degeneracies it snaps the splitting points that came nearest than a given
+   * threshold to the edge or the vertex of the closest triangle. When an input vertex
+   * is snapped the passed vertex position is modiried accordingly
+   * (this is the reason for having a vector of vertex pointers instead just a vector of points)
+   *
+   */
+  void SplitMeshWithPoints(MeshType &m, std::vector<VertexType *> &vec)
+  {
+    printf("Splitting with %i vertices\n",vec.size());
+    int faceToAdd=0;
+    int vertToAdd=0;
+  
+    // For each splitting point we save the index of the face to be splitten and the "kind" of split to do:
+    // 3 -> means classical 1 to 3 face split
+    // 2 -> means edge split.
+    // 0 -> means no need of a split (e.g. the point is coincident with a vertex)
+  
+    std::vector< std::pair<int,int> > toSplitVec(vec.size(), std::make_pair(0,0));
+    MeshGrid uniformGrid;
+    uniformGrid.Set(m.face.begin(), m.face.end());
+    const float eps=0.01f;
+    const float maxDist = m.bbox.Diag()/10.0;
+  
+    for(size_t i =0; i<vec.size();++i)
+    {
+      Point3f newP = vec[i]->P();
+      float minDist;
+      Point3f closestP,ip;
+      FaceType *f = vcg::tri::GetClosestFaceBase(m,uniformGrid,newP,maxDist, minDist, closestP);
+      assert(f);
+  
+  //    if(ip[0]<eps || ip[1]<eps || ip[2]<eps)
+  //    {
+  //      // TODO
+  //    }
+  //    else
+      {
+        toSplitVec[i].first = tri::Index(m,f);
+        toSplitVec[i].second = 3;
+        faceToAdd += 2;
+        vertToAdd += 1;
+      }
+    }
+    printf("adding %i faces and %i vertices\n",faceToAdd,vertToAdd);
+    FaceIterator newFi = tri::Allocator<MeshType>::AddFaces(m,faceToAdd);
+    VertexIterator newVi = tri::Allocator<MeshType>::AddVertices(m,vertToAdd);
+  
+    tri::UpdateColor<MeshType>::PerFaceConstant(m,Color4b::White);
+  
+    for(size_t i =0; i<vec.size();++i)
+    {
+      if(toSplitVec[i].second==3)
+      {
+       FaceType *fp0,*fp1,*fp2;
+       fp0 = &m.face[toSplitVec[i].first];
+       fp1 = &*newFi; newFi++;
+       fp2 = &*newFi; newFi++;
+       VertexType *vp = &*(newVi++);
+       vp->P() = vec[i]->P();
+       VertexType *vp0 = fp0->V(0);
+       VertexType *vp1 = fp0->V(1);
+       VertexType *vp2 = fp0->V(2);
+  
+       fp0->V(0) = vp0; fp0->V(1) = vp1; fp0->V(2) = vp;
+       fp1->V(0) = vp1; fp1->V(1) = vp2; fp1->V(2) = vp;
+       fp2->V(0) = vp2; fp2->V(1) = vp0; fp2->V(2) = vp;
+  
+       fp0->C() = Color4b::Green;
+       fp1->C() = Color4b::Green;
+       fp2->C() = Color4b::Green;
+      }
+    }
+  
+    }
+  
+  
+  void Init()
+  {
+    // Construction of the uniform grid
+    uniformGrid.Set(base.face.begin(), base.face.end());
+    UpdateNormal<MeshType>::PerFaceNormalized(base);
+    UpdateTopology<MeshType>::FaceFace(base);    
+    uniformGrid.Set(base.face.begin(), base.face.end());    
+  }
+  
+
+  
+//  Point3f GeodesicFlatten(PosType &p)
+//  {
+//    Point3f N0 = p.F()->N();
+//    Point3f N1 = p.FFlip()->N();
+//    PosType t=p;
+//    t.FlipF();
+//    t.FlipE();
+//    t.FlipV();
+//    // Now rotate the other point around the edge so that we get the two triangles on the same plane
+//    Eigen::Vector3f otherPoint,sharedEdgeDir;
+//    t.P().ToEigenVector(otherPoint);
+//    (p.V().P() - p.VFlip()->P()).Normalize().ToEigenVector(sharedEdgeDir);        
+//    ScalarType dihedralAngleRad = vcg::face::DihedralAngleRad(*p.F(),t.F());
+//    Eigen::Matrix3f mRot = Eigen::AngleAxisf(dihedralAngleRad,sharedEdgeDir).matrix();
+//    Point3f otherPointRot; 
+//    otherPointRot.FromEigenVector(mRot*otherPoint);    
+//    Plane3f facePlane; facePlane.Init(p->F()->P0(),p->F()->P1(),p->F()->P2());
+//    float dist = SignedDistancePlanePoint(facePlane,otherPointRot);    
+//  }
+
+  
+  void Simplify( MeshType &poly)
+  {
+    int startEn = poly.en;
+    Distribution<ScalarType> hist;
+    for(int i =0; i<poly.en;++i) 
+      hist.Add(edge::Length(poly.edge[i]));
+        
+    UpdateTopology<MeshType>::VertexEdge(poly);
+    
+    for(int i =0; i<poly.vn;++i)
+    {
+      std::vector<VertexPointer> starVecVp;
+      edge::VVStarVE<EdgeType>(&(poly.vert[i]),starVecVp);      
+      if(starVecVp.size()==2)
+      {      
+        float newSegLen = Distance(starVecVp[0]->P(), starVecVp[1]->P());
+        Segment3f seg(starVecVp[0]->P(),starVecVp[1]->P());
+        float segDist;
+        Point3f closestPSeg;
+        SegmentPointDistance(seg,poly.vert[i].cP(),closestPSeg,segDist);
+        Point3f fp,fn;
+        float maxSurfDist = MaxSegDist(starVecVp[0], starVecVp[1],fp,fn);
+        
+        if(maxSurfDist < par.surfDistThr && (newSegLen < par.maxSimpEdgeLen) )
+        {
+          edge::VEEdgeCollapse(poly,&(poly.vert[i]));
+        }
+      }
+    }
+    tri::Allocator<MeshType>::CompactEveryVector(poly);
+    printf("Simplify %i -> %i\n",startEn,poly.en);
+  }
+  
+  void EvaluateHausdorffDistance(MeshType &poly, Distribution<ScalarType> &dist)
+  {
+    dist.Clear();
+    tri::UpdateTopology<MeshType>::VertexEdge(poly);
+    tri::UpdateQuality<MeshType>::VertexConstant(poly,0);
+    for(int i =0; i<poly.edge.size();++i)
+    {      
+      Point3f farthestP, farthestN;      
+      float maxDist = MaxSegDist(poly.edge[i].V(0),poly.edge[i].V(1), farthestP, farthestN, &dist);      
+      poly.edge[i].V(0)->Q()+= maxDist;
+      poly.edge[i].V(1)->Q()+= maxDist;
+    }
+    for(int i=0;i<poly.vn;++i)
+    {
+      ScalarType deg = edge::VEDegree<EdgeType>(&poly.vert[i]);
+      poly.vert[i].Q()/=deg;
+    }
+    tri::UpdateColor<MeshType>::PerVertexQualityRamp(poly,0,dist.Max());    
+  }
+  // Given a segment find the maximum distance from it to the original surface. 
+  float MaxSegDist(VertexType *v0, VertexType *v1, Point3f &farthestPointOnSurf, Point3f &farthestN, Distribution<ScalarType> *dist=0)
+  {
+    float maxSurfDist = 0;
+    const float sampleNum = 10;
+    const float maxDist = base.bbox.Diag()/10.0;
+    for(float k = 1;k<sampleNum;++k)
+    {
+      float surfDist;
+      Point3f closestPSurf;
+      Point3f samplePnt = (v0->P()*k +v1->P()*(sampleNum-k))/sampleNum;          
+      FaceType *f = vcg::tri::GetClosestFaceBase(base,uniformGrid,samplePnt,maxDist, surfDist, closestPSurf);        
+      if(dist)
+        dist->Add(surfDist);
+      assert(f);
+      if(surfDist > maxSurfDist)
+      {
+        maxSurfDist = surfDist;
+        farthestPointOnSurf = closestPSurf;
+        farthestN = f->N();
+      }
+    }
+    return maxSurfDist;
+  }
+  
+  void Refine( MeshType &poly)
+  {
+    tri::Allocator<MeshType>::CompactEveryVector(poly);    
+    int startEdgeSize = poly.en;
+    for(int i =0; i<startEdgeSize;++i)
+    {
+      if(edge::Length(poly.edge[i])>par.minRefEdgeLen)  
+      {      
+        Point3f farthestP, farthestN;
+        float maxDist = MaxSegDist(poly.edge[i].V(0),poly.edge[i].V(1),farthestP, farthestN);
+        if(maxDist > par.surfDistThr)  
+        {
+          VertexIterator vi = Allocator<MeshType>::AddVertex(poly,farthestP, farthestN);        
+          Allocator<MeshType>::AddEdge(poly,poly.edge[i].V(0),&*vi);
+          Allocator<MeshType>::AddEdge(poly,&*vi,poly.edge[i].V(1));
+          Allocator<MeshType>::DeleteEdge(poly,poly.edge[i]);
+        }
+      }
+    }
+    tri::Allocator<MeshType>::CompactEveryVector(poly);
+    printf("Refine %i -> %i\n",startEdgeSize,poly.en);fflush(stdout);
+  }
+  
+  void SmoothProject(MeshType &poly, int iterNum, ScalarType smoothWeight, ScalarType projectWeight)
+  {
+    assert(poly.en>0 && base.fn>0);
+    for(int k=0;k<iterNum;++k)
+    {
+      std::vector<Point3f> posVec(poly.vn,Point3f(0,0,0));
+      std::vector<int>     cntVec(poly.vn,0);
+  
+      for(int i =0; i<poly.en;++i)
+      {
+        for(int j=0;j<2;++j)
+        {
+          int vertInd = tri::Index(poly,poly.edge[i].V(j));
+          posVec[vertInd] += poly.edge[i].V1(j)->P();
+          cntVec[vertInd] += 1;
+        }
+      }
+      
+      const float maxDist = base.bbox.Diag()/10.0;
+      for(int i=0; i<poly.vn; ++i)
+      {
+        Point3f smoothPos = (poly.vert[i].P() + posVec[i])/float(cntVec[i]+1);
+        
+        Point3f newP = poly.vert[i].P()*(1.0-smoothWeight) + smoothPos *smoothWeight;
+        
+        Point3f delta =  newP - poly.vert[i].P();
+        if(delta.Norm() > par.maxSmoothDelta) 
+        {
+          newP =  poly.vert[i].P() + ( delta / delta.Norm()) * maxDist*0.5;
+        }
+        
+        float minDist;
+        Point3f closestP;
+        FaceType *f = vcg::tri::GetClosestFaceBase(base,uniformGrid,newP,maxDist, minDist, closestP);
+        assert(f);
+        poly.vert[i].P() = newP*(1.0-projectWeight) +closestP*projectWeight;
+        poly.vert[i].N() = f->N();
+      }
+      
+      Refine(poly);      
+      Refine(poly);      
+      Simplify(poly);            
+    }
+  }  
+   
+  
+};
+} // end namespace tri
+} // end namespace vcg
+
+#endif // __VCGLIB_CURVE_ON_SURF_H