#include "dxstdafx.h"
#include "trianglemesh.h"
#include "Radion.hpp"

TriangleMesh::TriangleMesh(Material* material, unsigned short* indexBuffer, unsigned int nFaces,
						   D3DVERTEX* vertexBuffer, unsigned int nVertices)
{
	this->material = material;
	nMeshVertices = nVertices;
	nMeshPatches = nFaces;

	meshVertices = new Vector[nMeshVertices];
	normals = new Vector[nMeshVertices];
	texCoords = new Vector[nMeshVertices];
	for(int u=0; u < nMeshVertices; u++)
	{
		meshVertices[u] = vertexBuffer[u].pos;
		normals[u] = vertexBuffer[u].normal;
		texCoords[u] = vertexBuffer[u].tex0;
	}
	
	meshPatches = new Patch[nMeshPatches];

	int pup = 0;
	for(int p = 0; p < nMeshPatches; p++)
	{
		meshPatches[p].vertexIndices[0] = indexBuffer[pup++];
		meshPatches[p].vertexIndices[1] = indexBuffer[pup++];
		meshPatches[p].vertexIndices[2] = indexBuffer[pup++];
	}
	//calculate normals
	for(p = 0; p < nMeshPatches; p++)
	{
		meshPatches[p].flatNormal = (
			meshVertices[meshPatches[p].vertexIndices[0]] -
			meshVertices[meshPatches[p].vertexIndices[1]]
			)
			&&
			(
			meshVertices[meshPatches[p].vertexIndices[0]] -
			meshVertices[meshPatches[p].vertexIndices[2]]
			);
		meshPatches[p].flatNormal.normalize();
		normals[meshPatches[p].vertexIndices[0]] += meshPatches[p].flatNormal;
		normals[meshPatches[p].vertexIndices[1]] += meshPatches[p].flatNormal;
		normals[meshPatches[p].vertexIndices[2]] += meshPatches[p].flatNormal;
		meshPatches[p].hyperPlaneShiftOffset = 
			meshVertices[meshPatches[p].vertexIndices[0]] * 
						meshPatches[p].flatNormal;

		Vector A[3];
		A[0] = meshVertices[meshPatches[p].vertexIndices[0]];
		A[1] = meshVertices[meshPatches[p].vertexIndices[1]];
		A[2] = meshVertices[meshPatches[p].vertexIndices[2]];
		float		t4 = A[0][0]*A[1][1];
		float 	    t6 = A[0][0]*A[1][2];
		float 	    t8 = A[0][1]*A[1][0];
		float 	    t10 = A[0][2]*A[1][0];
		float 	    t12 = A[0][1]*A[2][0];
		float 	    t14 = A[0][2]*A[2][0];
		float 	    t17 = 1/(t4*A[2][2]-t6*A[2][1]-t8*A[2][2]+t10*A[2][1]+t12*A[1][2]-t14*A[1][1]);
	
//		if(_isnan (t17) )
//          AfxMessageBox("mtx inversion gbz");

		Vector* m = meshPatches[p].inverseVertexMatrix;
		
		m[0][0] = (A[1][1]*A[2][2]-A[1][2]*A[2][1])*t17;
		m[1][0] = -(A[0][1]*A[2][2]-A[0][2]*A[2][1])*t17;
		m[2][0] = -(-A[0][1]*A[1][2]+A[0][2]*A[1][1])*t17;
		m[0][1] = -(A[1][0]*A[2][2]-A[1][2]*A[2][0])*t17;
		m[1][1] = (A[0][0]*A[2][2]-t14)*t17;
		m[2][1] = -(t6-t10)*t17;
		m[0][2] = -(-A[1][0]*A[2][1]+A[1][1]*A[2][0])*t17;
		m[1][2] = -(A[0][0]*A[2][1]-t12)*t17;
		m[2][2] = (t4-t8)*t17;

		meshPatches[p].bbox.minPoint = meshVertices[meshPatches[p].vertexIndices[0]];
		meshPatches[p].bbox.minPoint <= meshVertices[meshPatches[p].vertexIndices[1]];
		meshPatches[p].bbox.minPoint <= meshVertices[meshPatches[p].vertexIndices[2]];

		meshPatches[p].bbox.maxPoint = meshVertices[meshPatches[p].vertexIndices[0]];
		meshPatches[p].bbox.maxPoint >= meshVertices[meshPatches[p].vertexIndices[1]];
		meshPatches[p].bbox.maxPoint >= meshVertices[meshPatches[p].vertexIndices[2]];
	}
	for(int n=0; n<nMeshVertices; n++)
	{
		normals[n].normalize();
	}
	Intersectable** objs = new Intersectable*[nMeshPatches];
	for(int t=0; t<nMeshPatches; t++)
		objs[t] = meshPatches + t;
	meshTree = new KDTree(objs, nMeshPatches);
	bbox = meshTree->getBoundingBox();
	delete objs;

	buildAreaTree();
}

TriangleMesh::TriangleMesh(std::istream& isc, Material** materialTable, int nMaterials)
{
	char keyword[100];
	isc >> keyword; //material
	isc >> keyword;
	for(int i=0; i < nMaterials; i++)
	{
		material = materialTable[i];
		if(strcmp(keyword, material->getName()) == 0) break;
	}
	isc >> keyword;
	isc >> nMeshVertices;
	meshVertices = new Vector[nMeshVertices];
	normals = new Vector[nMeshVertices];
	for(int v=0; v<nMeshVertices; v++)
	{
		float a, b, c;
		isc >> a >> b >> c;
//		meshVertices[v] = Vector(a + (0.1f * (float)rand() / RAND_MAX),
//			b + (0.1f * (float)rand() / RAND_MAX),
//			c + (0.1f * (float)rand() / RAND_MAX));
		meshVertices[v] = Vector(a , b, c);
		normals[v].clear();
	}
	bbox.minPoint = meshVertices[0];
	bbox.maxPoint = meshVertices[0];
	for(int w=1; w<nMeshVertices; w++)
	{
		bbox.minPoint <= meshVertices[w];
		bbox.maxPoint >= meshVertices[w];
	}
	isc >> keyword;
	isc >> nMeshPatches;
	meshPatches = new Patch[nMeshPatches];
	int p;
	for(p = 0; p < nMeshPatches; p++)
	{
		signed int ai, bi, ci;
		isc >> ai >> bi >> ci;
//		if(di < 0)
		{
			meshPatches[p].vertexIndices[0] = ai;
			meshPatches[p].vertexIndices[1] = bi;
			meshPatches[p].vertexIndices[2] = ci;
		}
/*		else	//quad, tessellate
		{
			meshPatches[p].vertexIndices[0] = ai;
			meshPatches[p].vertexIndices[1] = bi;
			meshPatches[p].vertexIndices[2] = ci;
			p++;
			meshPatches[p].vertexIndices[0] = ai;
			meshPatches[p].vertexIndices[1] = ci;
			meshPatches[p].vertexIndices[2] = di;
			isc >> ai; //final -1
		}*/
	}
	//calculate normals
	for(p = 0; p < nMeshPatches; p++)
	{
		meshPatches[p].flatNormal = (
			meshVertices[meshPatches[p].vertexIndices[0]] -
			meshVertices[meshPatches[p].vertexIndices[1]]
			)
			&&
			(
			meshVertices[meshPatches[p].vertexIndices[0]] -
			meshVertices[meshPatches[p].vertexIndices[2]]
			);
		meshPatches[p].flatNormal.normalize();
		normals[meshPatches[p].vertexIndices[0]] += meshPatches[p].flatNormal;
		normals[meshPatches[p].vertexIndices[1]] += meshPatches[p].flatNormal;
		normals[meshPatches[p].vertexIndices[2]] += meshPatches[p].flatNormal;
		meshPatches[p].hyperPlaneShiftOffset = 
			meshVertices[meshPatches[p].vertexIndices[0]] * 
						meshPatches[p].flatNormal;

		Vector A[3];
		A[0] = meshVertices[meshPatches[p].vertexIndices[0]];
		A[1] = meshVertices[meshPatches[p].vertexIndices[1]];
		A[2] = meshVertices[meshPatches[p].vertexIndices[2]];
		float		t4 = A[0][0]*A[1][1];
		float 	    t6 = A[0][0]*A[1][2];
		float 	    t8 = A[0][1]*A[1][0];
		float 	    t10 = A[0][2]*A[1][0];
		float 	    t12 = A[0][1]*A[2][0];
		float 	    t14 = A[0][2]*A[2][0];
		float 	    t17 = 1/(t4*A[2][2]-t6*A[2][1]-t8*A[2][2]+t10*A[2][1]+t12*A[1][2]-t14*A[1][1]);
	
//		if(_isnan (t17) )
//          AfxMessageBox("mtx inversion gbz");

		Vector* m = meshPatches[p].inverseVertexMatrix;
		
		m[0][0] = (A[1][1]*A[2][2]-A[1][2]*A[2][1])*t17;
		m[1][0] = -(A[0][1]*A[2][2]-A[0][2]*A[2][1])*t17;
		m[2][0] = -(-A[0][1]*A[1][2]+A[0][2]*A[1][1])*t17;
		m[0][1] = -(A[1][0]*A[2][2]-A[1][2]*A[2][0])*t17;
		m[1][1] = (A[0][0]*A[2][2]-t14)*t17;
		m[2][1] = -(t6-t10)*t17;
		m[0][2] = -(-A[1][0]*A[2][1]+A[1][1]*A[2][0])*t17;
		m[1][2] = -(A[0][0]*A[2][1]-t12)*t17;
		m[2][2] = (t4-t8)*t17;

		meshPatches[p].bbox.minPoint = meshVertices[meshPatches[p].vertexIndices[0]];
		meshPatches[p].bbox.minPoint <= meshVertices[meshPatches[p].vertexIndices[1]];
		meshPatches[p].bbox.minPoint <= meshVertices[meshPatches[p].vertexIndices[2]];

		meshPatches[p].bbox.maxPoint = meshVertices[meshPatches[p].vertexIndices[0]];
		meshPatches[p].bbox.maxPoint >= meshVertices[meshPatches[p].vertexIndices[1]];
		meshPatches[p].bbox.maxPoint >= meshVertices[meshPatches[p].vertexIndices[2]];
	}
	for(int n=0; n<nMeshVertices; n++)
	{
		normals[n].normalize();
	}
	Intersectable** objs = new Intersectable*[nMeshPatches];
	for(int t=0; t<nMeshPatches; t++)
		objs[t] = meshPatches + t;
	meshTree = new KDTree(objs, nMeshPatches);
	bbox = meshTree->getBoundingBox();
	delete objs;
}

void TriangleMesh::getTransformedBoundingBox(const Transformation& tf, BoundingBox& bb)
{
	Vector tfdvec;
	tf.transformPoint(meshVertices[0], tfdvec);
	bb.minPoint = tfdvec;
	bb.maxPoint = tfdvec;
	for(int w=1; w<nMeshVertices; w++)
	{
		tf.transformPoint(meshVertices[w], tfdvec);
		bb.minPoint <= tfdvec;
		bb.maxPoint >= tfdvec;
	}
	bb.minPoint.x -= 0.01;
	bb.minPoint.y -= 0.01;
	bb.minPoint.z -= 0.01;
	bb.maxPoint.x += 0.01;
	bb.maxPoint.y += 0.01;
	bb.maxPoint.z += 0.01;
}

bool TriangleMesh::Patch::intersectBackSide (const Ray& ray, float& depth, float rayMin, float rayMax)
{
	lastTestedRayId = ray.id;
	lastTestedRayResult.isIntersect = false;

	float cosa = flatNormal * ray.dir;
	if (cosa < 0.00001f)	// front facing triangle
		return false;

	float originDistOnNormal = -(flatNormal * ray.origin);
	depth = (hyperPlaneShiftOffset + originDistOnNormal) / cosa;
	if (depth < 0.01f)
		return false;

	Vector hitPoint = ray.origin;
	hitPoint.addScaled(depth, ray.dir);

	float baryA = hitPoint * inverseVertexMatrix[0];
	if(baryA < -0.0001f) return false;
	float baryB = hitPoint * inverseVertexMatrix[1];
	if(baryB < -0.0001f) return false;
	float baryC = hitPoint * inverseVertexMatrix[2];
	if(baryC < -0.0001f) return false;

	if(ray.isShadowRay)
	{
		//faster way to tell
		lastTestedRayResult.isIntersect	= true;
		lastTestedRayResult.depth		= depth;
		return true;
	}

	lastTestedRayResult.normal.setScaled(baryA, meshNormals[vertexIndices[0]]);
	lastTestedRayResult.normal.addScaled(baryB, meshNormals[vertexIndices[1]]);
	lastTestedRayResult.normal.addScaled(baryC, meshNormals[vertexIndices[2]]);

	lastTestedRayResult.texUV.setScaled(baryA, meshTexCoords[vertexIndices[0]]);
	lastTestedRayResult.texUV.addScaled(baryB, meshTexCoords[vertexIndices[1]]);
	lastTestedRayResult.texUV.addScaled(baryC, meshTexCoords[vertexIndices[2]]);

	lastTestedRayResult.isIntersect	= true;
	lastTestedRayResult.depth		= depth;
	lastTestedRayResult.point		= hitPoint;
	lastTestedRayResult.object		= this;
	return true;
}

bool TriangleMesh::Patch::intersect (const Ray& ray, float& depth, float rayMin, float rayMax)
{
	lastTestedRayId = ray.id;
	lastTestedRayResult.isIntersect = false;

	float cosa = flatNormal * ray.dir;
	if (cosa > -0.00001f)	// back facing triangle
		return false;

	float originDistOnNormal = -(flatNormal * ray.origin);
	depth = (hyperPlaneShiftOffset + originDistOnNormal) / cosa;
	if (depth < 0.0f)
		return false;

	Vector hitPoint = ray.origin;
	hitPoint.addScaled(depth, ray.dir);

	float baryA = hitPoint * inverseVertexMatrix[0];
	if(baryA < -0.1f) return false;
	float baryB = hitPoint * inverseVertexMatrix[1];
	if(baryB <  -0.1f) return false;
	float baryC = hitPoint * inverseVertexMatrix[2];
	if(baryC <  -0.1f) return false;

	if(ray.isShadowRay)
	{
		//faster way to tell
		lastTestedRayResult.isIntersect	= true;
		lastTestedRayResult.depth		= depth;
		return true;
	}

	lastTestedRayResult.normal.setScaled(baryA, meshNormals[vertexIndices[0]]);
	lastTestedRayResult.normal.addScaled(baryB, meshNormals[vertexIndices[1]]);
	lastTestedRayResult.normal.addScaled(baryC, meshNormals[vertexIndices[2]]);
//	lastTestedRayResult.normal = flatNormal;

	lastTestedRayResult.texUV.setScaled(baryA, meshTexCoords[vertexIndices[0]]);
	lastTestedRayResult.texUV.addScaled(baryB, meshTexCoords[vertexIndices[1]]);
	lastTestedRayResult.texUV.addScaled(baryC, meshTexCoords[vertexIndices[2]]);

	lastTestedRayResult.isIntersect	= true;
	lastTestedRayResult.depth		= depth;
	lastTestedRayResult.point		= hitPoint;
	return true;
}

void TriangleMesh::Patch::sampleSurface(Radion& radion)
{
	Vector u = meshVertices[vertexIndices[1]] - meshVertices[vertexIndices[0]];
	Vector v = meshVertices[vertexIndices[2]] - meshVertices[vertexIndices[0]];
	double r1 = (double)rand() / RAND_MAX;
	double r2 = (double)rand() / RAND_MAX;
	if(r1 + r2 > 1.0)
	{
		r1 = 1.0 - r1;
		r2 = 1.0 - r2;
	}
	radion.position = meshVertices[vertexIndices[0]] + u * r1 + v * r2;
	float baryA = (r1 + r2) * 0.5f;
	float baryB = (1.0f - r1) * 0.5f;
	float baryC = (1.0f - r2) * 0.5f;

	radion.normal.setScaled(baryA, meshNormals[vertexIndices[0]]);
	radion.normal.addScaled(baryB, meshNormals[vertexIndices[1]]);
	radion.normal.addScaled(baryC, meshNormals[vertexIndices[2]]);

	radion.radiance.setScaled(baryA, meshTexCoords[vertexIndices[0]]);
	radion.radiance.addScaled(baryB, meshTexCoords[vertexIndices[1]]);
	radion.radiance.addScaled(baryC, meshTexCoords[vertexIndices[2]]);

//	radion.radiance.z = getSurfaceArea();
}

Vector* TriangleMesh::Patch::meshVertices = 0x0;
Vector* TriangleMesh::Patch::meshNormals = 0x0;
Vector* TriangleMesh::Patch::meshTexCoords = 0x0;

bool TriangleMesh::intersect (const Ray& ray, float& depth, float rayMin, float rayMax)
{
	lastTestedRayId = ray.id;
	lastTestedRayResult.isIntersect = false;
	Patch::meshNormals = normals;
	Patch::meshTexCoords = texCoords;
	HitRec hitRec;
	meshTree->traverse(ray, hitRec, rayMin, rayMax);
	lastTestedRayResult = hitRec;
	lastTestedRayResult.material = this->material;
	depth = lastTestedRayResult.depth;
	lastTestedRayResult.object = this;
	return lastTestedRayResult.isIntersect;
}

bool TriangleMesh::intersectBackSide (const Ray& ray, float& depth, float rayMin, float rayMax)
{
	lastTestedRayId = ray.id;
	lastTestedRayResult.isIntersect = false;
	Patch::meshNormals = normals;
	HitRec hitRec;
	meshTree->traverseBackSide(ray, hitRec, rayMin, rayMax);
//	meshTree->forbidden = hitRec.object;
	lastTestedRayResult = hitRec;
	lastTestedRayResult.material = this->material;
	depth = lastTestedRayResult.depth;
	lastTestedRayResult.object = this;
	return lastTestedRayResult.isIntersect;
}

TriangleMesh::~TriangleMesh(void)
{
	delete [] areaTree;
	delete meshTree;
	delete meshVertices;
	delete [] meshPatches;
	delete normals;
	delete texCoords;
}


double TriangleMesh::getPatchArea(unsigned int index)
{
	if(index >= nMeshPatches)
		return 0;
	Vector a = meshVertices[meshPatches[index].vertexIndices[1]] -
			meshVertices[meshPatches[index].vertexIndices[0]];
	Vector b = meshVertices[meshPatches[index].vertexIndices[2]] -
			meshVertices[meshPatches[index].vertexIndices[0]];
	return (a && b).norm();
}

double TriangleMesh::buildAreaTree(unsigned int u)
{
	if(u >= nAreaTreeNodes)
	{
		u -= nAreaTreeNodes;
		return getPatchArea(u);
	}
	areaTree[u] = buildAreaTree(u * 2 + 1) + buildAreaTree(u * 2 + 2);
	return areaTree[u];
}

void TriangleMesh::buildAreaTree()
{
	nAreaTreeNodes = 1;
	while(nAreaTreeNodes < nMeshPatches - 1)
		nAreaTreeNodes <<= 1;
	nAreaTreeNodes--;
	areaTree = new double[nAreaTreeNodes];
	surfaceArea = buildAreaTree(0);
}

void TriangleMesh::sampleSurface(unsigned int u, double rnd, Radion& radion)
{
	if(u >= nAreaTreeNodes)
	{
		u -= nAreaTreeNodes;
		if(u >= nMeshPatches)
			u = 0;
		meshPatches[u].sampleSurface(radion);
		return;
	}
	float leftweight = 0.0f;
	if(u * 2 + 1 >= nAreaTreeNodes)
		leftweight = getPatchArea(u * 2 + 1 - nAreaTreeNodes);
	else
		leftweight = areaTree[u * 2 + 1];
	if(rnd <= leftweight)
		sampleSurface(u * 2 + 1, rnd, radion);
	else
		sampleSurface(u * 2 + 2, rnd - leftweight, radion);
}

void TriangleMesh::sampleSurface(Radion& radion)
{
	Patch::meshVertices = meshVertices;
	Patch::meshNormals = normals;
	Patch::meshTexCoords = texCoords;
	double rnd = surfaceArea * (double)rand() / RAND_MAX;
	sampleSurface(0, rnd, radion);
	radion.radiance.z = surfaceArea;
}

float TriangleMesh::getSurfaceArea()
{
	return surfaceArea;
}