#include <stdlib.h>
#include "SimplificationMethod.h"
#include <iostream>
#include <fstream>

using namespace std;
using namespace Geometry;

#include <gfx/geom/ProxGrid.h>
#include "quadrics.h"

//---------------------------------------------------------------------------
//	Constructor
//---------------------------------------------------------------------------
SimplificationMethod::SimplificationMethod(const Mesh *m) 
{
	objmesh							=	m;
	mGeoMesh						=	NULL;
	first_index_submesh	=	new unsigned int[objmesh->mSubMeshCount];
	indexMeshLeaves			=	-1;
}

//---------------------------------------------------------------------------
//	Destructor
//---------------------------------------------------------------------------
SimplificationMethod::~SimplificationMethod()
{
	delete heap;
	delete [] first_index_submesh;
}

//---------------------------------------------------------------------------
//	Recalculate the cost of the vertices of an edge
//---------------------------------------------------------------------------
void SimplificationMethod::compute_pair_info(simplif::pair_info *auxpair)
{
	simplif::Vertex *v0 = auxpair->v0;
	simplif::Vertex *v1 = auxpair->v1;

	simplif::vert_info& v0_info = vertex_info(v0);
	simplif::vert_info& v1_info = vertex_info(v1);

	simplif::Mat4 Q = v0_info.Q + v1_info.Q;
	simplif::real norm = v0_info.norm + v1_info.norm;

	auxpair->cost = simplif::quadrix_pair_target(Q, v0, v1, auxpair->candidate);

	if( simplif::will_weight_by_area )
		auxpair->cost /= norm;

	if( simplif::will_preserve_mesh_quality )
		auxpair->cost += pair_mesh_penalty(M0, v0, v1, auxpair->candidate);


	//
	// NOTICE:  In the heap we use the negative cost.  That's because
	//          the heap is implemented as a MAX heap.
	//
	if( auxpair->isInHeap() )
	{
		heap->update(auxpair, (float)-auxpair->cost);
	}
	else
	{
		heap->insert(auxpair, (float)-auxpair->cost);
	}
}

//---------------------------------------------------------------------------
//	Reasign the new vertex to the adequate face
//---------------------------------------------------------------------------
int SimplificationMethod::predict_face(simplif::Face& F, simplif::Vertex *v1, simplif::Vertex *v2, simplif::Vec3& vnew, simplif::Vec3& f1, simplif::Vec3& f2, simplif::Vec3& f3)
{
	int nmapped = 0;

	if( F.vertex(0) == v1 || F.vertex(0) == v2 )
	{ f1 = vnew;  nmapped++; }
	else f1 = *F.vertex(0);

	if( F.vertex(1) == v1 || F.vertex(1) == v2 )
	{ f2 = vnew;  nmapped++; }
	else f2 = *F.vertex(1);

	if( F.vertex(2) == v1 || F.vertex(2) == v2 )
	{ f3 = vnew;  nmapped++; }
	else f3 = *F.vertex(2);

	return nmapped;
}

#define MESH_INVERSION_PENALTY 1e9

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
simplif::real SimplificationMethod::pair_mesh_penalty(simplif::Model& M, simplif::Vertex *v1, simplif::Vertex *v2, simplif::Vec3& vnew)
{
	static simplif::face_buffer changed;
	changed.reset();

	M.contractionRegion(v1, v2, changed);

	// real Nsum = 0;
	simplif::real Nmin = 0;

	for(int i=0; i<changed.length(); i++)
	{
		simplif::Face& F = *changed(i);
		simplif::Vec3 f1, f2, f3;

		int nmapped = predict_face(F, v1, v2, vnew, f1, f2, f3);

		// Only consider non-degenerate faces
		if( nmapped < 2 )
		{
			simplif::Plane Pnew(f1, f2, f3);
			simplif::real delta =  Pnew.normal() * F.plane().normal();

			if( Nmin > delta ) Nmin = delta;
		}
	}

	//return (-Nmin) * MESH_INVERSION_PENALTY;
	if( Nmin < 0.0 )
		return MESH_INVERSION_PENALTY;
	else
		return 0.0;
}

//---------------------------------------------------------------------------
//	Returns true if the givens vertices are a valid pair
//---------------------------------------------------------------------------
bool SimplificationMethod::check_for_pair(simplif::Vertex *v0, simplif::Vertex *v1)
{
	const simplif::pair_buffer& pairs = vertex_info(v0).pairs;

	for(int i=0; i<pairs.length(); i++)
	{
		if( pairs(i)->v0==v1 || pairs(i)->v1==v1 )
			return true;
	}

	return false;
}

//---------------------------------------------------------------------------
//	Create a new pair with two given vertices
//---------------------------------------------------------------------------
simplif::pair_info *SimplificationMethod::new_pair(simplif::Vertex *v0, simplif::Vertex *v1)
{
	simplif::vert_info& v0_info = vertex_info(v0);
	simplif::vert_info& v1_info = vertex_info(v1);

	simplif::pair_info *pair = new simplif::pair_info(v0,v1);
	v0_info.pairs.add(pair);
	v1_info.pairs.add(pair);

	return pair;
}

//---------------------------------------------------------------------------
//	Remove a given pair
//---------------------------------------------------------------------------
void SimplificationMethod::delete_pair(simplif::pair_info *auxpair)
{
	simplif::vert_info& v0_info = vertex_info(auxpair->v0);
	simplif::vert_info& v1_info = vertex_info(auxpair->v1);

	v0_info.pairs.remove(v0_info.pairs.find(auxpair));
	v1_info.pairs.remove(v1_info.pairs.find(auxpair));

	if( auxpair->isInHeap() )
		heap->kill(auxpair->getHeapPos());

	delete auxpair;
}

//---------------------------------------------------------------------------
//	Contract a given pair
//---------------------------------------------------------------------------
void SimplificationMethod::do_contract(simplif::Model& m, simplif::pair_info *auxpair)
{
	simplif::Vertex *v0 = auxpair->v0;  simplif::Vertex *v1 = auxpair->v1;
	simplif::vert_info& v0_info = vertex_info(v0);
	simplif::vert_info& v1_info = vertex_info(v1);
	simplif::Vec3 vnew = auxpair->candidate;
	int i;

	// Make v0 be the new vertex
	v0_info.Q += v1_info.Q;
	v0_info.norm += v1_info.norm;

	simplif::Vec3 p(vnew);
	simplif::Vec3 v0antes(v0->operator[](simplif::X),v0->operator[](simplif::Y),v0->operator[](simplif::Z));

	// Find the edge to remove (before remove it (contract()))
	simplif::Edge *econ=NULL;
	bool finded = false;
	for (i=0; i<v0->edgeUses().length(); i++)
	{		
		econ = v0->edgeUses()(i);
		if (v1->vID == econ->org()->vID || v1->vID == econ->dest()->vID)
		{
			finded = true;
			break;
		}
	}

	assert(finded);


	// Perform the actual contraction
	static simplif::face_buffer changed;
	changed.reset();
	m.contract(v0, v1, vnew, changed);

	// Stores the simplification step
	MeshSimplificationSequence::Step simplifstep;
	simplif::Vec3 vdesp = p - v0antes;

	// Stores the displacement in the simplification step
	simplifstep.x=(float)vdesp[simplif::X];
	simplifstep.y = (float)vdesp[simplif::Y];
	simplifstep.z = (float)vdesp[simplif::Z];
	
	//	Submeshes which pertains each vertex
	simplifstep.mV1 = v0->validID(); 

	//	The simplification method returns as v1 the dead vertex.
	simplifstep.mV0 = v1->validID();

	//	Number of triangles that are removed in this simplification step.
	int _numDelTris = 0; // 1 or 2 triangles can be removed.

	for (i = 0;	i < changed.length();	i++)
	{
		if (!changed(i)->isValid())
		{
			_numDelTris++;
		}
	}

	int del_index	=	0;

	// mModfaces y mT0, mT1 of simplifstep stores the triangles id's -> geomesh has not triangles id's
	for (i=0; i<changed.length(); i++)
	{
		simplif::Face *auxface = changed(i);

		if (auxface->isValid())
		{
			// Modified triangles
			simplifstep.mModfaces.push_back(auxface->validID());
		}
		else
		{
			// Removed triangles
			if (_numDelTris==1) 
			{
				simplifstep.mT0=auxface->validID();
				simplifstep.mT1=auxface->validID();
			}
			else
			{
				if(del_index==0)
				{
					simplifstep.mT0=auxface->validID();
					del_index++;
				}
				else
				{
					simplifstep.mT1=auxface->validID();
				}
			}
		}
	}

	decim_data.push_back(simplifstep); // Simplification sequence that make reference to the simplification method sequence

#ifdef SUPPORT_VCOLOR
	//
	// If the vertices are colored, color the new vertex
	// using the average of the old colors.
	v0->props->color += v1->props->color;
	v0->props->color /= 2;
#endif

	// Remove the pair that we just contracted
	delete_pair(auxpair);

	// Recalculate pairs associated with v0
	for(i=0; i<v0_info.pairs.length(); i++)
	{
		simplif::pair_info *p = v0_info.pairs(i);
		compute_pair_info(p);
	}

	// Process pairs associated with now dead vertex
	static simplif::pair_buffer condemned(6); // collect condemned pairs for execution
	condemned.reset();

	for(i=0; i<v1_info.pairs.length(); i++)
	{
		simplif::pair_info *p = v1_info.pairs(i);

		simplif::Vertex *u;
		if( p->v0 == v1 )      
			u = p->v1;
		else 
			if( p->v1 == v1)  
				u = p->v0;
			else
				std::cerr << "YOW!  This is a bogus pair." << endl;

		if( !check_for_pair(u, v0) )
		{
			p->v0 = v0;
			p->v1 = u;
			v0_info.pairs.add(p);
			compute_pair_info(p);
		}
		else
			condemned.add(p);
	}

	for(i=0; i<condemned.length(); i++)
		// Do you have any last requests?
		delete_pair(condemned(i));
	v1_info.pairs.reset(); // safety precaution

}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
bool SimplificationMethod::decimate_quadric(simplif::Vertex *v, simplif::Mat4& Q)
{
	if( vinfo.length() > 0 )
	{
		Q = vinfo(v->uniqID).Q;
		return true;
	}
	else
		return false;
}

//---------------------------------------------------------------------------
//	Extract a pair from the heap, concract it and remove it form the heap
//---------------------------------------------------------------------------
void SimplificationMethod::decimate_contract(simplif::Model& m)
{
	simplif::heap_node *top;
	simplif::pair_info *pair;

	for(;;)
	{
		top = heap->extract();
		if( !top ) return;
		pair = (simplif::pair_info *)top->obj;

		//Remove all the vertices of the removed edges
		bool sharededge= false;
		for (int i=0; i<pair->v0->edgeUses().length(); i++)
		{		
			simplif::Edge *econ = pair->v0->edgeUses()(i);
			if (pair->v1->vID == econ->org()->vID || pair->v1->vID == econ->dest()->vID)
			{
				sharededge= true;
				break;
			}

		}

		if( pair->isValid() && sharededge)
			break;

		delete_pair(pair);
	}

	do_contract(m, pair);

	M0.validVertCount--;  // Attempt to maintain valid vertex information
}

//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
simplif::real SimplificationMethod::decimate_error(simplif::Vertex *v)
{
	simplif::vert_info& info = vertex_info(v);

	simplif::real err = simplif::quadrix_evaluate_vertex(*v, info.Q);

    if( simplif::will_weight_by_area )
		err /= info.norm;

    return err;
}

//---------------------------------------------------------------------------
//	Extract the minimum cost of the the valid nodes (not simplified)
//---------------------------------------------------------------------------
simplif::real SimplificationMethod::decimate_min_error()
{
	simplif::heap_node *top;
	simplif::pair_info *pair;

	for(;;)
	{
		top = heap->top();
		if( !top ) return -1.0;
		pair = (simplif::pair_info *)top->obj;

		if( pair->isValid() )
			break;

		top = heap->extract();
		delete_pair(pair);
	}

	return pair->cost;
}

//---------------------------------------------------------------------------
//	Returns the maximum error by vertex
//---------------------------------------------------------------------------
simplif::real SimplificationMethod::decimate_max_error(simplif::Model& m)
{
	simplif::real max_err = 0;

	for(int i=0; i<m.vertCount(); i++)
		if( m.vertex(i)->isValid() )
		{
			max_err = MAX(max_err, decimate_error(m.vertex(i)));
		}

	return max_err;
}

//---------------------------------------------------------------------------
//	Initializations
//---------------------------------------------------------------------------
void SimplificationMethod::decimate_init(simplif::Model& m, simplif::real limit)
{
	int i,j;

	vinfo.init(m.vertCount());

	if( simplif::will_use_vertex_constraint )
		for(i=0; i<m.vertCount(); i++)
		{
			simplif::Vertex *v = m.vertex(i);
			if( v->isValid() )
				vertex_info(v).Q = simplif::quadrix_vertex_constraint(*v);
		}

	for(i=0; i<m.faceCount(); i++)
		if( m.face(i)->isValid() )
		{
			if( simplif::will_use_plane_constraint )
			{
				simplif::Mat4 Q = simplif::quadrix_plane_constraint(*m.face(i));
				simplif::real norm = 0.0;

				if( simplif::will_weight_by_area )
				{
					norm = m.face(i)->area();
					Q *= norm;
				}

				for(j=0; j<3; j++)
				{
					vertex_info(m.face(i)->vertex(j)).Q += Q;
					vertex_info(m.face(i)->vertex(j)).norm += norm;

				}
			}
		}

	if( simplif::will_constrain_boundaries )
	{
		for(i=0; i<m.edgeCount(); i++)
			if( m.edge(i)->isValid() && simplif::check_for_discontinuity(m.edge(i)) )
			{
				simplif::Mat4 B = simplif::quadrix_discontinuity_constraint(m.edge(i));
				simplif::real norm = 0.0;

				if( simplif::will_weight_by_area )
				{
					norm = simplif::norm2(*m.edge(i)->org() - *m.edge(i)->dest());
					B *= norm;
				}

				B *= simplif::boundary_constraint_weight;
				vertex_info(m.edge(i)->org()).Q += B;
				vertex_info(m.edge(i)->org()).norm += norm;
				vertex_info(m.edge(i)->dest()).Q += B;
				vertex_info(m.edge(i)->dest()).norm += norm;
			}
	}

	heap = new simplif::Heap(m.validEdgeCount);

	int pair_count = 0;

	for(i=0; i<m.edgeCount(); i++)
		if( m.edge(i)->isValid() )
		{
			simplif::pair_info *pair = new_pair(m.edge(i)->org(), m.edge(i)->dest());
			//pointers_to_remove.push_back(pair);
			compute_pair_info(pair);
			pair_count++;
		}

	if( limit<0 )
	{
		limit = m.bounds.radius * 0.05;
	}
	proximity_limit = limit * limit;
	if( proximity_limit > 0 )
	{
		simplif::ProxGrid grid(m.bounds.min, m.bounds.max, limit);
		for(i=0; i<m.vertCount(); i++)
			grid.addPoint(m.vertex(i));

		simplif::buffer<simplif::Vec3 *> nearby(32);
		for(i=0; i<m.vertCount(); i++)
		{
			nearby.reset();
			grid.proximalPoints(m.vertex(i), nearby);

			for(j=0; j<nearby.length(); j++)
			{
				simplif::Vertex *v1 = m.vertex(i);
				simplif::Vertex *v2 = (simplif::Vertex *)nearby(j);

				if( v1->isValid() && v2->isValid() )
				{
#ifdef SAFETY
					assert(pair_is_valid(v1,v2));
#endif
					if( !check_for_pair(v1,v2) )
					{
						simplif::pair_info *pair = new_pair(v1,v2);
						compute_pair_info(pair);
						pair_count++;
					}
				}

			}
		}
	}
}

//---------------------------------------------------------------------------
//	Initilizations
//---------------------------------------------------------------------------
void SimplificationMethod::simplifmethod_init(void)
{
	int i;

	// Change mesh structure.
	geomesh2simplifModel();

	M0.bounds.complete();

	initialVertCount = M0.vertCount();
	initialEdgeCount = M0.edgeCount();
	initialFaceCount = M0.faceCount();

	// Get rid of degenerate faces
	for(i=0; i<M0.faceCount(); i++)
		if( !M0.face(i)->plane().isValid() )
			M0.killFace(M0.face(i));

	M0.removeDegeneracy(M0.allFaces());

	// Get rid of unused vertices
	for(i=0; i<M0.vertCount(); i++)
	{
		if( M0.vertex(i)->edgeUses().length() == 0 )
			M0.vertex(i)->kill();
	}
}

//---------------------------------------------------------------------------
//	Do the contractions till the required LOD (percentage option)
//---------------------------------------------------------------------------
void SimplificationMethod::simplifmethod_run(int finalfaces,TIPOFUNC upb)
{
	decimate_init(M0, simplif::pair_selection_tolerance);
	unsigned	int	contador	=	0;
	float					percent; // simplification percentage.

	percent	=	(float)(2.0 * 100.0) / (M0.validFaceCount - finalfaces);

	//	Update progress bar.
	if (upb)
	{
		upb(0.0);
	}

	for (std::vector<MeshSimplificationSequence::Step>::iterator it= decim_data.begin(); it!=decim_data.end(); it++)
	{
		it->mModfaces.clear();
	}
	decim_data.clear();

	while( M0.validFaceCount > finalfaces/*simplif::face_target*/
			&& M0.validFaceCount > 1
			&& decimate_min_error() < simplif::error_tolerance ) 
	{
		decimate_contract(M0);

		//	Update progress bar.
		if (upb)
		{
			upb(percent);
		}
	}

	//	Update progress bar.
	if (upb)
	{
		upb(100.0);
	}
}

//---------------------------------------------------------------------------
//	Do the contractions till the required LOD (number of vertices option)
//---------------------------------------------------------------------------
void SimplificationMethod::simplifmethod_runv(int numvertices,TIPOFUNC upb)
{
	decimate_init(M0, simplif::pair_selection_tolerance);
	unsigned	int	contador	=	0;
	int						previousValidVertCount;
	float					percent; // Simplification percentage.

	previousValidVertCount	=	0;

	percent	=	(float)(2.0 * 100.0) /(float)((M0.validFaceCount - numvertices / 3.0));

	//	Update progress bar.
	if (upb)
	{
		upb(0.0);
	}

	for (std::vector<MeshSimplificationSequence::Step>::iterator it= decim_data.begin(); it!=decim_data.end(); it++)
	{
		it->mModfaces.clear();
	}
	decim_data.clear();

	while(	M0.validVertCount > numvertices	/*simplif::face_target*/
			&&
			decimate_min_error() < simplif::error_tolerance
			&&
			previousValidVertCount	!=	M0.validVertCount)
	{
		previousValidVertCount	=	M0.validVertCount;

		decimate_contract(M0);

		//	Update progress bar.
		if (upb)
		{
			upb(percent);
		}
	}

	//	Update progress bar.
	if (upb)
	{
		upb(100.0);
	}
}

//---------------------------------------------------------------------------
// Class to create a map without repeated vertices
//---------------------------------------------------------------------------
class _float3_
{
	public:
		float x,y,z;
		_float3_(float x=0.0f, float y=0.0f, float z=0.0f){ this->x = x; this->y = y; this->z = z; }
		_float3_(const _float3_ &f){ x=f.x; y=f.y; z=f.z; }
		_float3_ & operator=(const _float3_ &f){ x=f.x; y=f.y; z=f.z; return *this; }
		bool operator<(const _float3_ &f) const 
		{ 
			if (x<f.x) return true;
			if (x>f.x) return false;
			if (y<f.y) return true;
			if (y>f.y) return false;
			if (z<f.z) return true;
			if (z>f.z) return false;
			return false;
		}
};


//---------------------------------------------------------------------------
//Unify the model to work with it
//The vertices that have the same space coordinates are stored as an unique vertex
//in order to avoid holes in the simplification
//The necessary information to store correctly the decimation data is stored
//---------------------------------------------------------------------------
void SimplificationMethod::geomesh2simplifModel(void)
{
	std::map<_float3_,int> vertices_map; //Unique vertices (without repetitions) from all the submeshes
	unsigned int num_final_verts = 0; //Number of unique vertices from all the submeshes

	int num_final_verts2=0;
	int num_repeated_verts=0;

	int simplifindex;
	for (size_t i=0; i<objmesh->mSubMeshCount; i++) //For all the submeshes
	{
		if (i	!=	indexMeshLeaves)
		{
			//If the model has texture information
			if (objmesh->mSubMesh[i].mVertexBuffer->mVertexInfo & VERTEX_TEXCOORDS) 
			{
				//For all the vertices of each submesh
				for (size_t j=0; j<objmesh->mSubMesh[i].mVertexBuffer->mVertexCount; j++) 
				{
					_float3_ auxvert(objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x, 
							objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y,
							objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z);

					std::map<_float3_,int>::iterator mit = vertices_map.find(auxvert);
					if (mit==vertices_map.end())
					{
						//If the vertex is not in vertices_map, then is added
						simplifindex=M0.in_Vertex(simplif::Vec3(objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x, 
									objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y,
									objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z));
						M0.in_Normal(simplif::Vec3(objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].x, 
									objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].y,
									objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].z));

						M0.in_TexCoord(simplif::Vec2(objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].x, 
									objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].y));
						vertices_map[auxvert] = num_final_verts;
						num_final_verts ++;

						std::vector<int> v;
						submeshmap[simplifindex]=v;
						submeshmap[simplifindex].push_back((int)i); //Add the index of the mesh that contains the vertex
						std::vector<int> v2;
						vertexbuffermap[simplifindex]=v2;
						vertexbuffermap[simplifindex].push_back((int)j); //Add the index of vertex in the submesh

						// Debug
						num_final_verts2++;
					}
					else
					{
						//the vertex is already in vertices_map, so store the reference of the first one with the same space coodinates
						simplifindex = mit->second;
						submeshmap[simplifindex].push_back((int)i); //Submeshes that contains the vertex
						vertexbuffermap[simplifindex].push_back((int)j); //Indices of the vertex buffers that contains the vertex

						// Debug
						num_repeated_verts++;
					}	
				}
			}
			else //the model has not texture information
			{ 
				//For all the vertices of each submesh
				for (size_t j=0; j<objmesh->mSubMesh[i].mVertexBuffer->mVertexCount; j++) 
				{
					_float3_ auxvert(objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x, 
							objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y,
							objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z);

					std::map<_float3_,int>::iterator mit = vertices_map.find(auxvert);
					if (mit==vertices_map.end())
					{
						//If the vertex is not in vertices_map, then is added
						simplifindex=M0.in_Vertex(simplif::Vec3(objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x, 
									objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y,
									objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z));
						M0.in_Normal(simplif::Vec3(objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].x, 
									objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].y,
									objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].z));
						M0.in_TexCoord(simplif::Vec2(0.0,0.0));
						vertices_map[auxvert] = num_final_verts;
						num_final_verts ++;

						std::vector<int> v;
						submeshmap[simplifindex]=v;
						submeshmap[simplifindex].push_back((int)i); //Add the index of the mesh that contains the vertex
						std::vector<int> v2;
						vertexbuffermap[simplifindex]=v2;
						vertexbuffermap[simplifindex].push_back((int)j); //Add the index of vertex in the submesh
					}
					else
					{
						//the vertex is already in vertices_map, so store the reference of the first one with the same space coodinates
						simplifindex = mit->second;
						submeshmap[simplifindex].push_back((int)i); //Submeshes that contains the vertex
						vertexbuffermap[simplifindex].push_back((int)j); //Indices of the vertex buffers that contains the vertex
					}	
				}
			}
			if (objmesh->mSubMesh[i].mSharedVertexBuffer)
			{
				printf("Shared break");
				break;
			}
		}
	}

	//Create the faces
	unsigned int base_index=0;

	number_of_triangles=0; //Number of triangles of the model
	int face_index=0;

	int index_in_vertices_map0,index_in_vertices_map1,index_in_vertices_map2;
	for (size_t i=0; i<objmesh->mSubMeshCount; i++) //For all the submeshes
	{ 
		if (i!=indexMeshLeaves)
		{
			for (Index j=0; j<objmesh->mSubMesh[i].mIndexCount; j+=3) 
			{
				//+1 is required because the first index in the simplification method is 1
				//Index of the vertex in vertices_map
				Index indice0=objmesh->mSubMesh[i].mIndex[j];
				_float3_ auxvert0(objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice0].x, 
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice0].y,
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice0].z);
				index_in_vertices_map0=vertices_map[auxvert0] + 1;

				Index indice1=objmesh->mSubMesh[i].mIndex[j+1];
				_float3_ auxvert1(objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice1].x, 
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice1].y,
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice1].z);
				index_in_vertices_map1=vertices_map[auxvert1] + 1;

				Index indice2=objmesh->mSubMesh[i].mIndex[j+2];
				_float3_ auxvert2(objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice2].x, 
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice2].y,
						objmesh->mSubMesh[i].mVertexBuffer->mPosition[indice2].z);
				index_in_vertices_map2=vertices_map[auxvert2] + 1;

				//Create a triangle with its indices in vertices_map
				face_index=M0.miin_Face(index_in_vertices_map0,index_in_vertices_map1,index_in_vertices_map2);

				//igeo allows to identify the submesh that contains the triangle
				M0.face(face_index)->igeo=(int)i;
			}
		}
	}
	number_of_triangles=face_index;
}


//---------------------------------------------------------------------------
//Wrties a file with the vertices, triangles and the simplification sequence
//---------------------------------------------------------------------------
void SimplificationMethod::WriteOBJ(void)
{
	MeshSimplificationSequence::Step stepaux;
	vector<int> v0,v1,submesh0,submesh1;
	std::ofstream obj("salida.lod");
	obj << "begin" << std::endl;
	//vertices
	for (size_t i=0; i<objmesh->mSubMeshCount; i++) 
	{
		for (size_t j=0; j<objmesh->mSubMesh[i].mVertexBuffer->mVertexCount; j++) 
		{	
			obj << "v " <<  objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x << " " <<
				objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y << " " <<
				objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z << " " << std::endl;
		}
	}


	//faces
	for (size_t i=0; i<objmesh->mSubMeshCount; i++) 
	{
		for (size_t j=0; j<objmesh->mSubMesh[i].mIndexCount; j=j+3) 
		{	
			obj << "f " <<  objmesh->mSubMesh[i].mIndex[j]+1 << " " << 
				objmesh->mSubMesh[i].mIndex[j+1]+1 << " " << 
				objmesh->mSubMesh[i].mIndex[j+2]+1 << std::endl;
		}
	}

	obj << "end" << std::endl;

	for (unsigned i=0; i<decim_data.size(); i++) 
	{
		stepaux = decim_data.operator [](i);
		v0=vertexbuffermap[stepaux.mV0];  //vertexbuffer index
		submesh0=submeshmap[stepaux.mV0]; //displacement for each submesh
		v1=vertexbuffermap[stepaux.mV1];
		submesh1=submeshmap[stepaux.mV1];
		for (unsigned j=0; j<v0.size(); j++)
		{
			for (unsigned k=0; k<v1.size(); k++) 
			{
				//obtain the index of the VertexBuffer of each submesh and the displacemente by submesh
				obj << "v% " << v1.operator [](k)+first_index_submesh[submesh1.operator [](k)] << 
					" " << v0.operator [](j)+first_index_submesh[submesh0.operator [](j)] << " " << 
					stepaux.x << " " << stepaux.y << " " << stepaux.z << " " << stepaux.mT0 << " " <<
					stepaux.mT1 << " & ";
				if (stepaux.mModfaces.size()>0) 
				{
					unsigned int ii=0;
					for (ii=0; ii<stepaux.mModfaces.size()-1; ii++) 
					{
						obj << stepaux.mModfaces.operator [](ii) << " ";
					}
					obj << stepaux.mModfaces.operator [](ii) << std::endl;
				}
				else 
				{ 
					obj << std::endl; 
				}					
			}
		}
	}
	obj.close();
}

//---------------------------------------------------------------------------
//simplifmethod_init is executed to do some initalizations
//simplifmethod_run is executed to do the decimation
//Then, the simpified model is created and the decimation information is returned
//---------------------------------------------------------------------------
MeshSimplificationSequence *SimplificationMethod::Decimate(	float lod,
																														int simpliftype,
																														TIPOFUNC upb)
{
	simplifmethod_init();

	if (simpliftype == 0)
	{
		//Percentage option
		simplifmethod_run((int)(number_of_triangles*lod),upb);
	}
	else
	{
		//Number of vertices option
		simplifmethod_runv((int)lod,upb);
	}

	//	Store unique vertices by submesh (without repetitions).
	typedef std::map<int,int> MAPAINDIND;

	//	Store all the indices by submesh (with repetitions).
	typedef std::vector<int> REPINDLIST;

	MAPAINDIND *unique_verts_inds_by_geo = new MAPAINDIND[objmesh->mSubMeshCount]; 
	REPINDLIST *ver_inds_rep_by_geo 	= new REPINDLIST[objmesh->mSubMeshCount]; 

	//	Index counter by submesh.
	int *inextvert = new int[objmesh->mSubMeshCount];

	for (unsigned int	i = 0;	i < objmesh->mSubMeshCount;	i++)
	{
		inextvert[i] = 0;
	}


	//	Construct the model.
	for (int	i = 0;	i < M0.faceCount();	i++)
	{
		if (M0.face(i)->isValid())
		{
			simplif::Face *auxface	=	M0.face(i);

			//	Determine to which submesh pertains the triangle.
			int igeo = auxface->igeo;

			//	Insert to each submesh its triangles.
			std::map<int,int>::iterator findit = unique_verts_inds_by_geo[igeo].find(auxface->vertex(0)->validID());
			if (findit == unique_verts_inds_by_geo[igeo].end())
			{
				//	If it is not added... add and map it.
				unique_verts_inds_by_geo[igeo][auxface->vertex(0)->validID()] = inextvert[igeo];
				inextvert[igeo]++;
			}			
			findit = unique_verts_inds_by_geo[igeo].find(auxface->vertex(1)->validID());
			if (findit == unique_verts_inds_by_geo[igeo].end())
			{
				//	If it is not added... add and map it.
				unique_verts_inds_by_geo[igeo][auxface->vertex(1)->validID()] = inextvert[igeo];
				inextvert[igeo]++;
			}			
			findit = unique_verts_inds_by_geo[igeo].find(auxface->vertex(2)->validID());
			if (findit == unique_verts_inds_by_geo[igeo].end())
			{
				//	If it is not added... add and map it.
				unique_verts_inds_by_geo[igeo][auxface->vertex(2)->validID()] = inextvert[igeo];
				inextvert[igeo]++;
			}			

			//	Total number of indices by submesh.
			ver_inds_rep_by_geo[igeo].push_back(auxface->vertex(0)->validID());
			ver_inds_rep_by_geo[igeo].push_back(auxface->vertex(1)->validID());
			ver_inds_rep_by_geo[igeo].push_back(auxface->vertex(2)->validID());
		}
	}

	//	Output simplified mesh.
	mGeoMesh								=	new Mesh();
	mGeoMesh->mVertexBuffer	=	new	VertexBuffer();

	mGeoMesh->mSubMeshCount = objmesh->mSubMeshCount;

	//	Memory allocation for the submeshes.
	mGeoMesh->mSubMesh = new SubMesh[objmesh->mSubMeshCount];

	//	Fill up bounding box settings.
	mGeoMesh->mMeshBounds.maxX					=	objmesh->mMeshBounds.maxX;
	mGeoMesh->mMeshBounds.maxY					=	objmesh->mMeshBounds.maxY;
	mGeoMesh->mMeshBounds.maxZ					=	objmesh->mMeshBounds.maxZ;
	mGeoMesh->mMeshBounds.minX					=	objmesh->mMeshBounds.minX;
	mGeoMesh->mMeshBounds.minY					=	objmesh->mMeshBounds.minY;
	mGeoMesh->mMeshBounds.minZ					=	objmesh->mMeshBounds.minZ;
	mGeoMesh->mMeshBounds.radius				=	objmesh->mMeshBounds.radius;
	mGeoMesh->mMeshBounds.scaleFactor		=	objmesh->mMeshBounds.scaleFactor;

	//	Copy skeleton name.
	if (objmesh->hasSkeleton)
	{
		mGeoMesh->hasSkeleton	=	true;

		strcpy(mGeoMesh->mSkeletonName,objmesh->mSkeletonName);
	}

	//	Copy mesh bones.
	if (!objmesh->mBones.empty())
	{
		for (unsigned int j = 0; j < objmesh->mBones.size(); j++)
		{
			mGeoMesh->mBones.push_back(objmesh->mBones[j]);
		}
	}

	bool copiedShared = false;

	for (size_t	i = 0;	i < objmesh->mSubMeshCount;	i++)
	{
		mGeoMesh->mSubMesh[i].mStripCount					=	0;
		mGeoMesh->mSubMesh[i].mStrip							=	NULL;

		mGeoMesh->
			mSubMesh[i].mSharedVertexBuffer	=	false;

		strcpy(	mGeoMesh->mSubMesh[i].mMaterialName,
				objmesh->mSubMesh[i].mMaterialName);

		//	Copy submesh bones.
		if (!objmesh->mSubMesh[i].mBones.empty())
		{
			for (unsigned int j = 0; j < objmesh->mSubMesh[i].mBones.size(); j++)
			{
				mGeoMesh->mSubMesh[i].mBones.push_back(objmesh->
						mSubMesh[i].mBones[j]);
			}
		}

		if (i != indexMeshLeaves)
		{
			//	Indices vectors.
			mGeoMesh->mSubMesh[i].mIndexCount	=	ver_inds_rep_by_geo[i].size();
			mGeoMesh->mSubMesh[i].mIndex			=	new Index[mGeoMesh->
				mSubMesh[i].mIndexCount];

			//	Store the indices.
			for (unsigned int	j = 0;	j < mGeoMesh->mSubMesh[i].mIndexCount;	j++)
			{
				//	Obtain the indices that point at VertexBuffer.
				mGeoMesh->mSubMesh[i].mIndex[j]	=	unique_verts_inds_by_geo[i].operator [](ver_inds_rep_by_geo[i].operator [](j));
			}

			//	Copy the vertices.
			if (mGeoMesh->mSubMesh[i].mSharedVertexBuffer)
			{
				//	Shared vertices.
				if (!copiedShared)
				{
					mGeoMesh->mVertexBuffer = new VertexBuffer();
					copiedShared 						= true;
				}

				mGeoMesh->mSubMesh[i].mVertexBuffer = mGeoMesh->mVertexBuffer;
			}
			else
			{
				//	There's no shared vertices.
				mGeoMesh->mSubMesh[i].mVertexBuffer	=	new VertexBuffer();
			}

			mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount	=	unique_verts_inds_by_geo[i].size();
			mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexInfo	=	objmesh->mSubMesh[i].mVertexBuffer->mVertexInfo;

			//	Allocate memory for position, normal and texutre coordinates.
			mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition	=	new Vector3[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];
			mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal		=	new Vector3[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];
			mGeoMesh->mSubMesh[i].mVertexBuffer->mTexCoords	=	new Vector2[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];

			for (	MAPAINDIND::iterator mapit = unique_verts_inds_by_geo[i].begin();
						mapit != unique_verts_inds_by_geo[i].end();
						mapit++)
			{
				//	Key and value.
				int isrc = mapit->first;
				int idst = mapit->second;

				//	Vertex coordinate.
				//	Construction of the submeshes.
				Vector3 v3((Real)M0.vertex(isrc)->operator [](0),
						(Real)M0.vertex(isrc)->operator [](1),
						(Real)M0.vertex(isrc)->operator [](2));
				mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition[idst]=v3;

				//	Normal coordinates.
				Vector3 v3n((Real)M0.normal(isrc)(0), (Real)M0.normal(isrc)(1), (Real)M0.normal(isrc)(2));
				mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal[idst]=v3n;

				//	Texture coordinates.
				Vector2 v2((Real)M0.texcoord(isrc)(0), (Real)M0.texcoord(isrc)(1));
				mGeoMesh->mSubMesh[i].mVertexBuffer->mTexCoords[idst]=v2;
			}
		}
		else
		{
			//	Leaves mesh.
			mGeoMesh->mSubMesh[i].mIndexCount=objmesh->mSubMesh[i].mIndexCount;
			mGeoMesh->mSubMesh[i].mIndex=new Index[mGeoMesh->mSubMesh[i].mIndexCount];

			//	Copy the leaves submesh indexes.
			for (unsigned int j=0; j<mGeoMesh->mSubMesh[i].mIndexCount; j++)
			{	
				mGeoMesh->mSubMesh[i].mIndex[j]=objmesh->mSubMesh[i].mIndex[j];
			}

			//	Copy the leaves submesh vertices.
			mGeoMesh->mSubMesh[i].mVertexBuffer=new VertexBuffer();
			mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount=objmesh->mSubMesh[i].mVertexBuffer->mVertexCount;
			mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexInfo=objmesh->mSubMesh[i].mVertexBuffer->mVertexInfo;

			//	Allocate memory for position, normal and texture coordinates.
			mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition=new Vector3[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];
			mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal=new Vector3[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];
			mGeoMesh->mSubMesh[i].mVertexBuffer->mTexCoords=new Vector2[mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount];

			for (	unsigned int j = 0;
					j < mGeoMesh->mSubMesh[i].mVertexBuffer->mVertexCount;
					j++)
			{
				//	Position.
				mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition[j].x=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x;
				mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition[j].y=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y;
				mGeoMesh->mSubMesh[i].mVertexBuffer->mPosition[j].z=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z;

				//	Normals.
				mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal[j].x=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].x;
				mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal[j].y=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].y;
				mGeoMesh->mSubMesh[i].mVertexBuffer->mNormal[j].z=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].z;

				//	Textures.
				mGeoMesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].x=objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].x;
				mGeoMesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].y=objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].y;
			}
		}
	}

	delete[] unique_verts_inds_by_geo;
	delete[] ver_inds_rep_by_geo;
	delete[] inextvert;

	//	Store the simplification steps in MeshSimplificationSequence.
	int acum	=	0;

	for (size_t i	=	0; i < objmesh->mSubMeshCount; i++)
	{
		if (objmesh->mSubMesh[i].mSharedVertexBuffer)
		{
			first_index_submesh[i]	=	0;
		}
		else
		{
			first_index_submesh[i]	=	acum;
		}

		acum	+=	(int)objmesh->mSubMesh[i].mVertexBuffer->mVertexCount;
	}

	MeshSimplificationSequence::Step stepaux, newstep;
	vector<int> v0,v1,submesh0,submesh1;
	MeshSimplificationSequence *msimpseq;

	msimpseq	=	new MeshSimplificationSequence();

	for (unsigned i = 0;	i < decim_data.size();	i++) 
	{
		stepaux		=	decim_data.operator [](i);
		v0				=	vertexbuffermap[stepaux.mV0];  // vertexbuffer index
		submesh0	=	submeshmap[stepaux.mV0]; // displacement by submesh
		v1				=	vertexbuffermap[stepaux.mV1];
		submesh1	=	submeshmap[stepaux.mV1];

		for (unsigned j = 0;	j < v0.size();	j++)
		{
			for (unsigned k = 0;	k < v1.size(); k++)
			{
				//	Obtain the submesh index in VertexBuffer
				//	and the displacement by submesh
				newstep.mV0	=	v0.operator [](j)+first_index_submesh[submesh0.operator [](j)];
				newstep.mV1	=	v1.operator [](k)+first_index_submesh[submesh1.operator [](k)];
				newstep.x		=	stepaux.x;
				newstep.y		=	stepaux.y;
				newstep.z		=	stepaux.z;

				// mT0 y mT1 are unique triangles identifiers
				// returned by the simplification method
				newstep.mT0	=	stepaux.mT0;
				newstep.mT1	=	stepaux.mT1;

				if (j == v0.size() - 1 && k == v1.size() - 1)
				{
					newstep.obligatorio	=	0;
				}
				else
				{
					newstep.obligatorio	=	1;
				}

				if (j == 0 && k == 0)
				{
					//	First simplification step.
					//	unique triangle identifiers 
					//	returned by the simplification method.
					newstep.mModfaces	=	stepaux.mModfaces;
				}
				else 
				{
					std::vector<Index>	vvacio;

					newstep.mModfaces	=	vvacio;
				}

				msimpseq->mSteps.push_back(newstep);
			}
		}
	}

	// WriteOBJ();
	return msimpseq;
}

//---------------------------------------------------------------------------
//	Set submesh leaves.
//---------------------------------------------------------------------------
void SimplificationMethod::setMeshLeaves(int meshLeaves)
{
	indexMeshLeaves	=	meshLeaves;
}

//	Gets mesh simplified.
Mesh	*	SimplificationMethod::GetMesh()
{
	return	mGeoMesh;
}