diff --git a/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.cpp b/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.cpp
new file mode 100644
index 00000000..523e6389
--- /dev/null
+++ b/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.cpp
@@ -0,0 +1,77 @@
+/****************************************************************************
+* VCGLib                                                            o o     *
+* Visual and Computer Graphics Library                            o     o   *
+*                                                                _   O  _   *
+* Copyright(C) 2004-2009                                           \/)\/    *
+* 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.                                                         *
+*                                                                           *
+****************************************************************************/
+#include<vcg/complex/complex.h>
+#include<wrap/io_trimesh/import.h>
+#include<wrap/io_trimesh/export_ply.h>
+#include <vcg/complex/algorithms/voronoi_remesher.h>
+
+using namespace vcg;
+using namespace std;
+
+class MyVertex;
+class MyFace;
+class MyEdge;
+
+struct MyUsedTypes : public vcg::UsedTypes<
+        vcg::Use<MyVertex>::AsVertexType,
+        vcg::Use<MyFace>::AsFaceType,
+        vcg::Use<MyEdge>::AsEdgeType> {};
+
+class MyVertex: public vcg::Vertex<MyUsedTypes,
+         vcg::vertex::Coord3f,  vcg::vertex::Normal3f,
+         vcg::vertex::Color4b,  vcg::vertex::Qualityd,
+         vcg::vertex::VFAdj,    vcg::vertex::VEAdj,
+         vcg::vertex::BitFlags, vcg::vertex::Mark> {};
+class MyFace   : public vcg::Face<MyUsedTypes,
+         vcg::face::VertexRef, vcg::face::Normal3f,
+         vcg::face::Color4b, vcg::face::BitFlags,
+         vcg::face::VFAdj, vcg::face::FFAdj,
+         vcg::face::Mark> {};
+class MyEdge: public vcg::Edge<MyUsedTypes,
+         vcg::edge::VertexRef, vcg::edge::BitFlags,
+         vcg::edge::EEAdj,    vcg::edge::VEAdj> {};
+
+class MyMesh   : public vcg::tri::TriMesh< std::vector<MyVertex>, std::vector<MyFace> > {};
+
+int main( int argc, char **argv )
+{
+  MyMesh startMesh;
+  if(argc < 2 )
+  {
+    printf("Usage trimesh_voro mesh region_num\n");
+     return -1;
+  }
+  printf("Reading %s  \n",argv[1]);
+  int ret= tri::io::Importer<MyMesh>::Open(startMesh,argv[1]);
+  if(ret!=0)
+  {
+    printf("Unable to open %s for '%s'\n",argv[1],tri::io::ImporterPLY<MyMesh>::ErrorMsg(ret));
+    return -1;
+  }
+  float samplingRadius = 0.3;
+  auto remeshed = Remesher<MyMesh>::Remesh(startMesh, samplingRadius, 50.0);
+  
+  
+  tri::io::ExporterPLY<MyMesh>::Save(*remeshed,"Full.ply",tri::io::Mask::IOM_VERTCOLOR|tri::io::Mask::IOM_WEDGTEXCOORD );
+  return 0;
+}
diff --git a/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.pro b/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.pro
new file mode 100644
index 00000000..cba125e8
--- /dev/null
+++ b/apps/sample/trimesh_voronoi_remesh/trimesh_voronoi_remesh.pro
@@ -0,0 +1,3 @@
+include(../common.pri)
+TARGET = trimesh_voronoi_remesh
+SOURCES += trimesh_voronoi_remesh.cpp ../../../wrap/ply/plylib.cpp
diff --git a/vcg/complex/algorithms/voronoi_remesher.h b/vcg/complex/algorithms/voronoi_remesher.h
new file mode 100644
index 00000000..f586b2a9
--- /dev/null
+++ b/vcg/complex/algorithms/voronoi_remesher.h
@@ -0,0 +1,705 @@
+/****************************************************************************
+* VCGLib                                                            o o     *
+* Visual and Computer Graphics Library                            o     o   *
+*                                                                _   O  _   *
+* Copyright(C) 2004-2016                                           \/)\/    *
+* 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_VORONOI_REMESHER_H
+#define _VCGLIB_VORONOI_REMESHER_H
+
+#include <vcg/complex/complex.h>
+#include <vcg/complex/algorithms/update/topology.h>
+#include <vcg/complex/algorithms/refine.h>
+#include <vcg/complex/algorithms/clean.h>
+#include <vcg/simplex/face/pos.h>
+#include <vcg/complex/algorithms/voronoi_processing.h>
+#include <vcg/complex/algorithms/point_sampling.h>
+#include <vcg/complex/algorithms/crease_cut.h>
+
+#include <memory>
+#include <string>
+#include <vector>
+#include <map>
+#include <unordered_map>
+#include <unordered_set>
+#include <cmath>
+#include <array>
+#include <utility>
+
+//#define DEBUG_VORO 1
+namespace vcg {
+namespace tri {
+
+template <class MeshType, class EdgeMeshType = MeshType>
+class Remesher
+{
+public:
+  typedef Remesher                  ThisType;
+  
+  typedef MeshType                      Mesh;
+  typedef typename Mesh::ScalarType     ScalarType;
+  typedef typename Mesh::CoordType      CoordType;
+  typedef typename Mesh::FaceType       FaceType;
+  typedef typename Mesh::FacePointer    FacePointer;
+  typedef typename Mesh::VertexType     VertexType;
+  typedef typename Mesh::VertexPointer  VertexPointer;
+  typedef typename Mesh::FaceIterator   FaceIterator;
+  typedef typename Mesh::VertexIterator VertexIterator;
+  
+  typedef std::shared_ptr<Mesh>         MeshPtr;
+  
+protected:
+  typedef face::Pos<FaceType>                            PosType;
+  
+  
+  /// \brief splitCC split the provided mesh into connected components.
+  /// \param mesh the inputMesh.
+  /// \return the vector of connected components (meshes) for the input model
+  /// (if the input mesh is a single connected component returns an empty vector).
+  ///
+  inline static std::vector<MeshPtr> splitCC(MeshType & mesh)
+  {
+    std::vector<MeshPtr> ret;
+    
+    // find the connected components
+    std::vector<std::pair<int, typename MeshType::FacePointer> > CCV;
+    Clean<MeshType>::ConnectedComponents(mesh, CCV);
+    
+    if (CCV.size() == 1)
+      return ret;
+    
+    ConnectedComponentIterator<MeshType> ci;
+    for(size_t i=0; i<CCV.size(); ++i)
+    {
+      // clear selection
+      UpdateSelection<MeshType>::Clear(mesh);
+      for(ci.start(mesh, CCV[i].second); !ci.completed(); ++ci)
+      {
+        // select all faces for a CC
+        (*ci)->SetS();
+      }
+      
+      // create from selected
+      MeshPtr cc = std::make_shared<MeshType>();
+      Append<MeshType, MeshType>::MeshCopy(*cc, mesh, true);
+      ret.push_back(cc);
+    }
+    
+    return ret;
+  }
+  
+public:
+  static const int VoroRelaxationStep = 20;
+  
+  ///
+  /// \brief Remesh the main function that remeshes a mesh preserving creases.
+  /// \param original the mesh
+  /// \param samplingRadius is the sampling ragius for remeshing
+  /// \param borderCreaseAngleThreshold is the angle treshold for preserving corner points on the mesh boundary
+  /// \return the remeshed mesh
+  ///
+  static inline MeshPtr Remesh(Mesh & original, const ScalarType samplingRadius, const ScalarType borderCreaseAngleThreshold = 0.0)
+  {
+    UpdateTopology<Mesh>::FaceFace(original);
+    UpdateFlags<Mesh>::FaceBorderFromFF(original);
+    UpdateFlags<Mesh>::VertexBorderFromFaceAdj(original);
+    
+    if (Clean<Mesh>::CountNonManifoldEdgeFF(original) > 0)
+    {
+      std::cout << "Input mesh has non manifold edges" << std::endl;
+      return nullptr;
+    }
+    
+    // for closed watertight mesh try to split
+    if (Clean<Mesh>::CountHoles(original) < 1)
+    {
+      CreaseCut<Mesh>(original, vcg::math::ToRad(70.0));
+      Allocator<Mesh>::CompactEveryVector(original);
+      UpdateTopology<Mesh>::FaceFace(original);
+      UpdateFlags<Mesh>::FaceBorderFromFF(original);
+      UpdateFlags<Mesh>::VertexBorderFromFaceAdj(original);
+#ifdef DEBUG_VORO
+      io::Exporter<Mesh>::Save(original, "creaseSplit.ply", 0);
+#endif
+    }
+    
+    // One CC
+    std::vector<MeshPtr> ccs = splitCC(original);
+    if (ccs.empty())
+      return RemeshOneCC(original, samplingRadius, borderCreaseAngleThreshold);
+    
+    
+    // Multiple CCs
+    std::cout << "Remeshing " << ccs.size() << " components" << std::endl;
+    for (size_t i=0; i<ccs.size(); i++)
+    {
+      std::cout << "Remeshing component " << (i+1) << "/" << ccs.size() << std::endl;
+      ccs[i] = RemeshOneCC(*ccs[i], samplingRadius, borderCreaseAngleThreshold);
+    }
+    
+    MeshPtr ret = std::make_shared<Mesh>();
+    for (MeshPtr & mesh : ccs)
+    {
+      Append<Mesh,Mesh>::Mesh(*ret, *mesh);
+    }
+    Clean<Mesh>::RemoveDuplicateVertex(*ret, true);
+    return ret;
+  }
+  
+  ///
+  /// \brief RemeshOneCC the function that remeshes a single connected component mesh preserving its boundary (consistently for eventually adjacent meshes).
+  /// \param original the mesh
+  /// \param samplingRadius is the sampling ragius for remeshing
+  /// \param borderCreaseAngleThreshold is the angle treshold for preserving corner points on the mesh boundary
+  /// \return the remeshed mesh
+  ///
+  static inline MeshPtr RemeshOneCC(Mesh & original, const ScalarType samplingRadius, const ScalarType borderCreaseAngleThreshold = 0.0)
+  {
+    RequireCompactness(original);
+    RequirePerFaceFlags(original);
+    
+    UpdateTopology<Mesh>::FaceFace(original);
+    UpdateFlags<Mesh>::FaceBorderFromFF(original);
+    UpdateFlags<Mesh>::VertexBorderFromFaceAdj(original);
+    
+    
+    // Resample border
+    Mesh poissonEdgeMesh;
+    {
+      typedef typename EdgeMeshType::CoordType Coord;
+      
+      EdgeMeshType em;
+      //			ThisType::ExtractMeshSides(original, em);
+      
+      ThisType::ExtractMeshBorders(original, em);
+      
+      // wtf we should close the loops
+      Clean<EdgeMeshType>::RemoveDuplicateVertex(em);
+      Allocator<EdgeMeshType>::CompactVertexVector(em);
+      Allocator<EdgeMeshType>::CompactEdgeVector(em);
+      
+#ifdef DEBUG_VORO
+      io::ExporterOBJ<EdgeMeshType>::Save(em, QString("edgeMesh_%1.obj").arg(idx).toStdString().c_str(), io::Mask::IOM_EDGEINDEX);
+#endif
+      
+      // eventually split on 'creases'
+      if (borderCreaseAngleThreshold > 0.0)
+      {
+        UpdateFlags<EdgeMeshType>::VertexClearS(em);
+        UpdateFlags<EdgeMeshType>::VertexClearV(em);
+        Clean<EdgeMeshType>::SelectCreaseVertexOnEdgeMesh(em, vcg::math::ToRad(borderCreaseAngleThreshold));
+        std::cout << Clean<EdgeMeshType>::SplitSelectedVertexOnEdgeMesh(em) << " splits" << std::endl;
+      }
+#ifdef DEBUG_VORO
+      io::ExporterOBJ<EdgeMeshType>::Save(em, QString("edgeMesh_split_%1.obj").arg(idx).toStdString().c_str(), io::Mask::IOM_EDGEINDEX);
+#endif
+      
+      // Samples vector
+      std::vector<Coord> borderSamples;
+      TrivialSampler<EdgeMeshType> ps(borderSamples);
+      
+      // uniform sampling
+      // (use different sampling radius for the edges)
+      UpdateTopology<EdgeMeshType>::EdgeEdge(em);
+      SurfaceSampling<EdgeMeshType>::EdgeMeshUniform(em, ps, samplingRadius, true);
+      
+      // convert to mesh
+      auto vi = Allocator<Mesh>::AddVertices(poissonEdgeMesh, borderSamples.size());
+      for (auto p : borderSamples)
+      {
+        vi->P() = CoordType::Construct(p);
+        vi++;
+      }
+      UpdateBounding<Mesh>::Box(poissonEdgeMesh);
+      
+      // remove duplicate vertices
+      Clean<Mesh>::RemoveDuplicateVertex(poissonEdgeMesh, false);
+      Allocator<Mesh>::CompactVertexVector(poissonEdgeMesh);
+      
+      // select all vertices (to mark them fixed)
+      UpdateFlags<Mesh>::VertexSetS(poissonEdgeMesh);
+      
+#ifdef DEBUG_VORO
+      //			// temp remove
+      //			UpdateColor<Mesh>::PerVertexConstant(poissonEdgeMesh, vcg::Color4b::Gray);
+      
+      //			typedef typename vcg::SpatialHashTable<VertexType, ScalarType> HashVertexGrid;
+      //			HashVertexGrid HG;
+      //			HG.Set(poissonEdgeMesh.vert.begin(),poissonEdgeMesh.vert.end());
+      //			for (size_t i=0; i<creases.size(); i++)
+      //			{
+      //				const float dist_upper_bound=FLT_MAX;
+      //				ScalarType dist;
+      //				VertexType * vp = GetClosestVertex<MeshType,HashVertexGrid>(poissonEdgeMesh, HG, creases[i], dist_upper_bound, dist);
+      //				assert(vp);
+      //				vp->C() = vcg::Color4b::Red;
+      //			}
+      io::ExporterPLY<MeshType>::Save(poissonEdgeMesh, QString("borderMesh_%1.ply").arg(idx).toStdString().c_str(), io::Mask::IOM_VERTCOLOR);
+#endif
+    }
+    
+    typedef VoronoiProcessing<Mesh>            Voronoi;
+    typedef TrivialSampler<Mesh>               BaseSampler;
+    typedef SurfaceSampling<Mesh, BaseSampler> SurfaceSampler;
+    typedef SurfaceSampling<Mesh, FixSampler>  SurfaceFixSampler;
+    
+    // copy original mesh
+    Mesh baseMesh;
+    Append<Mesh, Mesh>::MeshCopy(baseMesh, original, false, true);
+    
+    // refine to obtain a base mesh
+    VoronoiProcessingParameter vpp;
+    vpp.refinementRatio = 4.0f;
+    Voronoi::PreprocessForVoronoi(baseMesh, samplingRadius, vpp);
+    
+    // Poisson sampling preserving border
+    Mesh poissonMesh;
+    std::vector<CoordType> seedPointVec;
+    std::vector<bool>      seedFixedVec;
+    FixSampler fix_sampler(seedPointVec, seedFixedVec);
+    
+    // montecarlo sampler
+    std::vector<CoordType> sampleVec;
+    BaseSampler mps(sampleVec);
+    
+    // NOTE in order to make results always the same the random sampling generator is reinitialized for
+    // for each patch resampling
+    SurfaceSampler::SamplingRandomGenerator().initialize(5489u);
+    
+    // Montecarlo oversampling
+    Mesh montecarloMesh;
+    int poissonCount = SurfaceSampler::ComputePoissonSampleNum(original, samplingRadius) * 0.7;
+    
+    std::cout << "poisson Count: " << poissonCount << std::endl;
+    if (poissonCount <= 0)
+    {
+      // no need for internal sampling
+      Append<Mesh, Mesh>::MeshCopy(poissonMesh, poissonEdgeMesh);
+      for (auto vi = poissonEdgeMesh.vert.begin(); vi != poissonEdgeMesh.vert.end(); vi++)
+      {
+        fix_sampler.AddVert(*vi);
+      }
+    }
+    else
+    {
+      // Montecarlo poisson sampling
+      SurfaceSampler::MontecarloPoisson(original, mps, poissonCount * 20);
+      BuildMeshFromCoordVector(montecarloMesh,sampleVec);
+      
+#ifdef DEBUG_VORO
+      io::ExporterPLY<MeshType>::Save(montecarloMesh, QString("montecarloMesh_%1.ply").arg(idx).toStdString().c_str());
+#endif
+      
+      // Poisson disk pruning initialized with edges
+      typename SurfaceFixSampler::PoissonDiskParam pp;
+      pp.preGenMesh = &poissonEdgeMesh;
+      pp.preGenFlag = true;
+      SurfaceFixSampler::PoissonDiskPruning(fix_sampler, montecarloMesh, samplingRadius, pp);
+      
+#ifdef DEBUG_VORO
+      BuildMeshFromCoordVector(poissonMesh,seedPointVec);
+      io::ExporterPLY<MeshType>::Save(poissonMesh, QString("poissonMesh_%1.ply").arg(idx).toStdString().c_str());
+#endif
+    }
+    
+    
+    std::cout << "poisson samples " << seedPointVec.size() << std::endl;
+    
+    // restricted relaxation with fixed points
+    vpp.seedPerturbationProbability = 0.0f;
+    // TODO check preserveFixedSeed flag (NO)
+    Voronoi::RestrictedVoronoiRelaxing(baseMesh, seedPointVec, seedFixedVec, VoroRelaxationStep, vpp);
+    
+#ifdef DEBUG_VORO
+    BuildMeshFromCoordVector(poissonMesh,seedPointVec);
+    io::ExporterPLY<MeshType>::Save(poissonMesh, QString("relaxedMesh_%1.ply").arg(idx).toStdString().c_str());
+#endif
+    
+    // FAIL?
+    MeshPtr finalMeshPtr = std::make_shared<Mesh>();
+    std::vector<VertexType *> seedVertexVec;
+    //		Voronoi::SeedToVertexConversion(baseMesh, seedPointVec, seedVertexVec, false);
+    ThisType::SeedToFixedBorderVertexConversion(baseMesh, seedPointVec, seedFixedVec, seedVertexVec);
+    EuclideanDistance<Mesh> dd;
+    std::cout << "BEGIN compute vertex sources (basemesh vn:" << baseMesh.VN() << " fn:" << baseMesh.FN() << ")" << std::endl;
+    
+    Voronoi::ComputePerVertexSources(baseMesh, seedVertexVec, dd);
+    std::cout << "END   compute vertex sources" << std::endl;
+    //		Voronoi::ConvertDelaunayTriangulationToMesh(baseMesh, *finalMeshPtr, seedVertexVec, false); // traditional
+    ThisType::ConvertDelaunayTriangulationExtendedToMesh(baseMesh, *finalMeshPtr, seedVertexVec); // border-preserving
+    
+#ifdef DEBUG_VORO
+    io::ExporterPLY<MeshType>::Save(*finalMeshPtr, QString("voroMesh_%1.ply").arg(idx).toStdString().c_str());
+    io::ExporterPLY<MeshType>::Save(baseMesh, QString("baseMesh_%1.ply").arg(idx).toStdString().c_str(), io::Mask::IOM_VERTCOLOR);
+#endif
+    
+    return finalMeshPtr;
+  }
+  
+protected:
+  static inline void ExtractMeshBorders(Mesh & mesh, EdgeMeshType & sides)
+  {
+    RequireFFAdjacency(mesh);
+    RequireVFAdjacency(mesh);
+    
+    // clean the edge mesh containing the borders
+    sides.Clear();
+    
+    // gather into separate vertices lists
+    std::vector<std::vector<VertexType *> > edges;
+    
+    for (auto fi = mesh.face.begin(); fi != mesh.face.end(); fi++)
+    {
+      for (int e=0; e<fi->VN(); e++)
+      {
+        if (vcg::face::IsBorder(*fi, e))
+        {
+          std::vector<VertexType *> tmp;
+          tmp.push_back(fi->V(e));
+          tmp.push_back(fi->V((e+1)%fi->VN()));
+          edges.push_back(tmp);
+        }
+      }
+    }
+    
+    // convert to edge mesh
+    for (auto & e : edges)
+    {
+      assert(e.size() >= 2);
+      
+      std::vector<typename EdgeMeshType::VertexType *> newVtx;
+      
+      // insert new vertices and store their pointer
+      auto vi = Allocator<EdgeMeshType>::AddVertices(sides, e.size());
+      for (const auto & v : e)
+      {
+        vi->ImportData(*v);
+        newVtx.push_back(&(*vi++));
+      }
+      
+      auto ei = Allocator<EdgeMeshType>::AddEdges(sides, e.size() - 1);
+      for (int i=0; i<static_cast<int>(e.size() - 1); i++)
+      {
+        ei->V(0) = newVtx[i];
+        ei->V(1) = newVtx[i+1];
+        ei++;
+      }
+    }
+    
+    Clean<EdgeMeshType>::RemoveDuplicateVertex(sides);
+  }
+  
+  
+  ///
+  /// \brief ExtractMeshSides
+  /// \param mesh the mesh (topology already computed)
+  /// \param sides the edge mesh filled with the extracted borders
+  ///
+  static inline void ExtractMeshSides(Mesh & mesh, EdgeMeshType & sides)
+  {
+    // TODO change this.... maybe wrong
+    RequireFFAdjacency(mesh);
+    RequireVFAdjacency(mesh);
+    
+    // clean the edge mesh containing the borders
+    sides.Clear();
+    
+    // find a border edge
+    assert(Clean<Mesh>::CountHoles(mesh) >= 1);
+    PosType pos;
+    for (auto fi = mesh.face.begin(); fi != mesh.face.end() && pos.IsNull(); fi++)
+    {
+      for (int e=0; e<fi->VN(); e++)
+      {
+        if (vcg::face::IsBorder(*fi, e))
+        {
+          pos = PosType(&(*fi), e);
+          break;
+        }
+      }
+    }
+    assert(!pos.IsNull());
+    assert(pos.IsBorder());
+    
+    // navigate to a corner
+    VertexType * v = pos.V();
+    do
+    {
+      pos.NextB();
+    } while(!pos.V()->IsV() && pos.V() != v);
+    // if it's a loop mark the initial point as a corner
+    pos.V()->SetV();
+    v = pos.V();
+    
+    // gather into separate vertices lists
+    std::vector<std::vector<VertexType *> > edges;
+    std::vector<VertexType *> edgePtrVec;
+    do
+    {
+      edgePtrVec.push_back(pos.V());
+      pos.NextB();
+      if (pos.V()->IsV())
+      {
+        edgePtrVec.push_back(pos.V());
+        edges.push_back(edgePtrVec);
+        edgePtrVec.clear();
+      }
+    } while (pos.V() != v);
+    
+    // convert to edge mesh
+    for (auto & e : edges)
+    {
+      assert(e.size() >= 2);
+      
+      std::vector<typename EdgeMeshType::VertexType *> newVtx;
+      
+      // insert new vertices and store their pointer
+      auto vi = Allocator<EdgeMeshType>::AddVertices(sides, e.size());
+      for (const auto & v : e)
+      {
+        vi->ImportData(*v);
+        newVtx.push_back(&(*vi++));
+      }
+      
+      auto ei = Allocator<EdgeMeshType>::AddEdges(sides, e.size() - 1);
+      for (int i=0; i<static_cast<int>(e.size() - 1); i++)
+      {
+        ei->V(0) = newVtx[i];
+        ei->V(1) = newVtx[i+1];
+        ei++;
+      }
+    }
+  }
+  
+  static void SeedToFixedBorderVertexConversion(MeshType & m,
+                                                const std::vector<CoordType> & seedPVec,
+                                                const std::vector<bool> & seedFixed,
+                                                std::vector<VertexType *> & seedVVec)
+  {
+    typedef typename vcg::SpatialHashTable<VertexType, ScalarType> HashVertexGrid;
+    seedVVec.clear();
+    
+    UpdateTopology<MeshType>::FaceFace(m);
+    UpdateFlags<MeshType>::VertexBorderFromFaceAdj(m);
+    
+    typename MeshType::BoxType bbox = m.bbox;
+    bbox.Offset(bbox.Diag()/4.0);
+    
+    // internal vertices grid
+    HashVertexGrid HG;
+    HG.Set(m.vert.begin(),m.vert.end(), bbox);
+    
+    // boundary vertices grid
+    MeshType borderMesh;
+    HashVertexGrid borderHG;
+    {
+      // get border vertices and build another mesh
+      std::vector<CoordType> borderPts;
+      for (auto vit=m.vert.begin(); vit!=m.vert.end(); vit++)
+      {
+        if (!vit->IsD() && vit->IsB())
+          borderPts.push_back(vit->cP());
+      }
+      BuildMeshFromCoordVector(borderMesh,borderPts);
+      borderMesh.bbox = m.bbox;
+      borderHG.Set(borderMesh.vert.begin(), borderMesh.vert.end(), bbox);
+    }
+    
+    const float dist_upper_bound=m.bbox.Diag()/4.0;
+    VertexType * vp = NULL;
+    
+    for( size_t i = 0; i < seedPVec.size(); i++)
+    {
+      const CoordType & p = seedPVec[i];
+      const bool fixed    = seedFixed[i];
+      
+      if (!fixed)
+      {
+        ScalarType dist;
+        vp = GetClosestVertex<MeshType,HashVertexGrid>(m, HG, p, dist_upper_bound, dist);
+      }
+      else
+      {
+        vp = NULL;
+        
+        ScalarType dist;
+        VertexType * borderVp = GetClosestVertex<MeshType,HashVertexGrid>(borderMesh, borderHG, p, dist_upper_bound, dist);
+        
+        if (borderVp)
+        {
+          vp = GetClosestVertex<MeshType,HashVertexGrid>(m, HG, borderVp->cP(), dist_upper_bound, dist);
+        }
+      }
+      
+      if (vp)
+      {
+        seedVVec.push_back(vp);
+      }
+    }
+  }
+  
+  static void ConvertDelaunayTriangulationExtendedToMesh(MeshType &m,
+                                                         MeshType &outMesh,
+                                                         std::vector<VertexType *> &seedVec)
+  {
+    typedef VoronoiProcessing<MeshType> Voronoi;
+    
+    RequirePerVertexAttribute(m ,"sources");
+    RequireCompactness(m);
+    RequireVFAdjacency(m);
+    
+    auto sources = Allocator<MeshType>::template GetPerVertexAttribute<VertexPointer> (m,"sources");
+    
+    outMesh.Clear();
+    UpdateTopology<MeshType>::FaceFace(m);
+    UpdateFlags<MeshType>::FaceBorderFromFF(m);
+    
+    std::map<VertexPointer, int> seedMap;  // It says if a given vertex of m is a seed (and its index in seedVec)
+    Voronoi::BuildSeedMap(m, seedVec, seedMap);
+    
+    std::vector<FacePointer> innerCornerVec,   // Faces adjacent to three different regions
+        borderCornerVec;  // Faces that are on the border and adjacent to at least two regions.
+    Voronoi::GetFaceCornerVec(m, sources, innerCornerVec, borderCornerVec);
+    
+    // First add all the needed vertices: seeds and corners
+    
+    for(size_t i=0;i<seedVec.size();++i)
+    {
+      Allocator<MeshType>::AddVertex(outMesh, seedVec[i]->P(), vcg::Color4b::White);
+    }
+    
+    // Now just add one face for each inner corner
+    for(size_t i=0; i<innerCornerVec.size(); ++i)
+    {
+      VertexPointer s0 = sources[innerCornerVec[i]->V(0)];
+      VertexPointer s1 = sources[innerCornerVec[i]->V(1)];
+      VertexPointer s2 = sources[innerCornerVec[i]->V(2)];
+      assert ( (s0!=s1) && (s0!=s2) && (s1!=s2) );
+      VertexPointer v0 = & outMesh.vert[seedMap[s0]];
+      VertexPointer v1 = & outMesh.vert[seedMap[s1]];
+      VertexPointer v2 = & outMesh.vert[seedMap[s2]];
+      Allocator<MeshType>::AddFace(outMesh,  v0, v1, v2 );
+    }
+    
+    // Now loop around the borders and find the missing delaunay triangles
+    // select border seed vertices only and pick one
+    UpdateFlags<Mesh>::VertexBorderFromFaceAdj(m);
+    UpdateFlags<Mesh>::VertexClearS(m);
+    UpdateFlags<Mesh>::VertexClearV(m);
+    
+    std::vector<VertexPointer> borderSeeds;
+    for (auto & s : seedVec)
+    {
+      if (s->IsB())
+      {
+        s->SetS();
+        borderSeeds.emplace_back(s);
+      }
+    }
+    
+    for (VertexPointer startBorderVertex : borderSeeds)
+    {
+      if (startBorderVertex->IsV())
+      {
+        continue;
+      }
+      
+      // unvisited border seed found
+      
+      // put the pos on the border
+      PosType pos(startBorderVertex->VFp(), startBorderVertex->VFi());
+      do {
+        pos.NextE();
+      } while (!pos.IsBorder() || (pos.VInd() != pos.E()));
+      
+      // check all border edges between each consecutive border seeds pair
+      do {
+        std::vector<VertexPointer> edgeVoroVertices(1, sources[pos.V()]);
+        //	among all sources found
+        do {
+          pos.NextB();
+          VertexPointer source = sources[pos.V()];
+          if (edgeVoroVertices.empty() || edgeVoroVertices.back() != source)
+          {
+            edgeVoroVertices.push_back(source);
+          }
+        } while (!pos.V()->IsS());
+        
+        pos.V()->SetV();
+        
+        //				assert(edgeVoroVertices.size() >= 2);
+        
+        // add face if 3 different voronoi regions are crossed by the edge
+        if (edgeVoroVertices.size() == 3)
+        {
+          VertexPointer v0 = & outMesh.vert[seedMap[edgeVoroVertices[0]]];
+          VertexPointer v1 = & outMesh.vert[seedMap[edgeVoroVertices[1]]];
+          VertexPointer v2 = & outMesh.vert[seedMap[edgeVoroVertices[2]]];
+          Allocator<MeshType>::AddFace(outMesh, v0,v1,v2);
+        }
+        
+      } while ((pos.V() != startBorderVertex));
+    }
+    
+    
+    Clean<MeshType>::RemoveUnreferencedVertex(outMesh);
+    Allocator<MeshType>::CompactVertexVector(outMesh);
+  }
+  
+  ///
+  /// \brief The FixSampler class is used with poisson disk pruning to preserve selected vertices and
+  /// keep an auxiliary vector indicating wether the sample is fixed or not
+  ///
+  class FixSampler
+  {
+  public:
+    typedef typename MeshType::CoordType  CoordType;
+    typedef typename MeshType::VertexType VertexType;
+    
+    FixSampler(std::vector<CoordType> & samples,
+               std::vector<bool>      & fixed)
+      : sampleVec(samples)
+      , fixedVec (fixed)
+    {
+      reset();
+    }
+    
+    void reset()
+    {
+      sampleVec.clear();
+      fixedVec .clear();
+    }
+    
+    void AddVert(const VertexType &p)
+    {
+      sampleVec.push_back(p.cP());
+      fixedVec .push_back(p.IsS());
+    }
+    
+  private:
+    std::vector<CoordType> & sampleVec;
+    std::vector<bool>      & fixedVec;
+  };
+};
+
+} // end namespace tri
+} // end namespace vcg
+
+#endif // REMESHER_H