/*

 * Pyrit Ray Tracer

 * file: raytracer.cc

 *

 * Radek Brich, 2006-2007

 */

#ifdef PTHREADS

#include <pthread.h>

#endif

#include <stdio.h>

#include <malloc.h>

#include "raytracer.h"

int pyrit_verbosity = 2;

// Hammersley spherical point distribution

// http://www.cse.cuhk.edu.hk/~ttwong/papers/udpoint/udpoints.html

Vector3 Raytracer::SphereDistribute(int i, int n, Float extent, Vector3 &normal)

{

	Float p, t, st, phi, phirad;

	int kk;

	t = 0;

	for (p=0.5, kk=i; kk; p*=0.5, kk>>=1)

	if (kk & 1)

		t += p;

	t = 1.0 + (t - 1.0)*extent;

	phi = (i + 0.5) / n;

	phirad =  phi * 2.0 * M_PI;

	st = sqrt(1.0 - t*t);

	Float x, y, z, xx, yy, zz, q;

	x = st * cos(phirad);

	y = st * sin(phirad);

	z = t;

	// rotate against Y axis

	q = acos(normal.z);

	zz = z*cos(q) - x*sin(q);

	xx = z*sin(q) + x*cos(q);

	yy = y;

	// rotate against Z axis

	q = atan2f(normal.y, normal.x);

	x = xx*cos(q) - yy*sin(q);

	y = xx*sin(q) + yy*cos(q);

	z = zz;

	return Vector3(x, y, z);

}

// ---- tyto dve funkce budou v budouci verzi metody objektu PhongShader

// calculate shader function

// P is point of intersection, N normal in this point

Colour PhongShader_ambient(Material &mat, Vector3 &P)

{

	Colour col = mat.texture.colour; //mat.texture.evaluate(P);

	// ambient

	return mat.ambient * col;

}

/*

 P is point of intersection,

 N normal in this point,

 R direction of reflected ray,

 V direction to the viewer

*/

Colour PhongShader_calculate(Material &mat, Vector3 &P, Vector3 &N, Vector3 &R, Vector3 &V,

	Light &light)

{

	Colour I = Colour();

	Vector3 L = light.pos - P;

	L.normalize();

	Float L_dot_N = dot(L, N);

	Float R_dot_V = dot(R, V);

	Colour col = mat.texture.colour; //mat.texture.evaluate(P);

	// diffuse

	I = mat.diffuse * col * light.colour * L_dot_N;

	// specular

	if (R_dot_V > 0)

		I += mat.specular * light.colour * powf(R_dot_V, mat.shininess);

	return I;

}

Colour Raytracer::raytrace(Ray &ray, int depth, Shape *origin_shape)

{

	Float nearest_distance = Inf;

	Shape *nearest_shape = top->nearest_intersection(origin_shape, ray, nearest_distance);

	if (nearest_shape == NULL) {

		return bg_colour;

	} else {

		Colour col = Colour();

		Vector3 P = ray.o + ray.dir * nearest_distance; // point of intersection

		Vector3 normal = nearest_shape->normal(P);

		bool from_inside = false;

		// make shapes double sided

		if (dot(normal, ray.dir) > 0.0)

		{

			normal = - normal;

			from_inside = true;

		}

		col = PhongShader_ambient(*nearest_shape->material, P);

		vector<Light*>::iterator light;

		for (light = lights.begin(); light != lights.end(); light++) {

			Vector3 jo, L = (*light)->pos - P; // direction vector to light

			L.normalize();

			Float L_dot_N = dot(L, normal);

			if (L_dot_N > 0) {

				// test if this light is occluded (sharp shadows)

				if ((*light)->cast_shadows) {

					Ray shadow_ray = Ray(P, L);

					Float dist = FLT_MAX;

					if (top->nearest_intersection(nearest_shape, shadow_ray, dist))

						continue;

				}

				// shading function

				Vector3 R = L - 2.0 * L_dot_N * normal;

				col += PhongShader_calculate(*nearest_shape->material,

					P, normal, R, ray.dir, **light);

			}

		}

		if (depth < max_depth)

		{

			Colour trans_col, refl_col;

			Float trans = nearest_shape->material->transmissivity;

			Float refl = nearest_shape->material->reflectivity;

			const Float cos_i = - dot(normal, ray.dir);

			// reflection

			if (refl > 0.01)

			{

				Vector3 newdir = ray.dir + 2.0 * cos_i * normal;

				Ray newray = Ray(P, newdir);

				refl_col = raytrace(newray, depth + 1, nearest_shape);

			}

			// refraction

		 	if (trans > 0.01)

			{

				Float n, n1, n2;

				if (from_inside)

				{

					n1 = nearest_shape->material->refract_index;

					n2 = 1.0;

					n = n1;

				}

				else

				{

					n1 = 1.0;

					n2 = nearest_shape->material->refract_index;

					n = 1.0 / n2;

				}

				const Float sin2_t = n*n * (1 - cos_i*cos_i);

				if (sin2_t >= 1.0)

				{

					 // totally reflected

					 refl += trans;

					 trans = 0;

				}

				else

				{

					const Float cos_t = sqrtf(1 - sin2_t);

					const Float Rdiv = 1.0/(n1*cos_i + n2*cos_t);

					const Float Rper = (n1*cos_i - n2*cos_t)*Rdiv;

					const Float Rpar = (n2*cos_i - n1*cos_t)*Rdiv;

					const Float R = (Rper*Rper + Rpar*Rpar)/2;

					refl += R*trans;

					trans = (1-R)*trans;

					Vector3 newdir = n * ray.dir + (n*cos_i - cos_t) * normal;

					Ray newray = Ray(P, newdir);

					trans_col = raytrace(newray, depth + 1, nearest_shape);

				}

			}

			col = (1-refl-trans)*col + refl*refl_col + trans*trans_col;

		}

		// ambient occlusion

		if (ao_samples)

		{

			Float miss = 0;

			for (int i = 0; i < ao_samples; i++) {

				Vector3 dir = SphereDistribute(i, ao_samples, ao_angle, normal);

				Ray ao_ray = Ray(P, dir);

				Float dist = ao_distance;

				Shape *shape_in_way = top->nearest_intersection(nearest_shape, ao_ray, dist);

				if (shape_in_way == NULL)

					miss += 1.0;

				else

					miss += dist / ao_distance;

			}

			Float ao_intensity = miss / ao_samples;

			col = col * ao_intensity;

		}

		return col;

	}

}

static inline void samplepixel(Colour &c, Vector3 &dir, RenderrowData* d, const int &oversample)

{

	if (oversample <= 0)

	{

		// no oversampling

		dir.normalize();

		Ray ray(d->eye, dir);

		c = d->rt->raytrace(ray, 0, NULL);

	}

	else

	if (oversample <= 3)

	{

		// grid oversampling

		c = Colour(0,0,0);

		static const int gridsamples[] = {5,9,16};

		static const Float osa5x[] = {0.0, -0.4, +0.4, +0.4, -0.4};

		static const Float osa5y[] = {0.0, -0.4, -0.4, +0.4, +0.4};

		static const Float osa9x[] = {-0.34,  0.00, +0.34,

			-0.34,  0.00, +0.34, -0.34,  0.00, +0.34};

		static const Float osa9y[] = {-0.34, -0.34, -0.34,

				0.00,  0.00,  0.00, +0.34, +0.34, +0.34};

		static const Float osa16x[] = {-0.375, -0.125, +0.125, +0.375,

			-0.375, -0.125, +0.125, +0.375, -0.375, -0.125, +0.125, +0.375,

			-0.375, -0.125, +0.125, +0.375};

		static const Float osa16y[] = {-0.375, -0.375, -0.375, -0.375,

			-0.125, -0.125, -0.125, -0.125, +0.125, +0.125, +0.125, +0.125,

			+0.375, +0.375, +0.375, +0.375};

		static const Float *osaSx[] = {osa5x, osa9x, osa16x};

		static const Float *osaSy[] = {osa5y, osa9y, osa16y};

		const int samples = gridsamples[oversample-1];

		const Float *osax = osaSx[oversample-1];

		const Float *osay = osaSy[oversample-1];

		for (int i = 0; i < samples; i++)

		{

			Vector3 tmpdir = dir + osax[i]*d->dx + osay[i]*d->dy;

			tmpdir.normalize();

			Ray ray(d->eye, tmpdir);

			c += d->rt->raytrace(ray, 0, NULL);

		}

		c = c * (1.0/samples);

	}

	else

	{

		// stochastic oversampling

		// ...todo

	}

}

static void *renderrow(void *data)

{

	RenderrowData *d = (RenderrowData*) data;

	const int subsample = d->rt->getSubsample();

	const Float subsample2 = 1.0/(subsample*subsample);

	const int oversample = d->rt->getOversample();

	const int ww = d->w*3;

	Vector3 dir = d->dfix;

	for (int x = 0; x < d->w; x += subsample) {

		// generate a ray from eye passing through this pixel

		if (subsample > 1)

		{

			Colour ic;

			// top-left

			dir.normalize();

			Ray ray(d->eye, dir);

			Colour c1 = d->rt->raytrace(ray, 0, NULL);

			// top-right

			Vector3 tmpdir = dir + (subsample-1)*d->dx;

			tmpdir.normalize();

			ray.dir = tmpdir;

			Colour c2 = d->rt->raytrace(ray, 0, NULL);

			// bottom right

			tmpdir += (subsample-1)*d->dy;

			tmpdir.normalize();

			ray.dir = tmpdir;

			Colour c4 = d->rt->raytrace(ray, 0, NULL);

			// bottom left

			tmpdir = dir + (subsample-1)*d->dy;

			tmpdir.normalize();

			ray.dir = tmpdir;

			Colour c3 = d->rt->raytrace(ray, 0, NULL);

			// are the colors similar?

			Float m = (c1-c2).mag2();

			m = max(m, (c2-c3).mag2());

			m = max(m, (c3-c4).mag2());

			m = max(m, (c4-c1).mag2());

			if (m < 0.001)

			{

				// interpolate

				Float *i = d->iter;

				for (int x = 0; x < subsample; x++)

				{

					for (int y = 0; y < subsample; y++)

					{

						ic = c1*(subsample-x)*(subsample-y)*subsample2

							+ c2*(x)*(subsample-y)*subsample2

							+ c3*(subsample-x)*(y)*subsample2

							+ c4*(x)*(y)*subsample2;

						*(i + ww*y) = ic.r;

						*(i + ww*y + 1) = ic.g;

						*(i + ww*y + 2) = ic.b;

					}

					i += 3;

				}

				d->iter = i;

			}

			else

			{

				// render all pixels

				Vector3 tmpdir = dir;

				if (oversample)

				{

					for (int x = 0; x < subsample; x++)

					{

						for (int y = 0; y < subsample; y++)

						{

							Vector3 tmp2dir = tmpdir + y*d->dy;

							samplepixel(ic, tmp2dir, d, oversample);

							*(d->iter + ww*y) = ic.r;

							*(d->iter + ww*y + 1) = ic.g;

							*(d->iter + ww*y + 2) = ic.b;

						}

						d->iter += 3;

						tmpdir += d->dx;

					}

				}

				else

				{

					/* this is so complex because it tries to reuse

					   already computed corner pixels

					   though, above code will also work for non-oversampling... */

					// first column

					*(d->iter) = c1.r;

					*(d->iter + 1) = c1.g;

					*(d->iter + 2) = c1.b;

					for (int y = 1; y < subsample-1; y++)

					{

						Vector3 tmp2dir = tmpdir + y*d->dy;

						tmp2dir.normalize();

						ray.dir = tmp2dir;

						ic = d->rt->raytrace(ray, 0, NULL);

						*(d->iter + ww*y) = ic.r;

						*(d->iter + ww*y + 1) = ic.g;

						*(d->iter + ww*y + 2) = ic.b;

					}

					*(d->iter + ww*(subsample-1)) = c3.r;

					*(d->iter + ww*(subsample-1) + 1) = c3.g;

					*(d->iter + ww*(subsample-1) + 2) = c3.b;

					d->iter += 3;

					tmpdir += d->dx;

					// middle

					for (int x = 1; x < subsample-1; x++)

					{

						for (int y = 0; y < subsample; y++)

						{

							Vector3 tmp2dir = tmpdir + y*d->dy;

							tmp2dir.normalize();

							ray.dir = tmp2dir;

							ic = d->rt->raytrace(ray, 0, NULL);

							*(d->iter + ww*y) = ic.r;

							*(d->iter + ww*y + 1) = ic.g;

							*(d->iter + ww*y + 2) = ic.b;

						}

						d->iter += 3;

						tmpdir += d->dx;

					}

					// last column

					*(d->iter) = c2.r;

					*(d->iter + 1) = c2.g;

					*(d->iter + 2) = c2.b;

					for (int y = 1; y < subsample-1; y++)

					{

						Vector3 tmp2dir = tmpdir + y*d->dy;

						tmp2dir.normalize();

						ray.dir = tmp2dir;

						ic = d->rt->raytrace(ray, 0, NULL);

						*(d->iter + ww*y) = ic.r;

						*(d->iter + ww*y + 1) = ic.g;

						*(d->iter + ww*y + 2) = ic.b;

					}

					*(d->iter + ww*(subsample-1)) = c4.r;

					*(d->iter + ww*(subsample-1) + 1) = c4.g;

					*(d->iter + ww*(subsample-1) + 2) = c4.b;

					d->iter += 3;

				}

			}

		}

		else // subsample <= 1

		{

			Colour c;

			samplepixel(c, dir, d, oversample);

			// write color to buffer

			*d->iter++ = c.r;

			*d->iter++ = c.g;

			*d->iter++ = c.b;

		}

		dir += d->dx*subsample;

	}

#ifdef PTHREADS

	pthread_exit((void *)d);

#endif

	return (void *)d;

}

void Raytracer::render(int w, int h, Float *buffer)

{

	if (!camera || !top || !buffer)

		return;

	RenderrowData *d;

	Float S = 0.5/w;

	Vector3 dfix = camera->u*(-w/2.0*S/camera->f)

		+ camera->v*(h/2.0*S/camera->f) + camera->p;

	Vector3 dx = camera->u * (S/camera->f);

	Vector3 dy = camera->v * (-S/camera->f);

#ifdef PTHREADS

	dbgmsg(1, "* pthreads enabled, using %d threads\n", num_threads);

	pthread_t threads[num_threads];

	for (int t = 0; t < num_threads; t++)

		threads[t] = pthread_self();

	int t = 0;

#endif

	/* for each pixel... */

	dbgmsg(1, "* raytracing...\n");

	dbgmsg(2, "- row   0 (  0%% done)");

	for (int y = 0; y < h; y += subsample)

	{

		d = (RenderrowData*) malloc(sizeof(RenderrowData));

		d->rt = this;

		d->w = w;

		d->eye = camera->eye;

		d->dfix = dfix;

		d->dx = dx;

		d->dy = dy;

		d->iter = buffer + y*3*w;

#ifdef PTHREADS

		/* create new thread and increase 't' */

		int rc = pthread_create(&threads[t++], NULL, renderrow, (void *)d);

		if (rc) {

			dbgmsg(0, "\nE pthread_create unsuccessful, return code was %d\n", rc);

			exit(1);

		}

		/* when 't' owerflows, reset it */

		if (t >= num_threads)

			t = 0;

		/* wait for next thread in fifo queue, so the descriptor can be reused;

		   this also limits number of running threads */

		if (!pthread_equal(threads[t], pthread_self()))

			if (pthread_join(threads[t], (void**)&d) == 0)

				free(d);

#else

		renderrow((void *)d);

		free(d);

#endif

		dfix += dy*subsample;

		dbgmsg(2, "\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%4d (%2d%% done)", y, y*100/(h-1));

	}

	dbgmsg(2, "\n");

#ifdef PTHREADS

	dbgmsg(2, "- waiting for threads to finish\n");

	for (t = 0; t < num_threads; t++)

		if (pthread_join(threads[t], (void**)&d) == 0)

			free(d);

#endif

}

void Raytracer::addlight(Light *light)

{

	lights.push_back(light);

}

void Raytracer::ambientocclusion(int samples, Float distance, Float angle)

{

	ao_samples = samples;

	ao_distance = distance;

	ao_angle = angle;

	if (ao_distance == 0)

		/* 0 ==> Inf */

		ao_distance = FLT_MAX;