/*
 * pbrt source code Copyright(c) 1998-2005 Matt Pharr and Greg Humphreys
 *
 * All Rights Reserved.
 * For educational use only; commercial use expressly forbidden.
 * NO WARRANTY, express or implied, for this software.
 * (See file License.txt for complete license)
 */

// parametricshape_skeleton.cpp*
// This file includes the PBRT implementation of TriangleMesh and
// Triangle to aid in the implementation of a parametric surface.
// Unfortunately PBRT has no way of linking to already written code in
// another plugin, so that's why this wasteful cut and paste has
// occured.
// Skeleton by:  Richard Sharp sharpr@cse.ohio-state.edu

#include "shape.h"
#include "paramset.h"
#include "mc.h"

// TriangleMesh Declarations
class TriangleMesh : public Shape {
public:
	// TriangleMesh Public Methods
	TriangleMesh(const Transform &o2w, bool ro,
	             int ntris, int nverts, const int *vptr,
				 const Point *P, const Normal *N,
				 const Vector *S, const float *uv);
	~TriangleMesh();
	BBox ObjectBound() const;
	BBox WorldBound() const;
	bool CanIntersect() const { return false; }
	void Refine(vector<Reference<Shape> > &refined) const;
	friend class Triangle;
	template <class T> friend class VertexTexture;
protected:
	// TriangleMesh Data
	int ntris, nverts;
	int *vertexIndex;
	Point *p;
	Normal *n;
	Vector *s;
	float *uvs;
};
class Triangle : public Shape {
public:
	// Triangle Public Methods
	Triangle(const Transform &o2w, bool ro,
	         TriangleMesh *m, int n)
			: Shape(o2w, ro) {
		mesh = m;
		v = &mesh->vertexIndex[3*n];
		// Update created triangles stats
		static StatsCounter trisMade("Geometry",
		                             "Triangles created");
		++trisMade;
	}
	BBox ObjectBound() const;
	BBox WorldBound() const;
	bool Intersect(const Ray &ray, float *tHit,
	               DifferentialGeometry *dg) const;
	bool IntersectP(const Ray &ray) const;
	void GetUVs(float uv[3][2]) const;
	float Area() const;
	virtual void GetShadingGeometry(const Transform &obj2world,
			const DifferentialGeometry &dg,
			DifferentialGeometry *dgShading) const {
		if (!mesh->n && !mesh->s) {
			*dgShading = dg;
			return;
		}
		// Initialize _Triangle_ shading geometry with _n_ and _s_
		// Compute barycentric coordinates for point
		float b[3];
		// Initialize _A_ and _C_ matrices for barycentrics
		float uv[3][2];
		GetUVs(uv);
		float A[2][2] =
		    { { uv[1][0] - uv[0][0], uv[2][0] - uv[0][0] },
		      { uv[1][1] - uv[0][1], uv[2][1] - uv[0][1] } };
		float C[2] = { dg.u - uv[0][0], dg.v - uv[0][1] };
		if (!SolveLinearSystem2x2(A, C, &b[1])) {
			// Handle degenerate parametric mapping
			b[0] = b[1] = b[2] = 1.f/3.f;
		}
		else
			b[0] = 1.f - b[1] - b[2];
		// Use _n_ and _s_ to compute shading tangents for triangle, _ss_ and _ts_
		Normal ns;
		Vector ss, ts;
		if (mesh->n) ns = Normalize(obj2world(b[0] * mesh->n[v[0]] +
			b[1] * mesh->n[v[1]] + b[2] * mesh->n[v[2]]));
		else   ns = dg.nn;
		if (mesh->s) ss = Normalize(obj2world(b[0] * mesh->s[v[0]] +
			b[1] * mesh->s[v[1]] + b[2] * mesh->s[v[2]]));
		else   ss = Normalize(dg.dpdu);
		ts = Normalize(Cross(ss, ns));
		ss = Cross(ts, ns);
		Vector dndu, dndv;
		if (mesh->n) {
			// Compute \dndu and \dndv for triangle shading geometry
			float uvs[3][2];
			GetUVs(uvs);
			// Compute deltas for triangle partial derivatives of normal
			float du1 = uvs[0][0] - uvs[2][0];
			float du2 = uvs[1][0] - uvs[2][0];
			float dv1 = uvs[0][1] - uvs[2][1];
			float dv2 = uvs[1][1] - uvs[2][1];
			Vector dn1 = Vector(mesh->n[v[0]] - mesh->n[v[2]]);
			Vector dn2 = Vector(mesh->n[v[1]] - mesh->n[v[2]]);
			float determinant = du1 * dv2 - dv1 * du2;
			if (determinant == 0)
				dndu = dndv = Vector(0,0,0);
			else {
				float invdet = 1.f / determinant;
				dndu = ( dv2 * dn1 - dv1 * dn2) * invdet;
				dndv = (-du2 * dn1 + du1 * dn2) * invdet;
			}
		}
		else
			dndu = dndv = Vector(0,0,0);
		*dgShading = DifferentialGeometry(dg.p, ss, ts,
			dndu, dndv, dg.u, dg.v, dg.shape);
		dgShading->dudx = dg.dudx;  dgShading->dvdx = dg.dvdx; // NOBOOK
		dgShading->dudy = dg.dudy;  dgShading->dvdy = dg.dvdy; // NOBOOK
		dgShading->dpdx = dg.dpdx;  dgShading->dpdy = dg.dpdy; // NOBOOK
	}
	Point Sample(float u1, float u2, Normal *Ns) const;
private:
	// Triangle Data
	Reference<TriangleMesh> mesh;
	int *v;
};
// TriangleMesh Method Definitions
TriangleMesh::TriangleMesh(const Transform &o2w, bool ro,
		int nt, int nv, const int *vi, const Point *P,
		const Normal *N, const Vector *S, const float *uv)
	: Shape(o2w, ro) {
	ntris = nt;
	nverts = nv;
	vertexIndex = new int[3 * ntris];
	memcpy(vertexIndex, vi, 3 * ntris * sizeof(int));
	// Copy _uv_, _N_, and _S_ vertex data, if present
	if (uv) {
		uvs = new float[2*nverts];
		memcpy(uvs, uv, 2*nverts*sizeof(float));
	}
	else uvs = NULL;
	p = new Point[nverts];
	if (N) {
		n = new Normal[nverts];
		memcpy(n, N, nverts*sizeof(Normal));
	}
	else n = NULL;
	if (S) {
		s = new Vector[nverts];
		memcpy(s, S, nverts*sizeof(Vector));
	}
	else s = NULL;
	// Transform mesh vertices to world space
	for (int i  = 0; i < nverts; ++i)
		p[i] = ObjectToWorld(P[i]);
}
TriangleMesh::~TriangleMesh() {
	delete[] vertexIndex;
	delete[] p;
	delete[] s;
	delete[] n;
	delete[] uvs;
}
BBox TriangleMesh::ObjectBound() const {
	BBox bobj;
	for (int i = 0; i < nverts; i++)
		bobj = Union(bobj, WorldToObject(p[i]));
	return bobj;
}
BBox TriangleMesh::WorldBound() const {
	BBox worldBounds;
	for (int i = 0; i < nverts; i++)
		worldBounds = Union(worldBounds, p[i]);
	return worldBounds;
}
void
TriangleMesh::Refine(vector<Reference<Shape> > &refined)
const {
	for (int i = 0; i < ntris; ++i)
		refined.push_back(new Triangle(ObjectToWorld,
		                               reverseOrientation,
                                       (TriangleMesh *)this,
									   i));
}
BBox Triangle::ObjectBound() const {
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	return Union(BBox(WorldToObject(p1), WorldToObject(p2)),
		WorldToObject(p3));
}
BBox Triangle::WorldBound() const {
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	return Union(BBox(p1, p2), p3);
}
bool Triangle::Intersect(const Ray &ray, float *tHit,
		DifferentialGeometry *dg) const {
	// Initialize triangle intersection statistics
	static
	StatsPercentage triangleHits("Geometry",
	                             "Triangle Ray Intersections");
	// Update triangle tests count
	triangleHits.Add(0, 1);
	// Compute $\VEC{s}_1$
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	Vector e1 = p2 - p1;
	Vector e2 = p3 - p1;
	Vector s1 = Cross(ray.d, e2);
	float divisor = Dot(s1, e1);
	if (divisor == 0.)
		return false;
	float invDivisor = 1.f / divisor;
	// Compute first barycentric coordinate
	Vector d = ray.o - p1;
	float b1 = Dot(d, s1) * invDivisor;
	if (b1 < 0. || b1 > 1.)
		return false;
	// Compute second barycentric coordinate
	Vector s2 = Cross(d, e1);
	float b2 = Dot(ray.d, s2) * invDivisor;
	if (b2 < 0. || b1 + b2 > 1.)
		return false;
	// Compute _t_ to intersection point
	float t = Dot(e2, s2) * invDivisor;
	if (t < ray.mint || t > ray.maxt)
		return false;
	triangleHits.Add(1, 0); //NOBOOK
	// Fill in _DifferentialGeometry_ from triangle hit
	// Compute triangle partial derivatives
	Vector dpdu, dpdv;
	float uvs[3][2];
	GetUVs(uvs);
	// Compute deltas for triangle partial derivatives
	float du1 = uvs[0][0] - uvs[2][0];
	float du2 = uvs[1][0] - uvs[2][0];
	float dv1 = uvs[0][1] - uvs[2][1];
	float dv2 = uvs[1][1] - uvs[2][1];
	Vector dp1 = p1 - p3, dp2 = p2 - p3;
	float determinant = du1 * dv2 - dv1 * du2;
	if (determinant == 0.f) {
		// Handle zero determinant for triangle partial derivative matrix
		CoordinateSystem(Normalize(Cross(e2, e1)), &dpdu, &dpdv);
	}
	else {
		float invdet = 1.f / determinant;
		dpdu = ( dv2 * dp1 - dv1 * dp2) * invdet;
		dpdv = (-du2 * dp1 + du1 * dp2) * invdet;
	}
	// Interpolate $(u,v)$ triangle parametric coordinates
	float b0 = 1 - b1 - b2;
	float tu = b0*uvs[0][0] + b1*uvs[1][0] + b2*uvs[2][0];
	float tv = b0*uvs[0][1] + b1*uvs[1][1] + b2*uvs[2][1];
	*dg = DifferentialGeometry(ray(t), dpdu, dpdv,
	                           Vector(0,0,0), Vector(0,0,0),
							   tu, tv, this);
	*tHit = t;
	return true;
}
bool Triangle::IntersectP(const Ray &ray) const {
	// Initialize triangle intersection statistics
	static
	StatsPercentage triangleHits("Geometry",
	                             "Triangle Ray Intersections");
	// Update triangle tests count
	triangleHits.Add(0, 1);
	// Compute $\VEC{s}_1$
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	Vector e1 = p2 - p1;
	Vector e2 = p3 - p1;
	Vector s1 = Cross(ray.d, e2);
	float divisor = Dot(s1, e1);
	if (divisor == 0.)
		return false;
	float invDivisor = 1.f / divisor;
	// Compute first barycentric coordinate
	Vector d = ray.o - p1;
	float b1 = Dot(d, s1) * invDivisor;
	if (b1 < 0. || b1 > 1.)
		return false;
	// Compute second barycentric coordinate
	Vector s2 = Cross(d, e1);
	float b2 = Dot(ray.d, s2) * invDivisor;
	if (b2 < 0. || b1 + b2 > 1.)
		return false;
	// Compute _t_ to intersection point
	float t = Dot(e2, s2) * invDivisor;
	if (t < ray.mint || t > ray.maxt)
		return false;
	triangleHits.Add(1, 0); //NOBOOK
	return true;
}
void Triangle::GetUVs(float uv[3][2]) const {
	if (mesh->uvs) {
		uv[0][0] = mesh->uvs[2*v[0]];
		uv[0][1] = mesh->uvs[2*v[0]+1];
		uv[1][0] = mesh->uvs[2*v[1]];
		uv[1][1] = mesh->uvs[2*v[1]+1];
		uv[2][0] = mesh->uvs[2*v[2]];
		uv[2][1] = mesh->uvs[2*v[2]+1];
	} else {
		uv[0][0] = 0.; uv[0][1] = 0.;
		uv[1][0] = 1.; uv[1][1] = 0.;
		uv[2][0] = 1.; uv[2][1] = 1.;
	}
}
float Triangle::Area() const {
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	return 0.5f * Cross(p2-p1, p3-p1).Length();
}
Point Triangle::Sample(float u1, float u2,
		Normal *Ns) const {
	float b1, b2;
	UniformSampleTriangle(u1, u2, &b1, &b2);
	// Get triangle vertices in _p1_, _p2_, and _p3_
	const Point &p1 = mesh->p[v[0]];
	const Point &p2 = mesh->p[v[1]];
	const Point &p3 = mesh->p[v[2]];
	Point p = b1 * p1 + b2 * p2 + (1.f - b1 - b2) * p3;
	Normal n = Normal(Cross(p2-p1, p3-p1));
	*Ns = Normalize(n);
	if (reverseOrientation) *Ns *= -1.f;
	return p;
}

class ParametricShape: public Shape {
public:
	// ParametricShape Public Methods
	ParametricShape(const Transform &o2w, bool ro, TriangleMesh *TM,
			Point PMIN, Point PMAX);
	BBox ObjectBound() const;
	bool ParametricShape::CanIntersect() const { return false; }
	void Refine(vector<Reference<Shape> > &refined) const;
	
private:
	// ParametricShape Private Data
	TriangleMesh *tm;
	Point pMin, pMax;
};

// ParametricShape Method Definitions
ParametricShape::ParametricShape(const Transform &o2w, bool ro,
				 TriangleMesh *TM, Point PMIN, Point PMAX)
	: Shape(o2w, ro) {
	tm = TM;
	pMin = PMIN;
	pMax = PMAX;
}

BBox ParametricShape::ObjectBound() const {
	return BBox(pMin, pMax);
}

void ParametricShape::Refine(vector<Reference<Shape> > &refined) const {
    refined.push_back(tm);
}


extern "C" DLLEXPORT
	Shape *CreateShape(const Transform &o2w,
					   bool reverseOrientation,
					   const ParamSet &params) {
//This is where you want to create the triangle mesh and then return a
//new parametric shape
	return new ParametricShape(o2w, reverseOrientation, tm, pMin, pMax);
}