/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
* version 0.1                                                                  *
* Copyright (C) 2009-2011, IGG Team, LSIIT, University of Strasbourg           *
*                                                                              *
* This library is free software; you can redistribute it and/or modify it      *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 of the License, or (at your     *
* option) any later version.                                                   *
*                                                                              *
* This library 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 PURVEC3E. See the GNU Lesser General Public License  *
* for more details.                                                            *
*                                                                              *
* You should have received a copy of the GNU Lesser General Public License     *
* along with this library; if not, write to the Free Software Foundation,      *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
*                                                                              *
* Web site: http://cgogn.u-strasbg.fr/                                         *
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#include "Algo/Import/importPlyData.h"
#include "Algo/Geometry/boundingbox.h"
#include "Topology/generic/autoAttributeHandler.h"

#include "openctm.h"

#include "assimp.h"
#include "aiPostProcess.h"
#include "aiScene.h"


namespace CGoGN
{

namespace Algo
{

namespace Import
{

template<typename PFP>
ImportSurfacique::ImportType MeshTablesSurface<PFP>::getFileType(const std::string& filename)
{
	if ((filename.rfind(".trianbgz")!=std::string::npos) || (filename.rfind(".TRIANBGZ")!=std::string::npos))
		return ImportSurfacique::TRIANBGZ;

	if ((filename.rfind(".trian")!=std::string::npos) || (filename.rfind(".TRIAN")!=std::string::npos))
		return ImportSurfacique::TRIAN;

	if ((filename.rfind(".plyptm")!=std::string::npos) || (filename.rfind(".PLYGEN")!=std::string::npos))
		return ImportSurfacique::PLYPTM;

	if ((filename.rfind(".ply")!=std::string::npos) || (filename.rfind(".PLY")!=std::string::npos))
		return ImportSurfacique::PLY;

	if ((filename.rfind(".off")!=std::string::npos) || (filename.rfind(".OFF")!=std::string::npos))
		return ImportSurfacique::OFF;

	if ((filename.rfind(".obj")!=std::string::npos) || (filename.rfind(".OBJ")!=std::string::npos))
		return ImportSurfacique::OBJ;

	if ((filename.rfind(".ctm")!=std::string::npos) || (filename.rfind(".OBJ")!=std::string::npos))
		return ImportSurfacique::CTM;

	return ImportSurfacique::UNKNOWNSURFACE;
}

template<typename PFP>
bool MeshTablesSurface<PFP>::importMesh(const std::string& filename, std::vector<std::string>& attrNames, ImportSurfacique::ImportType kind)
{
	if (kind == ImportSurfacique::UNKNOWNSURFACE)
		kind = getFileType(filename);

	attrNames.clear() ;

	switch (kind)
	{
	case ImportSurfacique::TRIAN:
		CGoGNout << "TYPE: TRIAN" << CGoGNendl;
		return importTrian(filename, attrNames);
		break;
	case ImportSurfacique::TRIANBGZ:
		CGoGNout << "TYPE: TRIANBGZ" << CGoGNendl;
		return importTrianBinGz(filename, attrNames);
		break;
	case ImportSurfacique::CTM:
		CGoGNout << "TYPE: CTM" << CGoGNendl;
		return importCTM(filename, attrNames);
		break;
	case ImportSurfacique::OFF:
		CGoGNout << "TYPE: OFF" << CGoGNendl;
		return importOff(filename, attrNames);
		break;
	case ImportSurfacique::PLY:
		CGoGNout << "TYPE: PLY" << CGoGNendl;
		return importPly(filename, attrNames);
		break;
	case ImportSurfacique::PLYPTM:
		CGoGNout << "TYPE: PLYPTM" << CGoGNendl;
		return importPlyPTM(filename, attrNames);
		break;
	case ImportSurfacique::PLYPTMgeneric:
		CGoGNout << "TYPE: PLYPTMgeneric" << CGoGNendl;
		return importPlyPTMgeneric(filename, attrNames);
		break;
	case ImportSurfacique::OBJ:
		CGoGNout << "TYPE: OBJ" << CGoGNendl;
		return importObj(filename, attrNames);
		break;
	default:
		CGoGNout << "TYPE: ASSIMP" << CGoGNendl;
		return importASSIMP(filename, attrNames);
		break;
	}
	return false;
}

template<typename PFP>
bool MeshTablesSurface<PFP>::importTrian(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	// open file
	std::ifstream fp(filename.c_str(), std::ios::in);
	if (!fp.good())
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl;
		return false;
	}

	// read nb of points
	fp >> m_nbVertices;

	// read points
	std::vector<unsigned int> verticesID;
	verticesID.reserve(m_nbVertices);

	for (unsigned int i = 0; i < m_nbVertices; ++i)
	{
		VEC3 pos;
		fp >> pos[0];
		fp >> pos[1];
		fp >> pos[2];
		unsigned int id = container.insertLine();
		positions[id] = pos;
		verticesID.push_back(id);
	}

	// read nb of faces
	fp >> m_nbFaces;
	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(3*m_nbFaces);

	// read indices of faces
	for (unsigned int i = 0; i < m_nbFaces; ++i)
	{
		m_nbEdges.push_back(3);
		// read the three vertices of triangle
		int pt;
		fp >> pt;
		m_emb.push_back(verticesID[pt]);
		fp >> pt;
		m_emb.push_back(verticesID[pt]);
		fp >> pt;
		m_emb.push_back(verticesID[pt]);

		// neighbour not always good in files !!
		int neigh;
		fp >> neigh;
		fp >> neigh;
		fp >> neigh;
	}

	fp.close();
	return true;
}

template<typename PFP>
bool MeshTablesSurface<PFP>::importTrianBinGz(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	// open file
	igzstream fs(filename.c_str(), std::ios::in|std::ios::binary);

	if (!fs.good())
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl;
		return false;
	}

	// read nb of points
	fs.read(reinterpret_cast<char*>(&m_nbVertices), sizeof(int));

	// read points
	std::vector<unsigned int> verticesID;
	{	// juste pour limiter la portee des variables
		verticesID.reserve(m_nbVertices);
		float* buffer = new float[m_nbVertices*3];
		fs.read(reinterpret_cast<char*>(buffer), 3*m_nbVertices*sizeof(float));
		float *ptr = buffer;
		for (unsigned int i = 0; i < m_nbVertices; ++i)
		{
			VEC3 pos;
			pos[0]= *ptr++;
			pos[1]= *ptr++;
			pos[2]= *ptr++;

			unsigned int id = container.insertLine();
			positions[id] = pos;

			verticesID.push_back(id);
		}
		delete[] buffer;
	}

	// read nb of faces
	fs.read(reinterpret_cast<char*>(&m_nbFaces), sizeof(int));
	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(3*m_nbFaces);

	// read indices of faces
	{	// juste pour limiter la portee des variables
		int* buffer = new int[m_nbFaces*6];
		fs.read(reinterpret_cast<char*>(buffer),6*m_nbFaces*sizeof(float));
		int *ptr = buffer;

		for (unsigned int i = 0; i < m_nbFaces; i++)
		{
			m_nbEdges.push_back(3);
			m_emb.push_back(verticesID[*ptr++]);
			m_emb.push_back(verticesID[*ptr++]);
			m_emb.push_back(verticesID[*ptr++]);
		}
	}
	
	fs.close();
	return true;
}

template<typename PFP>
bool MeshTablesSurface<PFP>::importOff(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	// open file
	std::ifstream fp(filename.c_str(), std::ios::in);
	if (!fp.good())
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl;
		return false;
	}

    std::string ligne;

    // lecture de OFF
    std::getline (fp, ligne);
    if (ligne.rfind("OFF") == std::string::npos)
    {
		CGoGNerr << "Problem reading off file: not an off file" << CGoGNendl;
		CGoGNerr << ligne << CGoGNendl;
		return false;
    }

    // lecture des nombres de sommets/faces/aretes
	int nbe;
    {
    	do
    	{
    		std::getline (fp, ligne);
    	} while (ligne.size()==0);

	    std::stringstream oss(ligne);
		oss >> m_nbVertices;
		oss >> m_nbFaces;
		oss >> nbe;
    }

	//lecture sommets
	std::vector<unsigned int> verticesID;
	verticesID.reserve(m_nbVertices);
	for (unsigned int i = 0; i < m_nbVertices;++i)
	{
    	do
    	{
    		std::getline (fp, ligne);
    	} while (ligne.size() == 0);

		std::stringstream oss(ligne);

		float x,y,z;
		oss >> x;
		oss >> y;
		oss >> z;
		// on peut ajouter ici la lecture de couleur si elle existe
		VEC3 pos(x,y,z);

		unsigned int id = container.insertLine();
		positions[id] = pos;

		verticesID.push_back(id);
	}

	// lecture faces
	// normalement nbVertices*8 devrait suffire largement
	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(m_nbVertices*8);

	for (unsigned int i = 0; i < m_nbFaces; ++i)
	{
    	do
    	{
    		std::getline (fp, ligne);
    	} while (ligne.size() == 0);

		std::stringstream oss(ligne);
		int n;
		oss >> n;
		m_nbEdges.push_back(n);
		for (int j=0;j<n; ++j)
		{
			int index; // index du plongement
			oss >> index;
			m_emb.push_back(verticesID[index]);
		}
		// on peut ajouter ici la lecture de couleur si elle existe
	}

	fp.close();
	return true;
}

template <typename PFP>
bool MeshTablesSurface<PFP>::importObj(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	// open file
	std::ifstream fp(filename.c_str(), std::ios::binary);
	if (!fp.good())
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl;
		return false;
	}

//	fp.seekg(0, std::ios::end);
//	int ab = fp.tellg();
//	fp.seekg(0, std::ios::beg);
//	int ac = fp.tellg();

    std::string ligne;
    std::string tag;

    do
    {
    	fp >> tag;
    	std::getline (fp, ligne);
    }while (tag != std::string("v"));

    // lecture des sommets
	std::vector<unsigned int> verticesID;
	verticesID.reserve(102400); // on tape large (400Ko wahouuuuu !!)

	unsigned int i = 0;
    do
    {
		if (tag == std::string("v"))
		{
			std::stringstream oss(ligne);
		
			float x,y,z;
			oss >> x;
			oss >> y;
			oss >> z;

			VEC3 pos(x,y,z);

			unsigned int id = container.insertLine();
			positions[id] = pos;

			verticesID.push_back(id);
			i++;
		}

		fp >> tag;
    	std::getline(fp, ligne);
    } while (!fp.eof());

	m_nbVertices = verticesID.size();

	// close/clear/open only way to go back to beginning of file
	fp.close();
	fp.clear();
	fp.open(filename.c_str());

	do
    {
    	fp >> tag;
    	std::getline (fp, ligne);
    } while (tag != std::string("f"));

	m_nbEdges.reserve(verticesID.size()*2);
	m_emb.reserve(verticesID.size()*8);

	std::vector<int> table;
	table.reserve(64); // NBV cotes pour une face devrait suffire
	m_nbFaces = 0;
    do
    {
    	if (tag == std::string("f")) // lecture d'une face
    	{
    		std::stringstream oss(ligne);
     		table.clear();
    		while (!oss.eof())  // lecture de tous les indices
    		{
    			std::string str;
    			oss >> str;

    			unsigned int ind = 0;
    			while ( (ind<str.length()) &&  (str[ind]!='/'))
    				ind++;

				if (ind > 0)
				{
    				long index;
					std::stringstream iss(str.substr(0, ind));
					iss >> index;
		   			table.push_back(index);
				}
    		}

    		unsigned int n = table.size();
			m_nbEdges.push_back(short(n));
    		for (unsigned int j = 0; j < n; ++j)
    		{
    			int index = table[j] - 1; // les index commencent a 1 (boufonnerie d'obj ;)
				m_emb.push_back(verticesID[index]);
    		}
			m_nbFaces++;
    	}
		fp >> tag;
    	std::getline(fp, ligne);
     } while (!fp.eof());

	fp.close ();
	return true;
}

template<typename PFP>
bool MeshTablesSurface<PFP>::importPly(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	PlyImportData pid;

	if (! pid.read_file(filename) )
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl;
		return false;
	}
	
    // lecture des nombres de sommets/aretes/faces
	m_nbVertices = pid.nbVertices();
	m_nbFaces = pid.nbFaces();
 
	//lecture sommets
	std::vector<unsigned int> verticesID;
	verticesID.reserve(m_nbVertices);
	for (unsigned int i = 0; i < m_nbVertices; ++i)
	{
		VEC3 pos;
		pid.vertexPosition(i, pos);

		unsigned int id = container.insertLine();
		positions[id] = pos;

		verticesID.push_back(id);
	}

	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(m_nbVertices*8);

	for (unsigned int i = 0 ; i < m_nbFaces ; ++i)
	{
		unsigned int n = pid.getFaceValence(i);
		m_nbEdges.push_back(n);
		int* indices = pid.getFaceIndices(i);
		for (unsigned int j = 0; j < n; ++j)
		{
			m_emb.push_back(verticesID[indices[j]]);
		}
	}

	return true;
}

/**
 * Import plyPTM (K Vanhoey generic format).
 * It can handle bivariable polynomials of any degree and returns the appropriate attrNames
 * @param filename the file to import;
 * @param attrNames reference that will be filled with the attribute names
 * the number of attrNames returned depends on the degree of the polynomials :
 *  - 1 attrName for geometric position (VEC3)
 *  - 3 attrNames for local frame (3xVEC3) : Tangent, Bitangent and Normal vector
 *  - N attrNames for the function coefficients (NxVEC3) : N RGB coefficients being successively the constants, the linears (v then u), the quadratics, etc. :  : a0 + a1*v + a2*u + a3*u*v + a4*v^2 + a5*u^2.
 * N = 1 for constant polynomial,
 * N = 3 for linear polynomial,
 * N = 6 for quadratic polynomial,
 * N = 10 for cubic degree polynomial,
 * N = 15 for 4th degree polynomial,
 * ...
 * Hint : N = attrNames.size() - 4 ;
 * @return bool : success.
 */
template <typename PFP>
bool MeshTablesSurface<PFP>::importPlyPTMgeneric(const std::string& filename, std::vector<std::string>& attrNames)
{
	// Open file
	std::ifstream fp(filename.c_str(), std::ios::in) ;
	if (!fp.good())
	{
		CGoGNerr << "Unable to open file " << filename << CGoGNendl ;
		return false ;
	}

	// Read quantities : #vertices, #faces, #properties, degree of polynomials
    std::string tag ;

    fp >> tag;
	if (tag != std::string("ply")) // verify file type
	{
		CGoGNerr << filename << " is not a ply file !" <<  CGoGNout ;
		return false ;
	}

	do // go to #vertices
	{
		fp >> tag ;
	} while (tag != std::string("vertex")) ;
	unsigned int nbVertices ;
	fp >> nbVertices ; // Read #vertices

	unsigned int nbProps = 0 ;
	do // go to #faces and count #properties
	{
		fp >> tag ;
		if (tag == std::string("property"))
			++nbProps ;
	} while (tag != std::string("face")) ;
	fp >> m_nbFaces ; // Read #vertices

	do // go to end of header
	{
		fp >> tag ;
	} while (tag != std::string("end_header")) ;

	unsigned int nbCoefsPerPol = (nbProps - 12) / 3 ; // get #coefficients per polynomial

	// Define containers
	AttributeHandler<typename PFP::VEC3> positions = m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;
	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;
	attrNames.push_back(positions.name()) ;

	AttributeHandler<typename PFP::VEC3> *frame = new AttributeHandler<typename PFP::VEC3>[3] ;
	frame[0] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_T") ; // Tangent
	frame[1] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_B") ; // Bitangent
	frame[2] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_N") ; // Normal
	attrNames.push_back(frame[0].name()) ;
	attrNames.push_back(frame[1].name()) ;
	attrNames.push_back(frame[2].name()) ;

	AttributeHandler<typename PFP::VEC3> *colorPTM = new AttributeHandler<typename PFP::VEC3>[nbCoefsPerPol] ;
	for (unsigned int i = 0 ; i < nbCoefsPerPol ; ++i)
	{
		std::stringstream name ;
		name << "colorPTM_a" << i ;
		colorPTM[i] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, name.str()) ;
		attrNames.push_back(colorPTM[i].name()) ;
	}

	// Read vertices
	std::vector<unsigned int> verticesID ;
	verticesID.reserve(nbVertices) ;

	float* properties = new float[nbProps] ;
	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;
	for (unsigned int i = 0 ; i < nbVertices ; ++i) // Read and store properties for current vertex
	{
		unsigned int id = container.insertLine() ;
		verticesID.push_back(id) ;

		for (unsigned int j = 0 ; j < nbProps ; ++j) // get all properties
			fp >> properties[j] ;

		positions[id] = VEC3(properties[0],properties[1],properties[2]) ; // position
		for (unsigned int k = 0 ; k < 3 ; ++k) // frame
			for (unsigned int l = 0 ; l < 3 ; ++l)
				frame[k][id][l] = properties[3+(3*k+l)] ;
		for (unsigned int k = 0 ; k < 3 ; ++k) // coefficients
			for (unsigned int l = 0 ; l < nbCoefsPerPol ; ++l)
				colorPTM[l][id][k] = properties[12+(nbCoefsPerPol*k+l)] ;
	}
	m_nbVertices = verticesID.size() ;
	delete[] properties ;

	// Read faces index
	m_nbEdges.reserve(m_nbFaces) ;
	m_emb.reserve(3*m_nbFaces) ;
	for (unsigned int i = 0 ; i < m_nbFaces ; ++i)
	{
		// read the indices of vertices for current face
		int nbEdgesForFace ;
		fp >> nbEdgesForFace ;
		m_nbEdges.push_back(nbEdgesForFace);

		int vertexID ;
		for (int j=0 ; j < nbEdgesForFace ; ++j)
		{
			fp >> vertexID ;
			m_emb.push_back(verticesID[vertexID]);
		}
	}

	// Close file
	fp.close() ;

	return true ;
}

/**
 * Import plyPTM (F Larue format).
 * It handles only 2nd degree polynomials
 * @param filename : the file to import;
 * @param attrNames : reference that will be filled with the attribute names ;
 *  - 1 attrName for geometric position (VEC3)
 *  - 3 attrNames for local frame (3xVEC3) : Tangent, Bitangent and Normal vector
 *  - 6 attrNames for the function coefficients (6xVEC3) : 6 RGB coefficients being successively the quadratic members, the linears and the constants (u then v) : a*u^2 + b*v^2 + c*uv + d*u + e*v +f.
  * @return bool : success.
 */
template <typename PFP>
bool MeshTablesSurface<PFP>::importPlyPTM(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeHandler<typename PFP::VEC3> frame[3] ;
	frame[0] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_T") ; // Tangent
	frame[1] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_B") ; // Bitangent
	frame[2] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "frame_N") ; // Normal
	for (unsigned int i = 0 ; i < 3 ; ++i)
		attrNames.push_back(frame[i].name()) ;

	AttributeHandler<typename PFP::VEC3> colorPTM[6] ;
	colorPTM[0] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_a") ;
	colorPTM[1] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_b") ;
	colorPTM[2] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_c") ;
	colorPTM[3] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_d") ;
	colorPTM[4] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_e") ;
	colorPTM[5] = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "colorPTM_f") ;

	for (unsigned int i = 0 ; i < 6 ; ++i)
		attrNames.push_back(colorPTM[i].name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	std::ifstream fp(filename.c_str(), std::ios::binary);
	if (!fp.good())
	{
		CGoGNerr << "Unable to open file " << filename<< CGoGNendl;
		return false;
	}

    std::string ligne;
    std::string tag;

	fp >> tag;
	if (tag != std::string("ply"))
	{
		CGoGNerr <<filename<< " is not a ply file !" <<  CGoGNendl;
		return false;
	}

	// va au nombre de sommets
	do
	{
		fp >> tag;
	} while (tag != std::string("vertex"));

	unsigned int nbp;
	fp >> nbp;
	// read points
	std::vector<unsigned int> verticesID;
	verticesID.reserve(nbp);

	// va au nombre de faces en comptant le nombre de "property"
	unsigned int nb_props = 0;
	do
	{
		fp >> tag;
		if (tag == std::string("property"))
			nb_props++;
	} while (tag != std::string("face"));

	fp >> m_nbFaces;
 	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(3*m_nbFaces);

	// lecture des sommets

	// saute à la fin du header
	do
	{
		fp >> tag;
	} while (tag != std::string("end_header"));

	float* properties = new float[nb_props];

	for (unsigned int i = 0; i < nbp; ++i)
	{
		unsigned int id = container.insertLine();
		verticesID.push_back(id);

		for (unsigned int j = 0; j < nb_props; ++j)
		{
			fp >> properties[j];
		}

		positions[id] = VEC3(properties[0],properties[1],properties[2]);

		for (unsigned int k = 0 ; k < 3 ; ++k)
			for (unsigned int l = 0 ; l < 3 ; ++l)
				frame[k][id][l] = properties[3+(3*k+l)] ;

		for (unsigned int k = 0 ; k < 3 ; ++k)
			for (unsigned int l = 0 ; l < 6 ; ++l)
				colorPTM[l][id][k] = properties[12+(6*k+l)];
	}

	m_nbVertices = verticesID.size();
	delete[] properties;

// read indices of faces
	for (unsigned int i = 0; i < m_nbFaces; i++)
	{
		// read the indices vertices of face
		int nbe;
		fp >> nbe;
		m_nbEdges.push_back(nbe);

		int pt;
		for (int j=0; j<nbe; ++j)
		{
			fp >> pt;
			m_emb.push_back(verticesID[pt]);
		}
	}

	fp.close();
	return true;
}

template <typename PFP>
bool MeshTablesSurface<PFP>::importCTM(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeHandler<typename PFP::VEC3> positions =  m_map.template getAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	if (!positions.isValid())
		positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;

	attrNames.push_back(positions.name()) ;

	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;

	// Load the file using the OpenCTM API
	CTMimporter ctm;
	// Load the file
	ctm.Load(filename.c_str());

 	m_nbVertices = ctm.GetInteger(CTM_VERTEX_COUNT);

 	// read points
	std::vector<unsigned int> verticesID;
	verticesID.reserve(m_nbVertices);

	const CTMfloat* vertices = ctm.GetFloatArray(CTM_VERTICES);
	for (unsigned int i = 0; i < m_nbVertices; ++i)
	{
		VEC3 pos(vertices[0], vertices[1], vertices[2]);
		vertices += 3; // next vertex in float*
		unsigned int id = container.insertLine();
		positions[id] = pos;
		verticesID.push_back(id);
	}

	// read nb of faces
	m_nbFaces = ctm.GetInteger(CTM_TRIANGLE_COUNT);
	m_nbEdges.reserve(m_nbFaces);
	m_emb.reserve(3*m_nbFaces);

	const CTMuint * indices = ctm.GetIntegerArray(CTM_INDICES);

	// read indices of faces
	for (unsigned i=0; i<m_nbFaces; i++)
	{
		m_nbEdges.push_back(3);
		// read the three vertices of triangle
		m_emb.push_back(verticesID[*indices++]);
		m_emb.push_back(verticesID[*indices++]);
		m_emb.push_back(verticesID[*indices++]);
	}

	return true;
}

template<typename PFP>
void MeshTablesSurface<PFP>::extractMeshRec(AttributeContainer& container, AttributeHandler<typename PFP::VEC3>& positions, const struct aiScene* scene, const struct aiNode* nd, struct aiMatrix4x4* trafo)
{
	struct aiMatrix4x4 prev;

	prev = *trafo;
	aiMultiplyMatrix4(trafo,&nd->mTransformation);

	std::vector<unsigned int> verticesID;

	// foreach mesh of node
	for (unsigned int n = 0; n < nd->mNumMeshes; ++n)
	{
		const struct aiMesh* mesh = scene->mMeshes[nd->mMeshes[n]];

		verticesID.clear();
		verticesID.reserve(mesh->mNumVertices);
		//read positions
		for (unsigned int t = 0; t < mesh->mNumVertices; ++t)
		{
			// transform position
			struct aiVector3D tmp = mesh->mVertices[t];
			aiTransformVecByMatrix4(&tmp, trafo);
			// now store it
			unsigned int id = container.insertLine();
			positions[id] = VEC3(tmp[0], tmp[1], tmp[2]);
			verticesID.push_back(id);
		}
		m_nbVertices += mesh->mNumVertices;

		// read faces indices
		for (unsigned int t = 0; t < mesh->mNumFaces; ++t)
		{
			const struct aiFace* face = &mesh->mFaces[t];
			m_nbEdges.push_back(face->mNumIndices);
			for(unsigned int i = 0; i < face->mNumIndices; i++)
			{
				unsigned int pt = face->mIndices[i];
				m_emb.push_back(verticesID[pt]);
			}
		}
		m_nbFaces += mesh->mNumFaces;
	}

	// recurse on all children of node
	for (unsigned int n = 0; n < nd->mNumChildren; ++n)
	{
//		CGoGNout << "Children "<<n<< CGoGNendl;
		extractMeshRec(container, positions, scene, nd->mChildren[n], trafo);
	}
	*trafo = prev;
}

template <typename PFP>
bool MeshTablesSurface<PFP>::importASSIMP(const std::string& filename, std::vector<std::string>& attrNames)
{
	AttributeContainer& container = m_map.getAttributeContainer(VERTEX) ;
	AttributeHandler<typename PFP::VEC3> positions = m_map.template addAttribute<typename PFP::VEC3>(VERTEX, "position") ;
	attrNames.push_back(positions.name()) ;

	m_nbVertices = 0;
	m_nbFaces = 0;

 	m_nbEdges.reserve(5000);
	m_emb.reserve(15000);

	struct aiMatrix4x4 trafo;
	aiIdentityMatrix4(&trafo);

	m_lab = 0;
	const struct aiScene* scene = aiImportFile(filename.c_str(), aiProcess_FindDegenerates | aiProcess_JoinIdenticalVertices);
	extractMeshRec(container, positions, scene, scene->mRootNode, &trafo);

	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

} // namespace CGoGN