/*
 * Pyrit Ray Tracer
 * file: scene.cc
 *
 * Radek Brich, 2006-2007
 */

#include <math.h>

#include "common.h"
#include "scene.h"

void Camera::rotate(const Quaternion &q)
{
	/*
	//non-optimized
	Quaternion res;
	res = q * Quaternion(u) * conjugate(q);
	u = res.toVector();
	res = q * Quaternion(v) * conjugate(q);
	v = res.toVector();
	res = q * Quaternion(p) * conjugate(q);
	p = res.toVector();
	*/

	// optimized
	Float t2 =   q.a*q.b;
	Float t3 =   q.a*q.c;
	Float t4 =   q.a*q.d;
	Float t5 =  -q.b*q.b;
	Float t6 =   q.b*q.c;
	Float t7 =   q.b*q.d;
	Float t8 =  -q.c*q.c;
	Float t9 =   q.c*q.d;
	Float t10 = -q.d*q.d;
	Float x,y,z;
	x = 2*( (t8 + t10)*p.x + (t6 -  t4)*p.y + (t3 + t7)*p.z ) + p.x;
	y = 2*( (t4 +  t6)*p.x + (t5 + t10)*p.y + (t9 - t2)*p.z ) + p.y;
	z = 2*( (t7 -  t3)*p.x + (t2 +  t9)*p.y + (t5 + t8)*p.z ) + p.z;
	p = Vector3(x,y,z);
	x = 2*( (t8 + t10)*u.x + (t6 -  t4)*u.y + (t3 + t7)*u.z ) + u.x;
	y = 2*( (t4 +  t6)*u.x + (t5 + t10)*u.y + (t9 - t2)*u.z ) + u.y;
	z = 2*( (t7 -  t3)*u.x + (t2 +  t9)*u.y + (t5 + t8)*u.z ) + u.z;
	u = Vector3(x,y,z);
	x = 2*( (t8 + t10)*v.x + (t6 -  t4)*v.y + (t3 + t7)*v.z ) + v.x;
	y = 2*( (t4 +  t6)*v.x + (t5 + t10)*v.y + (t9 - t2)*v.z ) + v.y;
	z = 2*( (t7 -  t3)*v.x + (t2 +  t9)*v.y + (t5 + t8)*v.z ) + v.z;
	v = Vector3(x,y,z);
	p.normalize();
	u.normalize();
	v.normalize();
}

void Camera::move(const Float fw, const Float left, const Float up)
{
	eye = eye + fw*p + left*u + up*v;
}

/* http://www.siggraph.org/education/materials/HyperGraph/raytrace/rtinter3.htm */
bool BBox::intersect(const Ray &ray, Float &a, Float &b)
{
	Float tnear = -FLT_MAX;
	Float tfar = FLT_MAX;
	Float t1, t2;

	for (int i = 0; i < 3; i++)
	{
		if (ray.dir.cell[i] == 0) {
			/* ray is parallel to these planes */
			if (ray.o.cell[i] < L.cell[i] || ray.o.cell[i] > H.cell[i])
				return false;
		} else
		{
			/* compute the intersection distance of the planes */
			t1 = (L.cell[i] - ray.o.cell[i]) / ray.dir.cell[i];
			t2 = (H.cell[i] - ray.o.cell[i]) / ray.dir.cell[i];

			if (t1 > t2)
				swap(t1, t2);

			if (t1 > tnear)
				tnear = t1; /* want largest Tnear */
			if (t2 < tfar)
				tfar = t2; /* want smallest Tfar */
			if (tnear > tfar)
				return false; /* box missed */
			if (tfar < 0)
				return false; /* box is behind ray */
		}
	}

	a = tnear;
	b = tfar;
	return true;
}

bool Sphere::intersect(const Ray &ray, Float &dist) const
{
	Vector3 V = ray.o - center;
	register Float d = -dot(V, ray.dir);
	register Float Det = d * d - (dot(V,V) - sqr_radius);
	if (Det > 0) {
		d -= sqrtf(Det);
		if (d > 0 && d < dist)
		{
			dist = d;
			return true;
		}
	}
	return false;
}

bool Sphere::intersect_all(const Ray &ray, Float dist, vector<Float> &allts) const
{
	//allts = new vector<Float>();

	Vector3 V = ((Ray)ray).o - center;
	Float Vd = - dot(V, ray.dir);
	Float Det = Vd * Vd - (dot(V,V) - sqr_radius);

	if (Det > 0) {
		Det = sqrtf(Det);
		Float t1 = Vd - Det;
		Float t2 = Vd + Det;
		if (t1 < 0)
		{
			if (t2 > 0)
			{
				// ray from inside of the sphere
				allts.push_back(0.0);
				allts.push_back(t2);
				return true;
			}
			else
				return false;
		}
		else
		{
			allts.push_back(t1);
			allts.push_back(t2);
			return true;
		}
	}
	return false;
}

bool Sphere::intersect_bbox(const BBox &bbox) const
{
	register float dmin = 0;
	for (int i = 0; i < 3; i++)
	{
		if (center[i] < bbox.L[i])
			dmin += (center[i] - bbox.L[i])*(center[i] - bbox.L[i]);
		else
		if (center[i] > bbox.H[i])
			dmin += (center[i] - bbox.H[i])*(center[i] - bbox.H[i]);
	}
	if (dmin <= sqr_radius)
		return true;
	return false;
};

BBox Sphere::get_bbox() const
{
	return BBox(center - radius, center + radius);
}

bool Box::intersect(const Ray &ray, Float &dist) const
{
	Float a,b;
	if (get_bbox().intersect(ray, a, b) && a < dist)
	{
		dist = a;
		return true;
	}
	else
		return false;
}

bool Box::intersect_bbox(const BBox &bbox) const
{
	return (
	H.x > bbox.L.x && L.x < bbox.H.x &&
	H.y > bbox.L.y && L.y < bbox.H.y &&
	H.z > bbox.L.z && L.z < bbox.H.z);
}

const Vector3 Box::normal(const Vector3 &P) const
{
	Vector3 N;
	for (int i = 0; i < 3; i++)
	{
		if (P.cell[i] >= L.cell[i]-Eps && P.cell[i] <= L.cell[i]+Eps)
		//if (P.cell[i] == L.cell[i])
		{
			N.cell[i] = -1.0;
			break;
		}
		if (P.cell[i] >= H.cell[i]-Eps && P.cell[i] <= H.cell[i]+Eps)
		//if (P.cell[i] == H.cell[i])
		{
			N.cell[i] = +1.0;
			break;
		}
	}
	return N;
}

#ifdef TRI_PLUCKER
inline void Plucker(const Vector3 &p, const Vector3 &q, Float* pl)
{
    pl[0] = p.x*q.y - q.x*p.y;
    pl[1] = p.x*q.z - q.x*p.z;
    pl[2] = p.x - q.x;
    pl[3] = p.y*q.z - q.y*p.z;
    pl[4] = p.z - q.z;
    pl[5] = q.y - p.y;
}

inline Float Side(const Float* pla, const Float* plb)
{
    return pla[0]*plb[4] + pla[1]*plb[5] + pla[2]*plb[3] + pla[4]*plb[0] + pla[5]*plb[1] + pla[3]*plb[2];
}
#endif

Triangle::Triangle(Vertex *aA, Vertex *aB, Vertex *aC, Material *amaterial)
	: smooth(false), A(aA), B(aB), C(aC)
{
	material = amaterial;

	const Vector3 c = B->P - A->P;
	const Vector3 b = C->P - A->P;

	N = cross(c, b);
	N.normalize();

#ifdef TRI_PLUCKER
	Plucker(B->P,C->P,pla);
	Plucker(C->P,A->P,plb);
	Plucker(A->P,B->P,plc);
#endif

#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
	if (fabsf(N.x) > fabsf(N.y))
	{
		if (fabsf(N.x) > fabsf(N.z))
			k = 0;
		else
			k = 2;
	}
	else
	{
		if (fabsf(N.y) > fabsf(N.z))
			k = 1;
		else
			k = 2;
	}
#endif
#ifdef TRI_BARI_PRE
	int u = (k + 1) % 3;
	int v = (k + 2) % 3;

	Float krec = 1.0 / N[k];
	nu = N[u] * krec;
	nv = N[v] * krec;
	nd = dot(N, A->P) * krec;

	// first line equation
	Float reci = 1.0f / (b[u] * c[v] - b[v] * c[u]);
	bnu = b[u] * reci;
	bnv = -b[v] * reci;

	// second line equation
	cnu = -c[u] * reci;
	cnv = c[v] * reci;
#endif
}

bool Triangle::intersect(const Ray &ray, Float &dist) const
{
#ifdef TRI_PLUCKER
	Float plr[6];
	Plucker(ray.o, ray.o+ray.dir, plr);
	const bool side0 = Side(plr, pla) >= 0.0;
	const bool side1 = Side(plr, plb) >= 0.0;
	if (side0 != side1)
		return false;
	const bool side2 = Side(plr, plc) >= 0.0;
	if (side0 != side2)
		return false;
	const Float t = - dot( (ray.o - A->P), N) / dot(ray.dir,N);
	if(t <= Eps || t >= dist)
		return false;
	dist = t;
	return true;
#endif

#if defined(TRI_BARI) || defined(TRI_BARI_PRE)
	static const int modulo3[5] = {0,1,2,0,1};
	const Vector3 &O = ray.o;
	const Vector3 &D = ray.dir;
	register const int u = modulo3[k+1];
	register const int v = modulo3[k+2];
#endif
#ifdef TRI_BARI_PRE
	const Float t = (nd - O[k] - nu * O[u] - nv * O[v]) / (D[k] + nu * D[u] + nv * D[v]);

	if (t >= dist || t < Eps)
		return false;

	const Float hu = O[u] + t * D[u] - A->P[u];
	const Float hv = O[v] + t * D[v] - A->P[v];
	const Float beta = hv * bnu + hu * bnv;

	if (beta < 0.)
		return false;

	const Float gamma = hu * cnv + hv * cnu;
	if (gamma < 0. || beta + gamma > 1.)
		return false;

	dist = t;
	return true;
#endif

#ifdef TRI_BARI
	// original barycentric coordinates based intesection
	// not optimized, just for reference
	const Vector3 c = B - A;
	const Vector3 b = C - A;
	// distance test
	const Float t = - dot( (O-A), N) / dot(D,N);
	if (t < Eps || t > dist)
		return false;

	// calc hitpoint
	const Float Hu = O[u] + t * D[u] - A[u];
	const Float Hv = O[v] + t * D[v] - A[v];
	const Float beta = (b[u] * Hv - b[v] * Hu) / (b[u] * c[v] - b[v] * c[u]);
	if (beta < 0)
		return false;
	const Float gamma = (c[v] * Hu - c[u] * Hv) / (b[u] * c[v] - b[v] * c[u]);
	if (gamma < 0)
		return false;
	if (beta+gamma > 1)
		return false;
	dist = t;
	return true;
#endif
}

bool Triangle::intersect_bbox(const BBox &bbox) const
{
	const Vector3 boxcenter = (bbox.L+bbox.H)*0.5;
	const Vector3 boxhalfsize = (bbox.H-bbox.L)*0.5;
	const Vector3 v0 = A->P - boxcenter;
	const Vector3 v1 = B->P - boxcenter;
	const Vector3 v2 = C->P - boxcenter;
	const Vector3 e0 = v1-v0;
	const Vector3 e1 = v2-v1;
	const Vector3 e2 = v0-v2;

	Float fex = fabsf(e0.x);
	Float fey = fabsf(e0.y);
	Float fez = fabsf(e0.z);

	Float p0,p1,p2,min,max,rad;

	p0 = e0.z*v0.y - e0.y*v0.z;
	p2 = e0.z*v2.y - e0.y*v2.z;
	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p0 = -e0.z*v0.x + e0.x*v0.z;
	p2 = -e0.z*v2.x + e0.x*v2.z;
	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p1 = e0.y*v1.x - e0.x*v1.y;
	p2 = e0.y*v2.x - e0.x*v2.y;
	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
	if(min>rad || max<-rad) return false;

	fex = fabsf(e1.x);
	fey = fabsf(e1.y);
	fez = fabsf(e1.z);

	p0 = e1.z*v0.y - e1.y*v0.z;
	p2 = e1.z*v2.y - e1.y*v2.z;
	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p0 = -e1.z*v0.x + e1.x*v0.z;
	p2 = -e1.z*v2.x + e1.x*v2.z;
	if(p0<p2) {min=p0; max=p2;} else {min=p2; max=p0;}
	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p0 = e1.y*v0.x - e1.x*v0.y;
	p1 = e1.y*v1.x - e1.x*v1.y;
	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
	if(min>rad || max<-rad) return false;

	fex = fabsf(e2.x);
	fey = fabsf(e2.y);
	fez = fabsf(e2.z);

	p0 = e2.z*v0.y - e2.y*v0.z;
	p1 = e2.z*v1.y - e2.y*v1.z;
	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
	rad = fez * boxhalfsize.y + fey * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p0 = -e2.z*v0.x + e2.x*v0.z;
	p1 = -e2.z*v1.x + e2.x*v1.z;
	if(p0<p1) {min=p0; max=p1;} else {min=p1; max=p0;}
	rad = fez * boxhalfsize.x + fex * boxhalfsize.z;
	if(min>rad || max<-rad) return false;

	p1 = e2.y*v1.x - e2.x*v1.y;
	p2 = e2.y*v2.x - e2.x*v2.y;
	if(p2<p1) {min=p2; max=p1;} else {min=p1; max=p2;}
	rad = fey * boxhalfsize.x + fex * boxhalfsize.y;
	if(min>rad || max<-rad) return false;

	/* test overlap in the {x,y,z}-directions */
	/* test in X-direction */
	min = v0.x;
	if (v1.x < min) min = v1.x;
	if (v2.x < min) min = v2.x;
	max = v0.x;
	if (v1.x > max) max = v1.x;
	if (v2.x > max) max = v2.x;
	if(min>boxhalfsize.x || max<-boxhalfsize.x) return false;

	/* test in Y-direction */
	min = v0.y;
	if (v1.y < min) min = v1.y;
	if (v2.y < min) min = v2.y;
	max = v0.y;
	if (v1.y > max) max = v1.y;
	if (v2.y > max) max = v2.y;
	if(min>boxhalfsize.y || max<-boxhalfsize.y) return false;

	/* test in Z-direction */
	min = v0.z;
	if (v1.z < min) min = v1.z;
	if (v2.z < min) min = v2.z;
	max = v0.z;
	if (v1.z > max) max = v1.z;
	if (v2.z > max) max = v2.z;
	if(min>boxhalfsize.z || max<-boxhalfsize.z) return false;

	/*  test if the box intersects the plane of the triangle */
	Vector3 vmin,vmax;
	Float v;
	for(int q=0;q<3;q++)
	{
		v=v0[q];
		if(N[q]>0.0f)
		{
			vmin.cell[q]=-boxhalfsize[q] - v;
			vmax.cell[q]= boxhalfsize[q] - v;
		}
		else
		{
			vmin.cell[q]= boxhalfsize[q] - v;
			vmax.cell[q]=-boxhalfsize[q] - v;
		}
	}
	if(dot(N,vmin)>0.0f) return false;
	if(dot(N,vmax)>=0.0f) return true;

	return false;
}

BBox Triangle::get_bbox() const
{
	BBox bbox = BBox();
	bbox.L = A->P;
	if (B->P.x < bbox.L.x)  bbox.L.x = B->P.x;
	if (C->P.x < bbox.L.x)  bbox.L.x = C->P.x;
	if (B->P.y < bbox.L.y)  bbox.L.y = B->P.y;
	if (C->P.y < bbox.L.y)  bbox.L.y = C->P.y;
	if (B->P.z < bbox.L.z)  bbox.L.z = B->P.z;
	if (C->P.z < bbox.L.z)  bbox.L.z = C->P.z;
	bbox.H = A->P;
	if (B->P.x > bbox.H.x)  bbox.H.x = B->P.x;
	if (C->P.x > bbox.H.x)  bbox.H.x = C->P.x;
	if (B->P.y > bbox.H.y)  bbox.H.y = B->P.y;
	if (C->P.y > bbox.H.y)  bbox.H.y = C->P.y;
	if (B->P.z > bbox.H.z)  bbox.H.z = B->P.z;
	if (C->P.z > bbox.H.z)  bbox.H.z = C->P.z;
	return bbox;
};