diff --git a/include/Algo/Import/import.h b/include/Algo/Import/import.h
index 628a435f8fff4983bfc111e12a667d33e841b16f..000b421f115885295d3487ef5178437ffd531eb2 100644
--- a/include/Algo/Import/import.h
+++ b/include/Algo/Import/import.h
@@ -51,7 +51,7 @@ namespace Import
 * @return a boolean indicating if import was successfull
 */
 template <typename PFP>
-bool importMesh(typename PFP::MAP& map, const std::string& filename, std::vector<std::string>& attrNames, ImportSurfacique::ImportType kind = ImportSurfacique::UNKNOWNSURFACE);
+bool importMesh(typename PFP::MAP& map, const std::string& filename, std::vector<std::string>& attrNames, ImportSurfacique::ImportType kind = ImportSurfacique::UNKNOWNSURFACE, bool mergeCloseVertices=false);
 
 /**
  * import a volumic mesh
diff --git a/include/Algo/Import/import2tables.h b/include/Algo/Import/import2tables.h
index 299e6dd5a42dc364565ece399c513697aaafed75..8a6d00cb3a02b65d5f2e2488a1cbb79aefb06d13 100644
--- a/include/Algo/Import/import2tables.h
+++ b/include/Algo/Import/import2tables.h
@@ -115,6 +115,8 @@ public:
 	bool importCTM(const std::string& filename, std::vector<std::string>& attrNames);
 
 	bool importASSIMP(const std::string& filename, std::vector<std::string>& attrNames);
+	
+	bool mergeCloseVertices();
 
 	/**
 	 * @param container container of vertex orbite
diff --git a/include/Algo/Import/import2tablesSurface.hpp b/include/Algo/Import/import2tablesSurface.hpp
index 0110d10f85cdef3b56d0cd46555af725a2745cce..1b2275c8644b8da5595b1ef5fd54bebbc1418fbb 100644
--- a/include/Algo/Import/import2tablesSurface.hpp
+++ b/include/Algo/Import/import2tablesSurface.hpp
@@ -23,6 +23,8 @@
 *******************************************************************************/
 
 #include "Algo/Import/importPlyData.h"
+#include "Algo/Geometry/boundingbox.h"
+#include "Topology/generic/autoAttributeHandler.h"
 
 #include "openctm.h"
 
@@ -433,7 +435,7 @@ bool MeshTablesSurface<PFP>::importObj(const std::string& filename, std::vector<
 	m_emb.reserve(verticesID.size()*8);
 
 	std::vector<int> table;
-	table.reserve(64); // 64 cotes pour une face devrait suffire
+	table.reserve(64); // NBV cotes pour une face devrait suffire
 	m_nbFaces = 0;
     do
     {
@@ -931,6 +933,158 @@ bool MeshTablesSurface<PFP>::importASSIMP(const std::string& filename, std::vect
 	return true;
 }
 
+template<typename PFP>
+bool MeshTablesSurface<PFP>::mergeCloseVertices()
+{
+	const unsigned int NBV=64; // seems to be good
+
+	const int NEIGH[27]={
+	-NBV*NBV - NBV - 1, 	-NBV*NBV - NBV, 	-NBV*NBV - NBV + 1,
+	-NBV*NBV - 1, 	-NBV*NBV, 	-NBV*NBV + 1,
+	-NBV*NBV + NBV - 1,	-NBV*NBV + NBV,	- NBV*NBV + NBV + 1,
+	-NBV - 1,	- NBV,	-NBV + 1,
+	-1,	0,	1,
+	NBV - 1,	NBV,	NBV + 1,
+	NBV*NBV - NBV - 1,	NBV*NBV - NBV,	NBV*NBV - NBV + 1,
+	NBV*NBV - 1,	NBV*NBV,	NBV*NBV + 1,
+	NBV*NBV + NBV - 1,	NBV*NBV + NBV,	NBV*NBV + NBV + 1};
+
+	std::vector<unsigned int>** grid;
+	grid = new std::vector<unsigned int>*[NBV*NBV*NBV];
+
+	// init grid with null ptrs	
+	for (unsigned int i=0; i<NBV*NBV*NBV; ++i)
+		grid[i]=NULL;
+	
+	AttributeHandler<typename PFP::VEC3> positions = m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position");
+	
+	// compute BB
+	Geom::BoundingBox<typename PFP::VEC3> bb(positions[ positions.begin() ]) ;
+	for (unsigned int i = positions.begin(); i != positions.end(); positions.next(i))
+	{
+		bb.addPoint(positions[i]) ;
+	}
+
+	// add one voxel around to avoid testing		
+	typename PFP::VEC3 bbsize = (bb.max() - bb.min());
+	typename PFP::VEC3 one = bbsize/(NBV-2);
+	one*= 1.001f;
+	bb.addPoint( bb.min() - one);
+	bb.addPoint( bb.max() + one);
+	bbsize = (bb.max() - bb.min());
+	
+
+	AutoAttributeHandler<unsigned int> gridIndex(m_map,VERTEX, "gridIndex");
+	AutoAttributeHandler<unsigned int> newIndices(m_map,VERTEX, "newIndices");
+	
+	
+	// Store each vertex in the grid and store voxel index in vertex attribute
+	for (unsigned int i = positions.begin(); i != positions.end(); positions.next(i))
+	{
+		typename PFP::VEC3 P = positions[i];
+		P -= bb.min();
+		float pz = floor((P[2]/bbsize[2])*NBV);
+		float py = floor((P[1]/bbsize[1])*NBV);
+		float px = floor((P[0]/bbsize[0])*NBV);
+
+		unsigned int index = NBV*NBV*pz + NBV*py + px;
+		
+		if (pz==63) 
+			std::cout << "z 63 bb:"<<bb<<"  P="<<positions[i]<< std::endl;
+		
+		std::vector<unsigned int>* vox = grid[index];
+		if (vox==NULL)
+		{
+			grid[index] = new std::vector<unsigned int>();
+			grid[index]->reserve(8);
+			vox = grid[index];
+		}
+		vox->push_back(i);
+		gridIndex[i] = index;
+		newIndices[i] = 0xffffffff;
+	}
+	
+	// compute EPSILON: average length of 50 of 100 first edges of first faces divide by 10000
+	int nb = 100;
+	if (m_nbEdges.size()< 100) 
+		nb = m_nbEdges.size();
+		
+	int k=0;
+	typename PFP::REAL d=0;
+	for (int i=0; i< nb; i+=2)
+	{
+		typename PFP::VEC3 e1 = positions[m_emb[k+1]] - positions[m_emb[k]];
+		d += e1.norm();
+		k += m_nbEdges[i];
+	}
+	d /= float(nb/2);
+	
+	typename PFP::REAL EPSILON = d/10000.0f;
+	
+	
+	// traverse vertices
+	for (unsigned int i = positions.begin(); i != positions.end(); positions.next(i))
+	{
+		if (newIndices[i] == 0xffffffff) 
+		{
+			const typename PFP::VEC3& P = positions[i];
+			
+			for (unsigned int n=0; n<27; ++n)
+			{
+				std::vector<unsigned int>* voxel = grid[gridIndex[i]+NEIGH[n]];
+				if (voxel != NULL)
+				{
+					for (std::vector<unsigned int>::iterator v = voxel->begin(); v != voxel->end(); ++v)
+					{
+						if ((*v != i) && (*v != 0xffffffff))
+						{
+							typename PFP::VEC3 Q = positions[*v];
+							Q -= P;
+							typename PFP::REAL d2 = Q*Q;
+							if (d2 < EPSILON*EPSILON)
+							{
+								newIndices[*v] = i;
+								*v = 0xffffffff;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	// update faces indices
+	for	(std::vector<unsigned int>::iterator it = m_emb.begin(); it != m_emb.end(); ++it)
+	{
+		if (newIndices[*it] != 0xffffffff)
+		{
+			*it = newIndices[*it];
+		}
+	}
+
+	// delete embeddings
+	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;
+
+	for (unsigned int i = positions.begin(); i != positions.end(); positions.next(i))
+	{
+		if (newIndices[i] != 0xffffffff)
+		{
+			container.removeLine(i);
+		}
+	}
+
+	// release grid memory	
+	for (unsigned int i=0; i<NBV*NBV*NBV; ++i)
+		if (grid[i]!=NULL)
+			delete grid[i];
+
+	delete[] grid;
+
+	return true;
+}
+
+
+
 } // namespace Import
 
 } // namespace Algo
diff --git a/include/Algo/Import/importMesh.hpp b/include/Algo/Import/importMesh.hpp
index 22bf913399de9cc08c8a07e1782a557498670e48..08638bf98eb50a4c3b7c1b6e04c3438cf253e456 100644
--- a/include/Algo/Import/importMesh.hpp
+++ b/include/Algo/Import/importMesh.hpp
@@ -425,12 +425,15 @@ bool importMesh(typename PFP::MAP& map, MeshTablesVolume<PFP>& mtv)
 }
 
 template <typename PFP>
-bool importMesh(typename PFP::MAP& map, const std::string& filename, std::vector<std::string>& attrNames, ImportSurfacique::ImportType kind)
+bool importMesh(typename PFP::MAP& map, const std::string& filename, std::vector<std::string>& attrNames, ImportSurfacique::ImportType kind, bool mergeCloseVertices)
 {
 	MeshTablesSurface<PFP> mts(map);
 
 	if(!mts.importMesh(filename, attrNames, kind))
 		return false;
+	
+	if (mergeCloseVertices)
+		mts.mergeCloseVertices();
 
 	return importMesh<PFP>(map, mts);
 }