#include "Mesh.h"
#include <set>
#include <deque>

//=======================================================================
//=======================================================================

UINT8 *Mesh::m_fileptr = NULL;
UINT8 *Mesh::m_nextchunkptr = NULL;
UINT8 *Mesh::m_nextsubchunkptr = NULL;
UINT8 *Mesh::m_fileend = NULL;
UINT8 *Mesh::m_readptr = NULL;

//=======================================================================
//=======================================================================

Mesh::Mesh()
{
	m_ambiguousNeighbours = 0;
	m_adjacencyCalculated = false;
	m_filename = "";
}

//=======================================================================
//=======================================================================

Mesh::Mesh( int vertices, int triangles )
{
	SS_ASSERT(vertices > 0);
	SS_ASSERT(triangles > 0);

	m_ambiguousNeighbours = 0;
	m_adjacencyCalculated = false;
	m_filename = "";
	m_vertices.resize(vertices);
	m_triangles.resize(triangles);
}

//=======================================================================
//=======================================================================

Mesh::Mesh( const char *fileName )
{
	load( fileName );
}

//=======================================================================
//=======================================================================

Mesh::Mesh( const Mesh &m )
{
	m_vertices = m.m_vertices;
	m_triangles = m.m_triangles;
	m_edges = m.m_edges;
	m_surfaces = m.m_surfaces;
	m_filename = m.m_filename;

	m_boundingBoxMin = m.m_boundingBoxMin;
	m_boundingBoxMax = m.m_boundingBoxMax;
	m_boundingRadius = m.m_boundingRadius;
	m_ambiguousNeighbours = m.m_ambiguousNeighbours;
	m_adjacencyCalculated = m.m_adjacencyCalculated;
}

//=======================================================================
//=======================================================================

Mesh::~Mesh()
{
}

//=======================================================================
//=======================================================================

void Mesh::load( const char *fileName )
{
	svPrint("loading %s\n",fileName);

	m_vertices.clear();
	m_triangles.clear();
	m_edges.clear();
	m_surfaces.clear();

	m_ambiguousNeighbours = 0;
	m_adjacencyCalculated = false;

	m_filename = fileName;

	FILE	*pF;
	if( (pF = fopen( fileName, "rb" )) != NULL )
	{
		fseek( pF, 0L, SEEK_END );
		UINT32 l = ftell( pF );
		fseek( pF, 0L, SEEK_SET );

		m_fileptr = new UINT8[l];

		if( (fread( m_fileptr, sizeof(char), l, pF )) != l)
		{
			svError( STR("error while reading file '%s'",fileName) );
		}
		fclose( pF );
	}
	else
	{
		svError( STR("couldn't open file '%s'",fileName) );
	}
	m_readptr = m_fileptr;

	if(readID4() != 'FORM')
	{
		svError("not a valid IFF file");
	}

	m_fileend = m_fileptr + readU4();

	UINT32 i = readID4();
	if(i == 'LWOB') loadLW5();
	else if(i == 'LWO2') loadLW6();
	else svError("not a LWO file");

	delete[] m_fileptr;

	svPrint("polys %d vertices %d surfaces %d\n",m_triangles.size(),m_vertices.size(),m_surfaces.size());
}

//=======================================================================
//=======================================================================

void Mesh::calculateAdjacency()
{
	if( m_adjacencyCalculated ) return;

	svPrint("calculating polygon adjacency\n");

	std::set<Edge> edgeset;

	m_edges.clear();

	//create a set of edges + calculate neighbourhood information
	//O(T log E)
	m_ambiguousNeighbours = 0;
	for(int i=0;i<m_triangles.size();++i)
	{
		for(int j=0;j<3;j++)
		{
			m_triangles[i].m_n[j] = -1;

			int k = (j+1)%3;
			Edge edge;
			edge.m_v[0] = min2(m_triangles[i].m_v[j],m_triangles[i].m_v[k]);
			edge.m_v[1] = max2(m_triangles[i].m_v[j],m_triangles[i].m_v[k]);
			edge.m_t[0] = i;
			edge.m_t[1] = -1;
			std::pair<std::set<Edge>::iterator,bool> ret = edgeset.insert( edge );	//O(log E)
			if( !ret.second )
			{	//already exists, check if it already has 2 triangles
				if( (*ret.first).m_t[1] != -1 ) ++m_ambiguousNeighbours;
				else
				{	//neighbour found, update information
					(*ret.first).m_t[1] = i;
					int n = (*ret.first).m_t[0];
					int ne;
					if( edge.m_v[0] == m_triangles[n].m_v[0] )
					{
						if( edge.m_v[1] == m_triangles[n].m_v[1] ) ne = 0;
						else if( edge.m_v[1] == m_triangles[n].m_v[2] ) ne = 2;
					}
					else if( edge.m_v[0] == m_triangles[n].m_v[1] )
					{
						if( edge.m_v[1] == m_triangles[n].m_v[0] ) ne = 0;
						else if( edge.m_v[1] == m_triangles[n].m_v[2] ) ne = 1;
					}
					else if( edge.m_v[0] == m_triangles[n].m_v[2] )
					{
						if( edge.m_v[1] == m_triangles[n].m_v[0] ) ne = 2;
						else if( edge.m_v[1] == m_triangles[n].m_v[1] ) ne = 1;
					}
					m_triangles[i].m_n[j] = n;
					m_triangles[i].m_nedge[j] = ne;
					m_triangles[n].m_n[ne] = i;
					m_triangles[n].m_nedge[ne] = j;
				}
			}
		}
	}

	//copy edges to vector
	//O( E )
	int i=0;
	int interior = 0,boundary = 0;
	m_edges.resize( edgeset.size() );
	std::set<Edge>::iterator it = edgeset.begin();
	while( it != edgeset.end() )
	{
		m_edges[i] = (*it);
		if( m_edges[i].m_t[1] == -1 ) ++boundary;
		else ++interior;
		++it;
		++i;
	}

	svPrint("edges %d interior edges %d boundary edges %d\n",m_edges.size(),interior,boundary);
	svPrint("edges shared by more than 2 triangles: %d\n",m_ambiguousNeighbours);
	m_adjacencyCalculated = true;
}

//=======================================================================
//=======================================================================

void Mesh::setParametrization()
{
	for(int t=0;t<m_triangles.size();t++)
	{
		m_triangles[t].m_uv[0].set(0,0);
		m_triangles[t].m_uv[1].set(1,0);
		m_triangles[t].m_uv[2].set(0,1);
		m_triangles[t].m_baseTriangle = t;
	}
}

//=======================================================================
//=======================================================================

void Mesh::calculateNormals()
{
	//recalculate triangle and vertex normals
	for(int i=0;i<m_vertices.size();i++)
		m_vertices[i].m_normal.set(0,0,0);
	float area;
	for(int t=0;t<m_triangles.size();t++)
	{
		Vector3 n,e1,e2;
		e1 = m_vertices[m_triangles[t].m_v[1]].m_position;
		e1 -= m_vertices[m_triangles[t].m_v[0]].m_position;
		e2 = m_vertices[m_triangles[t].m_v[2]].m_position;
		e2 -= m_vertices[m_triangles[t].m_v[0]].m_position;
		n = cross(e1, e2);
		area = n.length();
		m_triangles[t].m_area = 0.5f * area;
		n /= area;
		m_triangles[t].m_normal = n;

		m_vertices[m_triangles[t].m_v[0]].m_normal += n;
		m_vertices[m_triangles[t].m_v[1]].m_normal += n;
		m_vertices[m_triangles[t].m_v[2]].m_normal += n;
	}
	for(int i=0;i<m_vertices.size();i++)
		m_vertices[i].m_normal.normalize();
}

//=======================================================================
//=======================================================================

void Mesh::setTriangleUsing()
{
	for(int i=0;i<m_vertices.size();i++)
	{
		m_vertices[i].m_triangleUsing = -1;
	}
	for(int i=0;i<m_triangles.size();++i)
	{
		m_vertices[m_triangles[i].m_v[0]].m_triangleUsing = (i<<2) | 0;
		m_vertices[m_triangles[i].m_v[1]].m_triangleUsing = (i<<2) | 1;
		m_vertices[m_triangles[i].m_v[2]].m_triangleUsing = (i<<2) | 2;
	}
}

//=======================================================================
//=======================================================================

void Mesh::subdivide( int n )
{
	if( n == 0 ) return;
	SS_ASSERT(n > 0);

	setTriangleUsing();
	calculateAdjacency();
	if( m_ambiguousNeighbours ) svError("can't subdivide: '%s' contains %d edges which are shared by more than 2 polygons",m_filename.c_str(),m_ambiguousNeighbours);

	std::vector<Vertex>		vertices2;
	std::vector<Triangle>	triangles2;
	std::vector<Edge>		edges2;
	std::vector<Vertex>		*vIn,*vOut;
	std::vector<Triangle>	*tIn,*tOut;
	std::vector<Edge>		*eIn,*eOut;

	vIn = &m_vertices;
	tIn = &m_triangles;
	eIn = &m_edges;
	vOut = &vertices2;
	tOut = &triangles2;
	eOut = &edges2;

//t(n+1) = t(n) * 4
//e(n+1) = e(n) * 2 + t(n) * 3
//v(n+1) = v(n) + e(n)

//t(n) = t(0) * 2^2n   , n>=0
//e(n) = e(0) * 2^n + t(0) * 3 * 2^(n-1) * (2^n - 1)   , n>=0
//v(n) = v(0) + e(0) * (2^n - 1) + t(0) * 3 * (2^2n / 6 - 2^(n-1) + 1)  , n > 0

	int ct=m_triangles.size(),cv=m_vertices.size(),ce=m_edges.size();
	cv += ce * ((1<<n)-1);
	if(n > 1) cv += ct * 3 * ((1<<(2*n)) / 6 + 1 - (1<<(n-1)));
	ce *= (1<<n);
	ce += ct * 3 * (1<<(n-1)) * ((1<<n) - 1);
	ct *= (1<<(2*n));
	//TODO reserve (other array (n-1), other n)
//	svPrint("subdivision reserve      v:%d t:%d e:%d\n",cv,ct,ce);

	svPrint("subdivision before       v:%d t:%d e:%d\n",m_vertices.size(),m_triangles.size(),m_edges.size());

	for(int step=0;step<n;step++)
	{
		Edge en;

		*vOut = *vIn;
		tOut->clear();
		eOut->clear();

		for(int t=0;t<tIn->size();t++)
		{
			//calculate edge midpoints
			int nv[3];
			for(int e=0;e<3;e++)
			{
				int i1,i2,i3,a3,nbour,nedge;
				Vertex vn;

				i1 = (*tIn)[t].m_v[ (e+0)%3 ];
				i2 = (*tIn)[t].m_v[ (e+1)%3 ];
				i3 = (*tIn)[t].m_v[ (e+2)%3 ];

				nv[e] = vOut->size();
				vn.m_triangleUsing = ((t * 4 + e) << 2) | ((e+1)%3);
				vn.reset();

				nbour = (*tIn)[t].m_n[e];
				if( nbour != -1 )
				{
					nedge = (*tIn)[t].m_nedge[e];
					if( nbour < t )
					{
						//edge is split already
						nv[e] = (*tOut)[ nbour * 4 + nedge ].m_v[ (nedge+1)%3 ];
					}
					else
					{
						//new vertex
						a3 = (*tIn)[ nbour ].m_v[ (nedge+2)%3 ];

						vn.mad( (*vIn)[i1], 0.375f );
						vn.mad( (*vIn)[i2], 0.375f );
						vn.mad( (*vIn)[i3], 0.125f );
						vn.mad( (*vIn)[a3], 0.125f );
						vn.normalize();
						vOut->push_back( vn );

						en.m_t[0] = t * 4 + e;
						en.m_t[1] = nbour * 4 + ((nedge + 1) % 3);
						en.m_v[0] = i1;
						en.m_v[1] = nv[e];
						eOut->push_back( en );

						en.m_t[0] = t * 4 + ((e + 1) % 3);
						en.m_t[1] = nbour * 4 + nedge;
						en.m_v[0] = nv[e];
						en.m_v[1] = i2;
						eOut->push_back( en );
					}
				}
				else //boundary
				{
					vn.mad( (*vIn)[i1], 0.5f );
					vn.mad( (*vIn)[i2], 0.5f );
					vn.normalize();
					vOut->push_back( vn );

					en.m_t[0] = t * 4 + e;
					en.m_t[1] = -1;
					en.m_v[0] = i1;
					en.m_v[1] = nv[e];
					eOut->push_back( en );

					en.m_t[0] = t * 4 + ((e + 1) % 3);
					en.m_t[1] = -1;
					en.m_v[0] = nv[e];
					en.m_v[1] = i2;
					eOut->push_back( en );
				}
			}

			//make new triangles
			Triangle tn;
			int ct = tOut->size();
			tn.m_surface = (*tIn)[t].m_surface;
			tn.m_baseTriangle = (*tIn)[t].m_baseTriangle;

			tn.m_v[0] = (*tIn)[t].m_v[0];
			tn.m_v[1] = nv[0];
			tn.m_v[2] = nv[2];
			tn.m_n[0] = (*tIn)[t].m_n[0] != -1 ? (*tIn)[t].m_n[0] * 4 + ((*tIn)[t].m_nedge[0] + 1) % 3 : -1;
			tn.m_n[1] = ct + 3;
			tn.m_n[2] = (*tIn)[t].m_n[2] != -1 ? (*tIn)[t].m_n[2] * 4 + (*tIn)[t].m_nedge[2]: -1;
			tn.m_uv[0] = (*tIn)[t].m_uv[0];
			tn.m_uv[1] = 0.5f * ((*tIn)[t].m_uv[1] + (*tIn)[t].m_uv[0]);
			tn.m_uv[2] = 0.5f * ((*tIn)[t].m_uv[2] + (*tIn)[t].m_uv[0]);
			tn.m_nedge[0] = (*tIn)[t].m_nedge[0];
			tn.m_nedge[1] = 0;
			tn.m_nedge[2] = (*tIn)[t].m_nedge[2];
			tOut->push_back( tn );

			tn.m_v[0] = nv[0];
			tn.m_v[1] = (*tIn)[t].m_v[1];
			tn.m_v[2] = nv[1];
			tn.m_n[0] = (*tIn)[t].m_n[0] != -1 ? (*tIn)[t].m_n[0] * 4 + (*tIn)[t].m_nedge[0]: -1;
			tn.m_n[1] = (*tIn)[t].m_n[1] != -1 ? (*tIn)[t].m_n[1] * 4 + ((*tIn)[t].m_nedge[1] + 1) % 3 : -1;
			tn.m_n[2] = ct + 3;
			tn.m_uv[0] = 0.5f * ((*tIn)[t].m_uv[0] + (*tIn)[t].m_uv[1]);
			tn.m_uv[1] = (*tIn)[t].m_uv[1];
			tn.m_uv[2] = 0.5f * ((*tIn)[t].m_uv[2] + (*tIn)[t].m_uv[1]);
			tn.m_nedge[0] = (*tIn)[t].m_nedge[0];
			tn.m_nedge[1] = (*tIn)[t].m_nedge[1];
			tn.m_nedge[2] = 1;
			tOut->push_back( tn );

			tn.m_v[0] = nv[2];
			tn.m_v[1] = nv[1];
			tn.m_v[2] = (*tIn)[t].m_v[2];
			tn.m_n[0] = ct + 3;
			tn.m_n[1] = (*tIn)[t].m_n[1] != -1 ? (*tIn)[t].m_n[1] * 4 + (*tIn)[t].m_nedge[1]: -1;
			tn.m_n[2] = (*tIn)[t].m_n[2] != -1 ? (*tIn)[t].m_n[2] * 4 + ((*tIn)[t].m_nedge[2] + 1) % 3 : -1;
			tn.m_uv[0] = 0.5f * ((*tIn)[t].m_uv[0] + (*tIn)[t].m_uv[2]);
			tn.m_uv[1] = 0.5f * ((*tIn)[t].m_uv[1] + (*tIn)[t].m_uv[2]);
			tn.m_uv[2] = (*tIn)[t].m_uv[2];
			tn.m_nedge[0] = 2;
			tn.m_nedge[1] = (*tIn)[t].m_nedge[1];
			tn.m_nedge[2] = (*tIn)[t].m_nedge[2];
			tOut->push_back( tn );

			tn.m_v[0] = nv[2];
			tn.m_v[1] = nv[0];
			tn.m_v[2] = nv[1];
			tn.m_n[0] = ct;
			tn.m_n[1] = ct + 1;
			tn.m_n[2] = ct + 2;
			tn.m_uv[0] = 0.5f * ((*tIn)[t].m_uv[2] + (*tIn)[t].m_uv[0]);
			tn.m_uv[1] = 0.5f * ((*tIn)[t].m_uv[0] + (*tIn)[t].m_uv[1]);
			tn.m_uv[2] = 0.5f * ((*tIn)[t].m_uv[2] + (*tIn)[t].m_uv[1]);
			tn.m_nedge[0] = 1;
			tn.m_nedge[1] = 2;
			tn.m_nedge[2] = 0;
			tOut->push_back( tn );

			en.m_t[0] = ct;
			en.m_t[1] = ct + 3;
			en.m_v[0] = nv[0];
			en.m_v[1] = nv[2];
			eOut->push_back( en );

			en.m_t[0] = ct + 1;
			en.m_t[1] = ct + 3;
			en.m_v[0] = nv[1];
			en.m_v[1] = nv[0];
			eOut->push_back( en );

			en.m_t[0] = ct + 2;
			en.m_t[1] = ct + 3;
			en.m_v[0] = nv[2];
			en.m_v[1] = nv[1];
			eOut->push_back( en );
		}

		//move even vertices
		for(int v=0;v<vIn->size();v++)
		{
			if( (*vIn)[v].m_triangleUsing == -1 ) continue;	//not used

			Vertex sum;
			int spoly = (*vIn)[v].m_triangleUsing >> 2,cpoly,npoly,lB,rB;
			int sedge = (*vIn)[v].m_triangleUsing & 3,cedge,nedge,cv;
			SS_ASSERT( sedge >= 0 && sedge < 3 );
			SS_ASSERT( spoly >= 0 && spoly < tIn->size() );
			int valence = 0;
			bool boundary = false;

			sum.reset();

			//loop right
			npoly = spoly;
			nedge = sedge;
			do
			{
				cpoly = npoly;
				cedge = nedge;

				//loop code
				valence++;
				cv = (*tIn)[cpoly].m_v[(cedge+1) % 3];
				sum.add( (*vIn)[cv] );

				npoly = (*tIn)[cpoly].m_n[cedge];
				if( npoly == -1 )
				{
					boundary = true;
					rB = cv;
					break;
				}
				nedge = ((*tIn)[cpoly].m_nedge[cedge] + 1) % 3;
			} while( npoly != spoly );
			if( boundary )
			{
				npoly = spoly;
				nedge = (sedge + 2)%3;
				//loop left
				do
				{
					cpoly = npoly;
					cedge = nedge;

					//loop code
					valence++;
					cv = (*tIn)[cpoly].m_v[cedge];
					sum.add( (*vIn)[cv] );

					npoly = (*tIn)[cpoly].m_n[cedge];
					if( npoly == -1 )
					{
						lB = cv;
						break;
					}
					nedge = ((*tIn)[cpoly].m_nedge[cedge] + 2) % 3;
				} while( npoly != spoly );


				(*vOut)[v].reset();
				(*vOut)[v].mad( (*vIn)[v], 0.75f );
				(*vOut)[v].mad( (*vIn)[lB], 0.125f );
				(*vOut)[v].mad( (*vIn)[rB], 0.125f );
				(*vOut)[v].normalize();
				(*vOut)[v].m_triangleUsing = ((spoly * 4 + sedge) << 2) | sedge;
			}
			else	//no boundary
			{
				float beta;
				if( valence == 3 )
					beta = 3.0f / 16.0f;
				else
					beta = 3.0f / (float(valence) * 8.0f);

				(*vOut)[v].reset();
				(*vOut)[v].mad( (*vIn)[v], (1.0f - float(valence) * beta) );
				(*vOut)[v].mad( sum, beta );
				(*vOut)[v].normalize();
				(*vOut)[v].m_triangleUsing = ((spoly * 4 + sedge) << 2) | sedge;
			}
		}

    std::swap(vIn,vOut);
		std::swap(tIn,tOut);
		std::swap(eIn,eOut);

		svPrint("subdivision after step %d v:%d t:%d e:%d\n",(step+1),vIn->size(),tIn->size(),eIn->size());
	}

	if(vIn != &m_vertices)
	{
		m_vertices = *vIn;
		m_triangles = *tIn;
		m_edges = *eIn;
	}

	calculateNormals();
}

//=======================================================================
//=======================================================================

const std::vector<Mesh::ventry> &Mesh::get1Neighbourhood( int ver, int &valence, bool &boundary, int &left, int &right )
{
	static std::vector<ventry> vertexlist;
	valence = 0;
	boundary = false;
	vertexlist.clear();
	if( m_vertices[ver].m_triangleUsing == -1 ) return vertexlist;	//not used
	if( m_ambiguousNeighbours ) svError("can't traverse neighbourhood: '%s' contains %d edges which are shared by more than 2 polygons",m_filename.c_str(),m_ambiguousNeighbours);

	int spoly = m_vertices[ver].m_triangleUsing >> 2,cpoly,npoly;
	int sedge = m_vertices[ver].m_triangleUsing & 3,cedge,nedge;
	SS_ASSERT( sedge >= 0 && sedge < 3 );
	SS_ASSERT( spoly >= 0 && spoly < m_triangles.size() );
	int i;
	ventry ve;

	//these loops count on that there are no neighbour ambiguities (edges with >2 polygons)
	// if this condition doesn't hold, these may loop infinitely
	//loop right
	npoly = spoly;
	nedge = sedge;
	i = 0;
	do
	{
		cpoly = npoly;
		cedge = nedge;

		valence++;
		ve.m_vertex = m_triangles[cpoly].m_v[(cedge+1) % 3];
		ve.m_i = i;
		vertexlist.push_back( ve );
		i++;

		npoly = m_triangles[cpoly].m_n[cedge];
		if( npoly == -1 )
		{
			boundary = true;
			right = ve.m_vertex;
			break;
		}
		nedge = (m_triangles[cpoly].m_nedge[cedge] + 1) % 3;
	} while( npoly != spoly );
	if( boundary )
	{
		i = -1;
		npoly = spoly;
		nedge = (sedge + 2)%3;
		//loop left
		do
		{
			cpoly = npoly;
			cedge = nedge;

			valence++;
			ve.m_vertex = m_triangles[cpoly].m_v[cedge];
			ve.m_i = i;
			vertexlist.push_back( ve );
			i--;

			npoly = m_triangles[cpoly].m_n[cedge];
			if( npoly == -1 )
			{
				left = ve.m_vertex;
				break;
			}
			nedge = (m_triangles[cpoly].m_nedge[cedge] + 2) % 3;
		} while( npoly != spoly );
	}
	return vertexlist;
}

//=======================================================================
//=======================================================================

void Mesh::calculateBoundingBox()
{
	m_boundingBoxMin.set(FLT_MAX,FLT_MAX,FLT_MAX);
	m_boundingBoxMax.set(-FLT_MAX,-FLT_MAX,-FLT_MAX);
	m_boundingRadius = 0.0f;
	for(int i=0;i<m_vertices.size();i++)
	{
		Vector3 v = m_vertices[i].m_position;
		if( v.x > m_boundingBoxMax.x ) m_boundingBoxMax.x = v.x;
		if( v.y > m_boundingBoxMax.y ) m_boundingBoxMax.y = v.y;
		if( v.z > m_boundingBoxMax.z ) m_boundingBoxMax.z = v.z;
		if( v.x < m_boundingBoxMin.x ) m_boundingBoxMin.x = v.x;
		if( v.y < m_boundingBoxMin.y ) m_boundingBoxMin.y = v.y;
		if( v.z < m_boundingBoxMin.z ) m_boundingBoxMin.z = v.z;

		if( dot(v,v) > m_boundingRadius ) m_boundingRadius = dot(v,v);
	}
	m_boundingRadius = sqrt(m_boundingRadius);
}

//=======================================================================
//=======================================================================
