/*

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

	}

}

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_count = 0;

	int sampnum = 0, sampdone;

	bool more_samples;

	sampler->init();

	// create workers

	pthread_t threads[num_threads];

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

	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_cond_broadcast(&sample_queue_cond);

		pthread_mutex_unlock(&sample_queue_mutex);

	}

	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;

}