/*
 * raytracer.cc: Raytracer class
 *
 * This file is part of Pyrit Ray Tracer.
 *
 * Copyright 2006, 2007, 2008  Radek Brich
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

#include <pthread.h>
#include <stdio.h>
#include <stdlib.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;
	if (mat.texture)
		col = mat.texture->evaluate(P);
	else
		col = mat.colour;

	// 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;
	if (mat.texture)
		col = mat.texture->evaluate(P);
	else
		col = mat.colour;

	// 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 + 0.001*newdir, newdir);
					trans_col = raytrace(newray, depth + 1, NULL);
				}
			}
			col = (1-refl-trans)*col + refl*refl_col + trans*trans_col;
		}

		// ambient occlusion
		if (!from_inside && 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;
	}
}

#if 0
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;
}
#endif

void *Raytracer::raytrace_worker(void *d)
{
	static const int my_queue_size = 256;
	Raytracer *rt = (Raytracer*)d;
	Sample my_queue[my_queue_size];
	Colour my_colours[my_queue_size];
	int my_count;
	Ray ray;
	for (;;)
	{
		pthread_mutex_lock(&rt->sample_queue_mutex);
		while (rt->sample_queue_count == 0)
		{
			if (rt->end_of_samples)
			{
				dbgmsg(4, "T thread [%d] exiting\n", pthread_self());
				pthread_mutex_unlock(&rt->sample_queue_mutex);
				pthread_exit(NULL);
			}
			pthread_cond_wait(&rt->sample_queue_cond, &rt->sample_queue_mutex);
		}

		if (rt->sample_queue_count >= my_queue_size)
			my_count = my_queue_size;
		else
			my_count = rt->sample_queue_count;
		rt->sample_queue_count -= my_count;

		// copy samples to local queue
		if (rt->sample_queue_pos + my_count >= rt->sample_queue_size)
		{
			register int c = rt->sample_queue_size - rt->sample_queue_pos;
			memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos, c*sizeof(Sample));
			memcpy(my_queue + c, rt->sample_queue, (my_count - c)*sizeof(Sample));
			rt->sample_queue_pos = my_count - c;
		}
		else
		{
			memcpy(my_queue, rt->sample_queue + rt->sample_queue_pos,
				my_count*sizeof(Sample));
			rt->sample_queue_pos += my_count;
		}
		if (rt->sample_queue_count <= my_queue_size*2)
			pthread_cond_signal(&rt->worker_ready_cond);
		pthread_mutex_unlock(&rt->sample_queue_mutex);

		// do the work
		for (int i = 0; i < my_count; i++)
		{
			ray = rt->camera->makeRay(my_queue[i]);
			my_colours[i] = rt->raytrace(ray, 0, NULL);
		}

		// save the results
		pthread_mutex_lock(&rt->sampler_mutex);
		for (int i = 0; i < my_count; i++)
			rt->sampler->saveSample(my_queue[i], my_colours[i]);
		pthread_mutex_unlock(&rt->sampler_mutex);
	}
}

void Raytracer::render()
{
	if (!sampler || !camera || !top)
		return;

	static const int my_count_max = 256;
	sample_queue_size = my_count_max*2*num_threads;
	sample_queue = new Sample [sample_queue_size];
	sample_queue_pos = 0;
	sample_queue_count = 0;
	end_of_samples = false;
	int my_count, my_pos = 0;
	int sampnum = 0, sampdone;
	bool more_samples;

	sampler->init();

	// create workers
	dbgmsg(1, "* running %d threads\n", num_threads);
	pthread_t threads[num_threads];
	for (int t = 0; t < num_threads; t++)
	{
		int rc = pthread_create(&threads[t], NULL, raytrace_worker, (void*)this);
		if (rc) {
			dbgmsg(0, "\nE pthread_create unsuccessful, return code was %d\n", rc);
			exit(1);
		}
	}

	dbgmsg(1, "* raytracing...\n");
	dbgmsg(2, "-  0%% done");

	pthread_mutex_lock(&sampler_mutex);
	while ( (sampnum = sampler->initSampleSet()) > 0 )
	{
		sampdone = 0;
		for (;;)
		{
			my_count = 0;
			while ( more_samples = sampler->nextSample(&sample_queue[my_pos++]) )
			{
				my_count++;
				if (my_pos >= sample_queue_size)
					my_pos = 0;
				if (my_count >= my_count_max)
					break;
			}
			if (!more_samples && !my_count)
				break;

			pthread_mutex_unlock(&sampler_mutex);
			pthread_mutex_lock(&sample_queue_mutex);
			sample_queue_count += my_count;

			// wait for workers if there is enough samples ready on queue
			while (sample_queue_count > (2*num_threads-1)*my_count_max)
				pthread_cond_wait(&worker_ready_cond, &sample_queue_mutex);

			pthread_cond_signal(&sample_queue_cond);
			pthread_mutex_unlock(&sample_queue_mutex);

			sampdone += my_count;
			dbgmsg(2, "\b\b\b\b\b\b\b\b%2d%% done", (sampdone - sample_queue_count)*100/sampnum);

			pthread_mutex_lock(&sampler_mutex);
		}
		dbgmsg(2, "\b\b\b\b\b\b\b\b100%% done\n");
	}
	pthread_mutex_unlock(&sampler_mutex);

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

	pthread_mutex_lock(&sample_queue_mutex);
	end_of_samples = true;
	while (sample_queue_count)
	{
		pthread_cond_broadcast(&sample_queue_cond);
		pthread_mutex_unlock(&sample_queue_mutex);
		pthread_mutex_lock(&sample_queue_mutex);
	}
	pthread_cond_broadcast(&sample_queue_cond);
	pthread_mutex_unlock(&sample_queue_mutex);

	for (int t = 0; t < num_threads; t++)
		pthread_join(threads[t], NULL);

	delete[] sample_queue;
}

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