#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 Geometry::Mesh *m) 
{
	objmesh=m;
	first_index_submesh = new unsigned int[objmesh->mSubMeshCount];
	indexMeshLeaves=-1;
}

//Destructor
SimplificationMethod::~SimplificationMethod()
{
/*	for (std::vector<qslim::pair_info*>::iterator it=pointers_to_remove.begin(); it!=pointers_to_remove.end(); it++)
	{
		if (*it!=NULL)
		{
			delete *it;
			*it=NULL;
		}
	}
*/	delete heap;
	delete [] first_index_submesh;
}

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

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

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

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

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

    if( qslim::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(qslim::Face& F, qslim::Vertex *v1, qslim::Vertex *v2, qslim::Vec3& vnew, qslim::Vec3& f1, qslim::Vec3& f2, qslim::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

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

    M.contractionRegion(v1, v2, changed);

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

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

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

		// Only consider non-degenerate faces
		if( nmapped < 2 )
		{
			qslim::Plane Pnew(f1, f2, f3);
			qslim::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(qslim::Vertex *v0, qslim::Vertex *v1)
{
    const qslim::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
qslim::pair_info *SimplificationMethod::new_pair(qslim::Vertex *v0, qslim::Vertex *v1)
{
    qslim::vert_info& v0_info = vertex_info(v0);
    qslim::vert_info& v1_info = vertex_info(v1);

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

	return pair;
}

//Remove a given pair
void SimplificationMethod::delete_pair(qslim::pair_info *auxpair)
{
    qslim::vert_info& v0_info = vertex_info(auxpair->v0);
    qslim::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(qslim::Model& m, qslim::pair_info *auxpair)
{
    qslim::Vertex *v0 = auxpair->v0;  qslim::Vertex *v1 = auxpair->v1;
    qslim::vert_info& v0_info = vertex_info(v0);
    qslim::vert_info& v1_info = vertex_info(v1);
	qslim::Vec3 vnew = auxpair->candidate;
    int i;

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

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

	// Find the edge to remove (before remove it (contract()))
	qslim::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 qslim::face_buffer changed;
    changed.reset();
    m.contract(v0, v1, vnew, changed);

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

	// Stores the displacement in the simplification step
	simplifstep.x=(float)vdesp[qslim::X];
	simplifstep.y = (float)vdesp[qslim::Y];
	simplifstep.z = (float)vdesp[qslim::Z];
	simplifstep.mV1 = v0->validID(); // the simplification method returns as v0 the alive vertex
	simplifstep.mV0 = v1->validID(); // the simplification method returns as v1 the dead vertex

	// 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++)
	{
		qslim::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++)
    {
		qslim::pair_info *p = v0_info.pairs(i);
		compute_pair_info(p);
    }

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

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

		qslim::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(qslim::Vertex *v, qslim::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(qslim::Model& m)
{
	qslim::heap_node *top;
	qslim::pair_info *pair;

    for(;;)
    {
		top = heap->extract();
		if( !top ) return;
		pair = (qslim::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++)
		{		
			qslim::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
}

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

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

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

    return err;
}

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

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

		if( pair->isValid() )
			break;

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

    return pair->cost;
}

//Returns the maximum error by vertex
qslim::real SimplificationMethod::decimate_max_error(qslim::Model& m)
{
	qslim::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(qslim::Model& m, qslim::real limit)
{
    int i,j;

    vinfo.init(m.vertCount());

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

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

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

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

				B *= qslim::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 qslim::Heap(m.validEdgeCount);

    int pair_count = 0;

    for(i=0; i<m.edgeCount(); i++)
		if( m.edge(i)->isValid() )
		{
			qslim::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 )
    {
		qslim::ProxGrid grid(m.bounds.min, m.bounds.max, limit);
		for(i=0; i<m.vertCount(); i++)
			grid.addPoint(m.vertex(i));

		qslim::buffer<qslim::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++)
			{
				qslim::Vertex *v1 = m.vertex(i);
				qslim::Vertex *v2 = (qslim::Vertex *)nearby(j);

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

			}
		}
    }
}

//Initilizations
void SimplificationMethod::simplfmethod_init(void)
{
    int i;

	// Convertir la estructura geomesh a modelo qslim
	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::simplfmethod_run(int finalfaces,Geometry::TIPOFUNC upb)
{
	decimate_init(M0, qslim::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<Geometry::MeshSimplificationSequence::Step>::iterator it= decim_data.begin(); it!=decim_data.end(); it++)
	{
		it->mModfaces.clear();
	}
	decim_data.clear();

	while( M0.validFaceCount > finalfaces/*qslim::face_target*/
		&& M0.validFaceCount > 1
		&& decimate_min_error() < qslim::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::simplfmethod_runv(int numvertices,Geometry::TIPOFUNC upb)
{
	decimate_init(M0, qslim::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<Geometry::MeshSimplificationSequence::Step>::iterator it= decim_data.begin(); it!=decim_data.end(); it++)
	{
		it->mModfaces.clear();
	}
	decim_data.clear();

	while(	M0.validVertCount > numvertices	/*qslim::face_target*/
					&&
					decimate_min_error() < qslim::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 not produce 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(qslim::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(qslim::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(qslim::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(qslim::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(qslim::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(qslim::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\n");
			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)
{
	Geometry::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
Geometry::MeshSimplificationSequence *SimplificationMethod::Decimate(float lod, Geometry::Mesh **mesh, int simpliftype, Geometry::TIPOFUNC upb)
{
	simplfmethod_init();

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


	typedef std::map<int,int> MAPAINDIND;   //Store unique vertices by submesh (without repetitions)
	typedef std::vector<int> REPINDLIST; //Store all the indices by submesh (with repetitions)

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

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

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

	qslim::Vertex *auxvertex;

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

			//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_repes_x_geo[igeo].push_back(auxface->vertex(0)->validID());
			ver_inds_repes_x_geo[igeo].push_back(auxface->vertex(1)->validID());
			ver_inds_repes_x_geo[igeo].push_back(auxface->vertex(2)->validID());
		}
	}


	//Output simplified mesh

	*mesh=new Geometry::Mesh();
	(*mesh)->mVertexBuffer	=	new	Geometry::VertexBuffer();
	(*mesh)->mType=objmesh->mType;

	(*mesh)->mSubMeshCount = objmesh->mSubMeshCount;
	//Memory allocation for the submeshes
	(*mesh)->mSubMesh = new Geometry::SubMesh[objmesh->mSubMeshCount];

	bool copiedShared = false;
	for (size_t i=0; i<objmesh->mSubMeshCount; i++) 
	{
		(*mesh)->mSubMesh[i].mSharedVertexBuffer=false; //objmesh->mSubMesh[i].mSharedVertexBuffer;
		(*mesh)->mSubMesh[i].mStripCount=0;
		(*mesh)->mSubMesh[i].mStrip=NULL;

		if (i!=indexMeshLeaves)
		{
			//Indices vectors
			(*mesh)->mSubMesh[i].mIndexCount=ver_inds_repes_x_geo[i].size(); 
			(*mesh)->mSubMesh[i].mIndex=new Geometry::Index[(*mesh)->mSubMesh[i].mIndexCount];
			//Store the indices
			for (unsigned int j=0; j<(*mesh)->mSubMesh[i].mIndexCount; j++) 
			{	
				//Obtain the indices that point at VertexBuffer
				(*mesh)->mSubMesh[i].mIndex[j]=unique_verts_inds_by_geo[i].operator [](ver_inds_repes_x_geo[i].operator [](j));
			}
			

			//Copy the vertices
			if ((*mesh)->mSubMesh[i].mSharedVertexBuffer) 
			{
				//Shared vertices
				if (!copiedShared)
				{
					(*mesh)->mVertexBuffer = new Geometry::VertexBuffer();
					copiedShared = true;
				}
				(*mesh)->mSubMesh[i].mVertexBuffer = (*mesh)->mVertexBuffer;
			}
			else 
			{
				//There's no shared vertices
				(*mesh)->mSubMesh[i].mVertexBuffer=new Geometry::VertexBuffer();
			}

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

			//Allocate memory for position, normal and texutre coordinates
			(*mesh)->mSubMesh[i].mVertexBuffer->mPosition=new Geometry::Vector3[(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount];
			(*mesh)->mSubMesh[i].mVertexBuffer->mNormal=new Geometry::Vector3[(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount];
			(*mesh)->mSubMesh[i].mVertexBuffer->mTexCoords=new Geometry::Vector2[(*mesh)->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
				Geometry::Vector3 v3((Geometry::Real)M0.vertex(isrc)->operator [](0),	
											(Geometry::Real)M0.vertex(isrc)->operator [](1),
											(Geometry::Real)M0.vertex(isrc)->operator [](2));
				(*mesh)->mSubMesh[i].mVertexBuffer->mPosition[idst]=v3;

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

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

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

			//copy the leaves submesh vertices
			(*mesh)->mSubMesh[i].mVertexBuffer=new Geometry::VertexBuffer();
			(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount=objmesh->mSubMesh[i].mVertexBuffer->mVertexCount;
			(*mesh)->mSubMesh[i].mVertexBuffer->mVertexInfo=objmesh->mSubMesh[i].mVertexBuffer->mVertexInfo;

			//allocate memory for position, normal and texture coordinates
			(*mesh)->mSubMesh[i].mVertexBuffer->mPosition=new Geometry::Vector3[(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount];
			(*mesh)->mSubMesh[i].mVertexBuffer->mNormal=new Geometry::Vector3[(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount];
			(*mesh)->mSubMesh[i].mVertexBuffer->mTexCoords=new Geometry::Vector2[(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount];
			for (unsigned int j=0; j<(*mesh)->mSubMesh[i].mVertexBuffer->mVertexCount; j++) 
			{	
				//position
				(*mesh)->mSubMesh[i].mVertexBuffer->mPosition[j].x=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].x;
				(*mesh)->mSubMesh[i].mVertexBuffer->mPosition[j].y=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].y;
				(*mesh)->mSubMesh[i].mVertexBuffer->mPosition[j].z=objmesh->mSubMesh[i].mVertexBuffer->mPosition[j].z;

				//normals
				(*mesh)->mSubMesh[i].mVertexBuffer->mNormal[j].x=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].x;
				(*mesh)->mSubMesh[i].mVertexBuffer->mNormal[j].y=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].y;
				(*mesh)->mSubMesh[i].mVertexBuffer->mNormal[j].z=objmesh->mSubMesh[i].mVertexBuffer->mNormal[j].z;

				//textures
				(*mesh)->mSubMesh[i].mVertexBuffer->mTexCoords[j].x=objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].x;
				(*mesh)->mSubMesh[i].mVertexBuffer->mTexCoords[j].y=objmesh->mSubMesh[i].mVertexBuffer->mTexCoords[j].y;
			}
		}
	}
	delete[] unique_verts_inds_by_geo;
	delete[] ver_inds_repes_x_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;
	}
	Geometry::MeshSimplificationSequence::Step stepaux, newstep;
	vector<int> v0,v1,submesh0,submesh1;
	Geometry::MeshSimplificationSequence *msimpseq;
	msimpseq=new Geometry::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;
				newstep.mT0=stepaux.mT0; // mT0 y mT1 are unique triangles identifiers returned by the simplification method
				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
                    newstep.mModfaces=stepaux.mModfaces; //unique triangle identifiers returned by the simplification method
				else 
				{
					std::vector<Geometry::Index> vvacio;
					newstep.mModfaces=vvacio;
				}
					
				msimpseq->mSteps.push_back(newstep);
			}
		}
	}

	// WriteOBJ();
	return msimpseq;
}

void SimplificationMethod::setMeshLeaves(int meshLeaves)
{
	indexMeshLeaves=meshLeaves;
}