//-*-c-*-
#include <math_constants.h>

#include "linalg.h"
#include "matrix.h"
#include "rotate.h"

texture<float4, 1, cudaReadModeElementType> mesh;

texture<float4, 1, cudaReadModeElementType> upper_bound_arr;
texture<float4, 1, cudaReadModeElementType> lower_bound_arr;

texture<uint, 1, cudaReadModeElementType> child_map_arr;
texture<uint, 1, cudaReadModeElementType> child_len_arr;

__device__ float3 make_float3(const float4 &a)
{
	return make_float3(a.x, a.y, a.z);
}

__device__ bool intersect_triangle(const float3 &origin, const float3 &direction, const float3 &v0, const float3 &v1, const float3 &v2, float3 &intersection)
{
	Matrix m = make_matrix(v1-v0, v2-v0, -direction);
	
	float determinant = det(m);

	if (determinant == 0.0f)
		return false;

	float3 b = origin-v0;

	float u1 = ((m.a11*m.a22 - m.a12*m.a21)*b.x + (m.a02*m.a21 - m.a01*m.a22)*b.y + (m.a01*m.a12 - m.a02*m.a11)*b.z)/determinant;

	if (u1 < 0.0f)
		return false;

	float u2 = ((m.a12*m.a20 - m.a10*m.a22)*b.x + (m.a00*m.a22 - m.a02*m.a20)*b.y + (m.a02*m.a10 - m.a00*m.a12)*b.z)/determinant;

	if (u2 < 0.0f)
		return false;

	float u3 = ((m.a10*m.a21 - m.a11*m.a20)*b.x + (m.a01*m.a20 - m.a00*m.a21)*b.y + (m.a00*m.a11 - m.a01*m.a10)*b.z)/determinant;

	if (u3 < 0.0f || (1-u1-u2) < 0.0f)
		return false;

	intersection = origin + direction*u3;

	return true;
}

__device__ int get_color(const float3 &direction, const float3 &v0, const float3& v1, const float3 &v2)
{
	int rgb = floorf(255*dot(normalize(cross(v1-v0,v2-v0)),-direction));

	return rgb << 16 | rgb << 8 | rgb;
}


__device__ bool intersect_box(const float3 &origin, const float3 &direction, const float3 &lower_bound, const float3 &upper_bound)
{
	float kmin, kmax, kymin, kymax, kzmin, kzmax;

	if (direction.x >= 0.0f)
	{
		kmin = (lower_bound.x - origin.x)/direction.x;
		kmax = (upper_bound.x - origin.x)/direction.x;
	}
	else
	{
		kmin = (upper_bound.x - origin.x)/direction.x;
		kmax = (lower_bound.x - origin.x)/direction.x;
	}

	if (kmax < kmin)
		return false;

	if (direction.y >= 0.0f)
	{
		kymin = (lower_bound.y - origin.y)/direction.y;
		kymax = (upper_bound.y - origin.y)/direction.y;
	}
	else
	{
		kymin = (upper_bound.y - origin.y)/direction.y;
		kymax = (lower_bound.y - origin.y)/direction.y;
	}

	if (kymax < kymin)
		return false;

	if (kymin > kmin)
		kmin = kymin;

	if (kymax < kmax)
		kmax = kymax;

	if (kmin > kmax)
		return false;

	if (direction.z >= 0.0f)
	{
		kzmin = (lower_bound.z - origin.z)/direction.z;
		kzmax = (upper_bound.z - origin.z)/direction.z;
	}
	else
	{
		kzmin = (upper_bound.z - origin.z)/direction.z;
		kzmax = (lower_bound.z - origin.z)/direction.z;
	}

	if (kzmax < kzmin)
		return false;

	if (kzmin > kmin)
		kmin = kzmin;

	if (kzmax < kmax)
		kmax = kzmax;

	if (kmin > kmax)
		return false;

	if (kmax < 0.0f)
		return false;

	return true;
}

__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
{
	float3 lower_bound = make_float3(tex1Dfetch(lower_bound_arr, i));
	float3 upper_bound = make_float3(tex1Dfetch(upper_bound_arr, i));

	return intersect_box(origin, direction, lower_bound, upper_bound);
}

extern "C"
{

__global__ void translate(int max_idx, float3 *pt, float3 v)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	pt[idx] += v;
}

__global__ void rotate(int max_idx, float3 *pt, float phi, float3 axis)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	pt[idx] = rotate(pt[idx], phi, axis);
}

__global__ void intersect_mesh(int max_idx, float3 *origin_arr, float3 *direction_arr, int *pixel_arr, int first_leaf)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	float3 origin = origin_arr[idx];
	float3 direction = direction_arr[idx];

	int *pixel = pixel_arr+idx;

	bool hit = false;

	float distance, min_distance;
	float3 intersection, min_intersection;
	
	if (!intersect_node(origin, direction, 0))
	{
		*pixel = 0;
		return;
	}

	int stack[100];

	int *head = &stack[0];
	int *node = &stack[1];
	*node = 0;

	int i;

	do
	{
		int child_map = tex1Dfetch(child_map_arr, *node);
		int child_len = tex1Dfetch(child_len_arr, *node);

		if (*node < first_leaf)
		{
			for (i=0; i < child_len; i++)
				if (intersect_node(origin, direction, child_map+i))
					*node++ = child_map+i;

			if (node == head)
				break;
			
			node--;

		}
		else // node is a leaf
		{
			for (i=0; i < child_len; i++)
			{
				int mesh_idx = child_map + 3*i;

				float3 v0 = make_float3(tex1Dfetch(mesh, mesh_idx));
				float3 v1 = make_float3(tex1Dfetch(mesh, mesh_idx+1));
				float3 v2 = make_float3(tex1Dfetch(mesh, mesh_idx+2));
				
				if (intersect_triangle(origin, direction, v0, v1, v2, intersection))
				{
					if (!hit)
					{
						*pixel = get_color(direction, v0, v1, v2);
						
						min_distance = norm(intersection-origin);
						min_intersection = intersection;
						hit = true;
						continue;
					}
					
					distance = norm(intersection-origin);
					
					if (distance < min_distance)
					{
						*pixel = get_color(direction, v0, v1, v2);
						
						min_distance = distance;
						min_intersection = intersection;
					}
				}

			} // triangle loop

			node--;

		} // node is a leaf

	} // while loop
	while (node != head);

	if (!hit)
		*pixel = 0;

} // intersect mesh

} // extern "c"