#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <float.h>
#include <math.h>

#include "../include/metrics.h"
#include "../include/global.h"
#include "../include/simplify.h"
#include "../include/saliency.h"
#include "../include/histogram.h"

//#define RECOMPUTE_THE_WHOLE_HEAP

using namespace	VMI;
using	namespace	std;

//	Multimap to store current vertex in a simplification state.
typedef	multimap<_float3_,	int>	t_map;
typedef	t_map::iterator						t_map_str;
typedef	t_map::value_type					t_pair;
t_map															mVertexMap;
multimap<int,int>									mVertices;

vector<Geometry::Vector3>	VMI::vPositions;
vector<Geometry::Vector3>	VMI::vNormals;
vector<Geometry::Vector2>	VMI::vTexCoords;

///////////////////////////////////////////////////////////////////////////////

void VMI::bh_mydump(Mesh *mesh, bheap_t *h) {
	int i, d;

	printf("Heap status.\n");
	for(i = 1; i <= h->n; i++) {
		d = h->a[i].item;
		printf ("Heap [%d] = pri %f, e%d(%d, %d) d:%d\n",
				i,
				h->a[i].key,d,
				mesh->edges[d].u,
				mesh->edges[d].v, 
				h->a[i].dirty);
	}
	printf("\n");

	for(i = 2; i <= h->n; i++) {
		if(h->a[i].key < h->a[i/2].key) {
			printf("key error at entry %d, value %f\n", i, h->a[i].key);
			exit(1);
		}
	}
	for(i = 1; i <= h->n; i++) {
		if(h->p[h->a[i].item] != i) {
			printf("indexing error at entry %d", i);
			exit(1);
		}
	}    
}

bheap_t *VMI::initHeap(Mesh *mesh) {
	GLuint i;
	double cost;
	bheap_t *h = NULL;

	printf("Creating a heap for simplification...");

	h = bh_alloc(mesh->numEdges);

	//	Init percent update count.
	float	percent	=	60.0 / mesh->numEdges;

	for (i = 0; i <mesh->numEdges; i++) {

		//	Update progress bar.
		mUPB(percent);

		if (mesh->edges[i].enable == TRUE) {

			cost = computeEdgeCost(mesh, i);

			// Add only valid edges
			if (cost != FLT_MAX) {
				bh_insert(h, i, cost);
				/* set dirty flag to FALSE */
				bh_mark(h, i, FALSE);
			}
		}
	}

	printf("Ok\n");

#ifdef DEBUG 
	bh_mydump(mesh, h);
	//getchar();
#endif

	return h;
}

bheap_t *VMI::updateHeap(bheap_t *h, Mesh *mesh, Change *c) {
	int e, i;
#ifdef RECOMPUTE_THE_WHOLE_HEAP
	bheap_t *nh;
	double cost;
#else
	int *lv, *le, numV, numE;
#endif

#ifdef RECOMPUTE_THE_WHOLE_HEAP

	nh = bh_alloc(mesh->numEdges);

	// Recompute heap
	for(i = 1; i <= h->n; i++) {
		e = h->a[i].item;

		cost = computeEdgeCost(mesh, e);

		// Add only valid edges
		if (cost != FLT_MAX)
			bh_insert(nh, e, cost);
	}

	bh_free(h);

	h = nh;

#else

	// Compute vertices adjacent to a vertex
	lv = verticesAdjToVertex(mesh, c->v, &numV);

	// Compute edges adjacent to vertices
	le = edgesAdjToVertices(mesh, lv, numV, &numE);

	free(lv);

	// Set a dirty flag to adjacent edges from the heap
	for(i = 0; i <numE; i++) {
		e = le[i];

		bh_mark(h, e, TRUE);
	}

	free(le);
#endif

#ifdef DEBUG 
	bh_mydump(mesh, h);
	//getchar();
#endif

	return h;
}

GLdouble VMI::computeEdgeLength(Vertex *vertices, int u, int v) {
	GLfloat ux, uy, uz, vx, vy ,vz, rx, ry, rz;

	ux = vertices[u].x;
	uy = vertices[u].y;
	uz = vertices[u].z;

	vx = vertices[v].x;
	vy = vertices[v].y;
	vz = vertices[v].z;

	rx = ux - vx;
	ry = uy - vy;
	rz = uz - vz;

	return sqrt((rx * rx) + (ry * ry) + (rz * rz));
}

///////////////////////////////////////////////////////////////////////////////

void VMI::chooseBestEndPoints(Mesh *mesh, int e) {
	int v = mesh->edges[e].v;
	int u = mesh->edges[e].u;
	int numTu = mesh->vertices[u].numTriangles;
	int numTuv = 0;
	int *Tuv = (int *)malloc(numTu * sizeof(int));
	int i, j, f, n, n0, n1, n2;
	float cost_u_v, cost_v_u, mincurve, curve, dot;

	// Compute Tuv
	for (i=0; i<(int)numTu; i++) {

		n = mesh->vertices[u].triangles[i];

		n0 = mesh->triangles[n].indices[0];
		n1 = mesh->triangles[n].indices[1];
		n2 = mesh->triangles[n].indices[2];

		if ((n0 == v) || (n1 == v) || (n2 == v))

			Tuv[numTuv++] = n;
	}
	//printItemList(Tuv, numTuv);

	curve = 0.0f;
	for (i=0; i<numTu; i++) {
		n = mesh->vertices[u].triangles[i];
		mincurve = 1.0f;
		for (j=0; j<numTuv; j++) {
			f = Tuv[j];
			dot = glmDot(mesh->triangles[f].normal, mesh->triangles[n].normal);
			//printf("dot: %f \n", dot);
			mincurve = MIN(mincurve, (1.0f - dot) / 2.0f);
		}
		//printf("mincurve: %f \n", mincurve);
		curve = MAX(curve, mincurve);
	}
	cost_u_v = curve;

	curve = 0.0f;
	for (i=0; i<(int)mesh->vertices[v].numTriangles; i++) {
		n = mesh->vertices[v].triangles[i];
		mincurve = 1.0f;
		for (j=0; j<numTuv; j++) {
			f = Tuv[j];
			dot = glmDot(mesh->triangles[f].normal, mesh->triangles[n].normal);
			//printf("%f \n", dot);
			mincurve = MIN(mincurve, (1.0f - dot) / 2.0f);
		}
		curve = MAX(curve, mincurve);
	}

	cost_v_u = curve;

	//printf("e%d(%d,%d) cost: %f \n", e, u, v, cost_u_v);
	//printf("e%d(%d,%d) cost: %f \n", e, v, u, cost_v_u);

	if ( (cost_v_u + 0.000001) < cost_u_v) {
		//printf("Swap edge\n");
		//printf("e%d(%d,%d) cost: %f \n", e, u, v, cost_u_v);
		//printf("e%d(%d,%d) cost: %f \n", e, v, u, cost_v_u);
		
		swap((int *)&mesh->edges[e].u, (int *)&mesh->edges[e].v);
	}

	free(Tuv);
}
///////////////////////////////////////////////////////////////////////////////

void saveProjectedAreas(GLuint **histogram, GLuint numCameras, Change *c, GLuint *dest) {
	GLuint i, j, t, n = 0;

	for (j=0; j<(GLuint)c->numDel; j++) {

		t = c->deleted[j].id;

		for (i=0; i<numCameras; i++) {

			dest[n] = histogram[i][t + 1];
			n++;
		}
	}

	for (j=0; j<(GLuint)c->numMod; j++) { 

		t = c->modified[j].id;

		for (i=0; i<numCameras; i++) {

			dest[n] = histogram[i][t + 1];
			n++;
		}   
	}

	// Save background area
	for (i=0; i<numCameras; i++) {

		dest[n] = histogram[i][0];
		n++;
	}
}

void loadProjectedAreas(GLuint *src, GLuint numCameras, Change *c, GLuint **histogram) {
	GLuint i, j, t, n = 0;

	for (j=0; j<(GLuint)c->numDel; j++) { 

		t = c->deleted[j].id;

		for (i=0; i<numCameras; i++) {

			histogram[i][t + 1] = src[n];
			n++;
		}
	}

	for (j=0; j<(GLuint)c->numMod; j++) { 

		t = c->modified[j].id;

		for (i=0; i<numCameras; i++) {

			histogram[i][t + 1] = src[n];
			n++;
		}
	}

	// Load background area
	for (i=0; i<numCameras; i++) {

		histogram[i][0] = src[n];
		n++;
	}
}

///////////////////////////////////////////////////////////////////////////////

GLdouble VMI::computeEdgeCost(Mesh *mesh, int e) {
	GLdouble cost = 0.0, sal = 1.0;
#ifdef MI
	GLdouble *auxIs = NULL;
#else
	GLdouble newI;
#endif
	GLuint i, *histoAux = NULL;
	Change *c;

	chooseBestEndPoints(mesh, e);
	//getchar();

	c = createChange(mesh, e);

	// Compute cost only for boundary or manifold edges
	if (c->numDel == 2
			/*(c->numDel <3) && (c->numDel > 0)*/) {
		// Allocate memory
		histoAux = (GLuint *)malloc((c->numDel + c->numMod + 1) * numCameras * sizeof(GLuint));

#ifdef MI // Mutual Information
		///////////////////////////////////////////////////////////////////////////////
		// Allocate memory
		auxIs = (GLdouble *)malloc(numCameras * sizeof(GLdouble));

		//printFullHistogram(histogram, numCameras, mesh->numTriangles);

		for (i=0; i<numCameras; i++) {
			// Copy old informations
			auxIs[i] = initialIs[i];

			initialIs[i] = decMI(initialIs[i], histogram, numCameras, i, c);
		}

		saveProjectedAreas(histogram, numCameras, c, histoAux);
    // Apply the edge collapse
		doChange(mesh, c);

		// Compute the cost of an edge collapse
		getProjectedAreasWin(histogram, numCameras, c);
		//printFullHistogram(histogram, numCameras, mesh->numTriangles);
		//getchar();

#ifdef SALIENCY
		sal = computeEdgeSaliency(mesh, c, percentile);
		//printf("  e:%d s: %f\n", c->e, sal);
#endif

		for (i=0; i<numCameras; i++) {

			initialIs[i] = incMI(initialIs[i], histogram, numCameras, i, c);
			//printf("  I0:%f Is: %f\n", auxE[i], initialIs[i]);

			cost += (ABS(auxIs[i] - initialIs[i]) * cameras[i].weight * sal);

			// Restore old informations
			initialIs[i] = auxIs[i];
		}

		undoChange(mesh, c);

		loadProjectedAreas(histoAux, numCameras, c, histogram);
		//printFullHistogram(histogram, numCameras, mesh->numTriangles);
		//exit(1);

		// Free memory
		free(auxIs);
#else
		///////////////////////////////////////////////////////////////////////////////

		saveProjectedAreas(histogram, numCameras, c, histoAux);

		doChange(mesh, c);

		// Compute the cost of an edge collapse
		getProjectedAreasWin(histogram, numCameras, c);
		//printFullHistogram(histogram, numCameras, mesh->numTriangles);

#ifdef SALIENCY
		sal = computeEdgeSaliency(mesh, c, percentile);
		//printf("  e:%d s: %f\n", c->e, sal);
#endif
		for (i=0;i <numCameras; i++) {

#ifdef KL // Kullback-Leibler
			newI = computeKL(mesh, histogram[i]);
#endif
#ifdef HE // Hellinger
			newI = computeHE(mesh, histogram[i]);
#endif
#ifdef CS // Chi-Square
			newI = computeCS(mesh, histogram[i]);
#endif
			cost += (ABS(initialIs[i] - newI) * cameras[i].weight * sal);
		}
		undoChange(mesh, c);

		loadProjectedAreas(histoAux, numCameras, c, histogram);
		//printFullHistogram(histogram, numCameras, mesh->numTriangles);
		//exit(1);

		///////////////////////////////////////////////////////////////////////////////
#endif
		// Free memory
		free(histoAux);

	} else cost = FLT_MAX;

	deleteChange(c);

	return cost;
}
/* Get a valid edge from the heap */
int getMinValidEdge(Mesh *mesh, bheap_t *h) {
	int e;

	/* pick the minimum-cost valid edge from heap */
	e = bh_min(h);

	while (mesh->edges[e].enable == FALSE) {
		/* delete invalid edge */
		bh_delete(h, e);
		/* pick again*/
		e = bh_min(h);
	}

	return e;
}

/* Get the mininum edge cost from the heap using lazy evaluation */
int getMinEdge(Mesh *mesh, bheap_t *h) {
	double cost;
	int e;

	/* pick the minimum-cost edge from heap */
	e = getMinValidEdge(mesh, h);

	while (bh_is_dirty(h, e)) {
		/* delete edge from heap */
		bh_delete(h, e);
		mesh->edges[e].enable = FALSE;
		/* recompute cost*/
		cost = computeEdgeCost(mesh, e);

		// Add only valid edges
		if (cost != FLT_MAX) {
			/* reinsert into the heap */
			bh_insert(h, e, cost);
			mesh->edges[e].enable = TRUE;
			/* set dirty flag to FALSE */
			bh_mark(h, e, FALSE);
		}

		/* pick again */
		e = getMinValidEdge(mesh, h);
	}

	return e;
}

///////////////////////////////////////////////////////////////////////////
void VMI::simplifyModel(Mesh *mesh, GLuint numDemandedTri)
{
	int			e;
	Change	*c;
	FILE		*file 					= NULL;
	FILE		*file_simp_seq	=	NULL;
	char		s[MAX_CHAR];
	double	cost 						= 0.0;
	bheap_t	*h 							= NULL;
	float		percent	=	(40.0) / (float)(mesh->numTriangles - numDemandedTri);

	if (numDemandedTri > mesh->currentNumTriangles)
	{
		return;
	}

	// Save changes into a file
	if (bSaveLog == TRUE)
	{
#ifdef SALIENCY       
		sprintf(s,"%s_%s_S.log", filename, EXT);
#else
		sprintf(s,"%s_%s.log", filename, EXT);
#endif

		glmWriteOBJ(pmodel, s, GLM_NONE);

		/* open the file */
		file = fopen(s, "a+");

		if (!file)
		{
			fprintf(stderr, "simplifyModel() failed: can't open file \"%s\" to write.\n", s);
			exit(1);
		}
	}

	//	Open file of simplification sequence.
	file_simp_seq	=	fopen("SimplifSequence.txt", "w");

	if (!file_simp_seq)
	{
		fprintf(stderr,
				"simplifyModel() failed: ",
				"can't open file \"SimplifSequence\" to write.\n");
	}

	h = initHeap(mesh);

	printf("Starting simplification...\n");

	while (mesh->currentNumTriangles > numDemandedTri /*&& cost == 0.0*/)
	{
		// Get the edge that has the minimum cost
		//e = getMinEdge(mesh, h);
		e = extractValidEdge(mesh, h);

		c = createChange(mesh, e);

		// Apply the edge collapse
		doChange(mesh, c);
		//deleteVertexOfMap(mesh, c->u);

		// Write Simplification sequence.
		saveSimplificationSequence(c,0);

		//contractTwinVertices(mesh, c->u, c->v); // vertex u is removed

		if (bSaveLog == TRUE) writeChange(file, c);

		cost = h->a[h->p[e]].key;

		/*
		printf(	"Edge collapse e%d(%d,%d) %f MIN d %d m %d\n",
				c->e,
				c->u,
				c->v,
				cost,
				c->numDel,
				c->numMod);
		*/

		mesh->edges[e].enable = FALSE;
		bh_delete(h, e);

		// Get projected areas after the edge collapse
		getProjectedAreas(histogram, numCameras);

		computeCameraIs(histogram, numCameras, initialIs);

		deleteEdges(mesh, c);
		modifyEdges(mesh, c);

		// Update the heap according to the edge collapse
		h = updateHeap(h, mesh, c);

		mUPB((float)c->numDel * percent);

		deleteChange(c);

		//printf("t %d\n", mesh->currentNumTriangles);
	}

	//	Debug.
	cout	<<	"Number of vertices: "
		<<	mesh->currentNumVertices
		<<	endl;

	if (bSaveLog == TRUE)
	{
		fclose(file);
		printf("Log file written...Ok\n");
	}

	//	Close simplification sequence file.
	fclose(file_simp_seq);

	bh_free(h);
}

///////////////////////////////////////////////////////////////////////////////
//	Extract a correct edge of the heap.
///////////////////////////////////////////////////////////////////////////////
int VMI::extractValidEdge(Mesh *mesh,	bheap_t *h)
{
	int	e;

	while((e	=	isValidEdge(mesh,getMinEdge(mesh, h))) == -1);

	return	e;
}

///////////////////////////////////////////////////////////////////////////////
//	Indicates if an edge is valid.
///////////////////////////////////////////////////////////////////////////////
int	VMI::isValidEdge(Mesh	*mesh,	int edge)
{
	int	result;
	Vertex	*u;
	Vertex	*v;
	_float3_	fu;
	_float3_	fv;

	u	=	&mesh->vertices[mesh->edges[edge].u];
	v	=	&mesh->vertices[mesh->edges[edge].v];

	fu	=	_float3_(u->x,u->y,u->z);
	fv	=	_float3_(v->x,v->y,v->z);

	if ((mVertexMap.find(fu) != mVertexMap.end())
			||
			(mVertexMap.find(fv) != mVertexMap.end()))
	{
		result	= edge;
	}
	//	Edge is not valid.
	else
	{
		result	=	-1;
	}

	return	result;
}

//-------------------------------------------------------------------------
//	Inits the multimap of vertices.
//-------------------------------------------------------------------------
void	VMI::initVertexMultimap(Mesh *mesh,multimap<int,int> &vertexMultimap)
{
	Vertex	*vertex;
	float		x,y,z;

	//	Clears multimap of vertices.
	mVertexMap.clear();
	mVertices.clear();

	mVertices	=	vertexMultimap;

	//	Debug.
	cout	<<	"Vertex Map Elements: "
				<<	mVertices.size()
				<<	endl;

	//	For each vertex.
	for (int	i	=	0;	i < mesh->numVertices;	i++)
	{
		vertex	=	&mesh->vertices[i];

		x	=	vertex->x;
		y	=	vertex->y;
		z	=	vertex->z;

		mVertexMap.insert(t_pair(_float3_(x,y,z),i));
	}
}

//-------------------------------------------------------------------------
//	Deletes a vertex in the multimap.
//-------------------------------------------------------------------------
void	VMI::deleteVertexOfMap(Mesh	*mesh, int u)
{
	_float3_	removed_vert;
	t_map_str	lb;
	t_map_str	ub;

	//	Position of the vertex removed.
	removed_vert	=	_float3_(	mesh->vertices[u].x,
														mesh->vertices[u].y,
														mesh->vertices[u].z);

	//	If position of the removed vertex is found.
	if (mVertexMap.end() != mVertexMap.find(removed_vert))
	{
		lb	=	mVertexMap.lower_bound(removed_vert);
		ub	=	mVertexMap.upper_bound(removed_vert);

		//	For each vertex.
		while	(lb != ub)
		{
			//	If removed vertex is found.
			if ((*lb).second == u)
			{
				//	Debug.
				cout	<<	"Vertex erased: "
							<<	(*lb).second
							<<	endl;

				//	Delete vertex that disappears.
				mVertexMap.erase(lb);
				
				//	Break while.
				lb	=	ub;
			}
			else
			{
				//	Next iteration.
				lb++;
			}
		}
	}
}

//-------------------------------------------------------------------------
//	Compare the coordinates of two vertices.
//-------------------------------------------------------------------------
bool	VMI::compareVertices(Vertex *vertices,int	u,	int v)
{
	if ((vertices[u].x == vertices[v].x)
			&&
			(vertices[u].y == vertices[v].y)
			&&
			(vertices[u].z == vertices[v].z))
	{
		return	true;
	}
	else
	{
		return	false;
	}
}

//-------------------------------------------------------------------------
//	Find edge of the simplification sequence.
//-------------------------------------------------------------------------
void	VMI::contractInitialMesh(Mesh	*mesh)
{
	bool		edge_found;
	Geometry::MeshSimplificationSequence::Step	step;

	multimap<int,int>::iterator	it0;
	multimap<int,int>::iterator	lb0;
	multimap<int,int>::iterator	ub0;
	multimap<int,int>::iterator	lb1;
	multimap<int,int>::iterator	ub1;

	std::vector<Geometry::MeshSimplificationSequence::Step> steps;

	Edge		*econ;
	float		x,y,z;
	int			*lv1;
	int			v0;
	int			v1;
	int			edge;
	int			num_edges;
	Change	*c;

	//	Copy simplification steps of the joined mesh.
	steps	=	mSequence->mSteps;

	mSequence->mSteps.clear();

	//	Debug.
	cout	<<	"Simplification steps: "
				<<	steps.size()
				<<	" new vertices: "
				<<	mSequence->mNewVertices.size()
				<<	endl;

	for	(size_t	i	=	0;	i < steps.size(); i++)
	{
		step	=	steps[i];

		lb0	=	mVertices.lower_bound(steps[i].mV0);
		ub0	=	mVertices.upper_bound(steps[i].mV0);
		lb1	=	mVertices.lower_bound(steps[i].mV1);
		ub1	=	mVertices.upper_bound(steps[i].mV1);

		edge_found	=	false;

		//	Debug.
		cout	<<	"step "
					<<	i
					<<	" V0: "
					<<	steps[i].mV0
					<<	" V1: "
					<<	steps[i].mV1
					<<	endl;

		//	Removed vertex.
		while ((lb1 != ub1) && !edge_found)
		{
			//	Real index.
			v1	=	(*lb1).second;

			lv1	=	edgesAdjToVertex(mesh,v1,&num_edges);

			//	Edje iterator.
			edge	=	0;

			while	((edge < num_edges) &&	!edge_found)
			{
				econ	=	&mesh->edges[lv1[edge]];

				//	Begin of iteration v0.
				it0	=	lb0;

				//	Remaining vertex.
				while ((it0 != ub0) &&	!edge_found)
				{
					//	Real index.
					v0	=	(*it0).second;
					
					if (compareVertices(mesh->vertices,econ->v,v0))
					{
						c = newChange(mesh, econ->u, econ->v);

						edge_found	=	true;
					}
					else if (compareVertices(mesh->vertices,econ->u,v0))
					{
						c = newChange(mesh, econ->v, econ->u);

						edge_found	=	true;
					}
					else
					{
						it0++;
					}
				}

				edge++;
			}

			lb1++;
		}

		if (edge_found)
		{
			//	Debug.
			cout	<<	"Contracting edge of the initial mesh..."	
						<<	endl
						<<	"u: "
						<<	c->u
						<<	" v: "
						<<	c->v
						<<	endl;
	
			//	Collapse new edge.
			doChange(mesh, c); // the mesh has been updated.

			//	Write Simplification sequence.
			saveSimplificationSequence(c,0);

			deleteVertexOfMap(mesh, c->u);
			deleteEdges(mesh, c);
			modifyEdges(mesh, c);
			deleteChange(c);

			//	Contract twin vertices.
			contractTwinVertices(mesh,v1,v0);
		}
	}
}

//-------------------------------------------------------------------------
//	Find twin vertices and contract them.
//-------------------------------------------------------------------------
void VMI::contractTwinVertices(	Mesh	*mesh,
																int		u,
																int		v)
{
	bool					twin_found;
	int						edge;
	int						lonely_vert;
	int						new_vert;
	t_map_str			lb;
	t_map_str			ub;
	t_map_str			it;
	_float3_			fu;
	_float3_			fv;
	Edge					*econ;
	float					x,y,z;
	int						*le;
	int						num_edges;
	Change				*c;

	Geometry::GeoVertex	vertex_added;

	if (!compareVertices(mesh->vertices,u,v))
	{
		//take_bone_from_vert = v;

		fu	=	_float3_(	mesh->vertices[u].x,
										mesh->vertices[u].y,
										mesh->vertices[u].z);

		fv	=	_float3_(	mesh->vertices[v].x,
										mesh->vertices[v].y,
										mesh->vertices[v].z);

		//	Find vertices width u coordinates.
		while ((it = mVertexMap.find(fu)) != mVertexMap.end())
		{
			twin_found	=	false;

			lonely_vert	=	(*it).second;

			le	=	edgesAdjToVertex(mesh,lonely_vert,&num_edges);

			//	Find an edge width v coordinates.
			int	i	=	0;

			while((i < num_edges) && !twin_found)
			{
				econ	=	&mesh->edges[le[i]];

				if (compareVertices(mesh->vertices,econ->u,v)
						||
						compareVertices(mesh->vertices,econ->v,v))
				{
					lb	=	mVertexMap.lower_bound(fv);
					ub	=	mVertexMap.upper_bound(fv);

					//	For each vertex.
					while	(lb != ub)
					{
						//	If removed vertex is found.
						if (((*lb).second == econ->u)
								||
								((*lb).second == econ->v))
						{
							twin_found	=	true;

							//	Break while.
							lb	=	ub;
						}
						else
						{
							//	Next iteration.
							lb++;
						}
					}
				}
				i++;
			}

			//	If a twin edge has been found.
			if (twin_found)
			{
				//	Debug.
				cout	<<	"Twin Collapsed..."	<<	endl;

				//	Compare vertices coordinates.
				if (compareVertices(mesh->vertices,econ->u,v))
				{
					// New edge for the change
					c = newChange(mesh, econ->v, econ->u);

					//	Debug.
					cout	<<	"--Reverse--"	<<	endl;
				}
				else
				{
					// New edge for the change
					c = newChange (mesh, econ->u, econ->v);
				}

				// Collapse new edge.
				doChange(mesh, c); // the mesh has been updated.

				// Write Simplification sequence.
				saveSimplificationSequence(c,1);

				deleteVertexOfMap(mesh, c->u);
				deleteEdges(mesh, c);
				modifyEdges(mesh, c);
				deleteChange(c);
			}
			else
			{
				//	Debug.
				cout	<<	"Collapsing new edge..."	<<	endl;

				x	=	mesh->vertices[v].x;
				y	=	mesh->vertices[v].y;
				z	=	mesh->vertices[v].z;

				mesh->vertices = addVertex(	mesh->vertices,
																		(int *)&mesh->numVertices,
																		x,
																		y,
																		z,
																		&new_vert);
				
				// When a new vertex is added to the mesh, not only the
				// total number of vertices is increased but also current number
				mesh->currentNumVertices++;

				//	Adds new vertex to multimap.
				mVertexMap.insert(t_pair(_float3_(x,y,z),new_vert));

				//	Creates new edge.
				mesh->edges	=	addEdge(mesh->edges,
															(int *)&mesh->numEdges,
															lonely_vert,
															new_vert,
															&edge);

				/*
				//	Debug.
				cout	<<	"lonely_vert"
					<<	lonely_vert
					<<	"("
					<<	mesh->vertices[lonely_vert].x
					<<	","
					<<	mesh->vertices[lonely_vert].y
					<<	","
					<<	mesh->vertices[lonely_vert].z
					<<	")"
					<<	endl;
				//	Debug.
				cout	<<	"new_vert"
					<<	new_vert
					<<	"("
					<<	mesh->vertices[new_vert].x
					<<	","
					<<	mesh->vertices[new_vert].y
					<<	","
					<<	mesh->vertices[new_vert].z
					<<	")"
					<<	endl;
				*/

				//	We assume here there are the same number of vertices
				//	and texture coordinates and normals.

				//	Assigns the position of the vertex.
				vPositions.push_back(Geometry::Vector3(x,y,z));

				//	Assigns a texture coordinate for the vertex.
				vTexCoords.push_back(vTexCoords[lonely_vert]);

				//	Assigns a normal coordinate for the vertex.
				vNormals.push_back(vNormals[lonely_vert]);

				//	Adds new vertex information to simplification sequence.
				vertex_added.id				=	new_vert;
				vertex_added.bonefrom	=	v;

				vertex_added.position	=	vPositions[new_vert];

				vertex_added.texcoord	=	vTexCoords[new_vert];

				vertex_added.normal	=	vNormals[new_vert];

				mSequence->mNewVertices.push_back(vertex_added);

				// Collapse new edge
				c = newChange(mesh, lonely_vert, new_vert);

				doChange(mesh, c); // the mesh has been updated

				// Write Simplification sequence.
				saveSimplificationSequence(c,1);

				deleteVertexOfMap(mesh, c->u);
				deleteEdges(mesh, c);
				modifyEdges(mesh, c);
				deleteChange(c);
			}
		}
	}
}