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

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

#define STACK_SIZE 100

/* flattened triangle mesh */
texture<float4, 1, cudaReadModeElementType> mesh;

/* lower/upper bounds for the bounding box associated with each node/leaf */
texture<float4, 1, cudaReadModeElementType> upper_bounds;
texture<float4, 1, cudaReadModeElementType> lower_bounds;

/* map to child nodes/triangles and the number of child nodes/triangles */
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);
}

/* Test the intersection between a ray starting at `origin` traveling in the
   direction `direction` and the bounding box around node `i`. If the ray
   intersects the bounding box return true, else return false. */
__device__ bool intersect_node(const float3 &origin, const float3 &direction, const int &i)
{
	float3 lower_bound = make_float3(tex1Dfetch(lower_bounds, i));
	float3 upper_bound = make_float3(tex1Dfetch(upper_bounds, i));

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

/* Find the intersection between a ray starting at `origin` traveling in the
   direction `direction` and the global mesh texture. If the ray intersects
   the texture return the index of the triangle which the ray intersected,
   else return -1. */
__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int first_leaf)
{
	int triangle_idx = -1;

	float distance;
	float min_distance;

	if (!intersect_node(origin, direction, 0))
		return -1;

	int stack[STACK_SIZE];

	int *head = &stack[0];
	int *node = &stack[1];
	int *tail = &stack[STACK_SIZE-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 = 3*(child_map + 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, distance))
				{
					if (triangle_idx == -1)
					{
						triangle_idx = child_map + i;
						min_distance = distance;
						continue;
					}

					if (distance < min_distance)
					{
						triangle_idx = child_map + i;
						min_distance = distance;
					}
				}
			} // triangle loop

			node--;

		} // node is a leaf

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

	return triangle_idx;
}

extern "C"
{

/* Translate `points` by the vector `v` */
__global__ void translate(int max_idx, float3 *points, float3 v)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	*(points+idx) += v;
}

/* Rotate `points` through an angle `phi` counter-clockwise about the
   axis `axis` (when looking towards +infinity). */
__global__ void rotate(int max_idx, float3 *points, float phi, float3 axis)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	*(points+idx) = rotate(*(points+idx), phi, axis);
}

/* Trace the rays starting at `origins` traveling in the direction `directions`
   to their intersection with the global mesh. If the ray intersects the mesh
   set the pixel associated with the ray to a 32 bit color whose brightness is
   determined by the cosine of the angle between the ray and the normal of the
   triangle it intersected, else set the pixel to 0. */
__global__ void ray_trace(int max_idx, float3 *origins, float3 *directions, int first_leaf, int *pixels)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	float3 origin = *(origins+idx);
	float3 direction = *(directions+idx);
	direction /= norm(direction);

	int intersection_idx = intersect_mesh(origin, direction, first_leaf);

	if (intersection_idx == -1)
	{
		*(pixels+idx) = 0;
	}
	else
	{
		int mesh_idx = 3*intersection_idx;

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

		*(pixels+idx) = get_color(direction, v0, v1, v2);
	}
} // ray_trace

/* Propagate the photons starting at `origins` traveling in the direction
   `directions` to their intersection with the global mesh. If the ray
   intersects the mesh set the hit_solid array value associated with the
   photon to the triangle index of the triangle the photon intersected, else
   set the hit_solid array value to -1. */
__global__ void propagate(int max_idx, float3 *origins, float3 *directions, int first_leaf, int *hit_solids)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	float3 origin = *(origins+idx);
	float3 direction = *(directions+idx);
	direction /= norm(direction);

	*(hit_solids+idx) = intersect_mesh(origin, direction, first_leaf);
} // propagate

} // extern "c"