/*
 * C++ RayTracer
 * file: scene.cc
 *
 * Radek Brich, 2006
 */

#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)
{
	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)
{
	//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;
}

BBox Sphere::get_bbox()
{
	BBox bbox = BBox();
	bbox.L = center - radius;
	bbox.H = center + radius;
	return bbox;
}

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

Vector3 Box::normal(Vector3 &P)
{
	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;
}

// this initialization and following intersection methods implements
// Fast Triangle Intersection algorithm from
// http://www.mpi-inf.mpg.de/~wald/PhD/
Triangle::Triangle(const Vector3 &aA, const Vector3 &aB, const Vector3 &aC, Material *amaterial)
	: A(aA), B(aB), C(aC)
{
	material = amaterial;
	Vector3 c = B - A;
	Vector3 b = C - A;

	N = cross(c, b);

	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;
	}

	int u = (k + 1) % 3;
	int v = (k + 2) % 3;

	Float krec = 1.0f / N[k];
	nu = N[u] * krec;
	nv = N[v] * krec;
	nd = dot(N, A) * 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[v] * reci;
	cnv = -c[u] * reci;

	// finalize normal
	N.normalize();
}

// see comment for previous method
bool Triangle::intersect(const Ray &ray, Float &dist)
{
	Vector3 O = ray.o;
	Vector3 D = ray.dir;

	const int modulo3[5] = {0,1,2,0,1};
	const int ku = modulo3[k+1];
	const int kv = modulo3[k+2];
	const Float lnd = 1.0f / (D[k] + nu * D[ku] + nv * D[kv]);
	const Float t = (nd - O[k] - nu * O[ku] - nv * O[kv]) * lnd;

	if (!(t < dist && t > 0))
		return false;

	Float hu = O[ku] + t * D[ku] - A[ku];
	Float hv = O[kv] + t * D[kv] - A[kv];
	Float beta = hv * bnu + hu * bnv;

	if (beta < 0)
		return false;

	Float gamma = hu * cnu + hv * cnv;
	if (gamma < 0)
		return false;

	if ((beta + gamma) > 1)
		return false;

	dist = t;
	return true;
}

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