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

#ifdef PTHREADS
#include <pthread.h>
#endif

#include <stdio.h>
#include <malloc.h>
#include "raytracer.h"

// 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 acc = Colour();
		Vector3 P = ray.o + ray.dir * nearest_distance; // point of intersection
		Vector3 normal = nearest_shape->normal(P);
		acc = 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;
				acc += PhongShader_calculate(*nearest_shape->material,
					P, normal, R, ray.dir, **light);
			}
		}

		// reflection
		Vector3 newdir = ray.dir - 2.0 * dot(ray.dir, normal) * normal;
		if (depth < max_depth && nearest_shape->material->reflection > 0.01) {
			Ray newray = Ray(P, newdir);
			Colour refl_col = raytrace(newray, depth + 1, nearest_shape);
			acc += nearest_shape->material->reflection * refl_col;
		}

		// refraction
		/* ... */

		// 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;
			acc = acc * ao_intensity;
		}

		return acc;
	}
}

static void *renderrow(void *data)
{
	RenderrowData *d = (RenderrowData*) data;
	int subsample = d->rt->getSubsample();
	Float subsample2 = 1.0/(subsample*subsample);
	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 OVERSAMPLING
		// 5x oversampling
		Colour c = Colour();

		for (int i = 0; i < 5; i++)
		{
			Float osax[] = {0.0, -0.4, +0.4, +0.4, -0.4};
			Float osay[] = {0.0, -0.4, -0.4, +0.4, +0.4};
			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./5);
#else
		// no oversampling
		dir.normalize();
		Ray ray(d->eye, dir);
		Colour c = d->rt->raytrace(ray, 0, NULL);
		if (subsample > 1)
		{
			Colour ic;
			// top-left is 'c'
			// top-right
			Vector3 tmpdir = dir + (subsample-1)*d->dx;
			tmpdir.normalize();
			Ray ray(d->eye, 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 = (c-c2).mag2();
			m = max(m, (c2-c3).mag2());
			m = max(m, (c3-c4).mag2());
			m = max(m, (c4-c).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 = c*(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
				Vector3 tmpdir = dir;
				// first column
				*(d->iter) = c.r;
				*(d->iter + 1) = c.g;
				*(d->iter + 2) = c.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;
			}
		}
#endif
		if (subsample <= 1)
		{
			*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
	infomsg("* 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... */
	infomsg("* rendering 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) {
			infomsg("\nERROR: return code from pthread_create() is %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;
		infomsg("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%4d (%2d%% done)", y, y*100/(h-1));
	}
	infomsg("\n");

#ifdef PTHREADS
	infomsg("* 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;
}