#ifndef __MESH_H__
#define __MESH_H__

#include "intersect.h"
#include "geometry.h"

#define STACK_SIZE 500

/* Tests the intersection between a ray and a node in the bounding volume
   hierarchy. If the ray intersects the bounding volume and `min_distance`
   is less than zero or the distance from `origin` to the intersection is
   less than `min_distance`, return true, else return false. */
__device__ bool
intersect_node(const float3 &origin, const float3 &direction,
	       Geometry *geometry, const int &i, const float &min_distance)
{
    /* assigning these to local variables is faster for some reason */
    float3 lower_bound = geometry->lower_bounds[i];
    float3 upper_bound = geometry->upper_bounds[i];

    float distance_to_box;

    if (intersect_box(origin, direction, lower_bound, upper_bound,
		      distance_to_box)) {
	if (min_distance < 0.0f)
	    return true;

	if (distance_to_box > min_distance)
	    return false;

	return true;
    }
    else {
	return false;
    }

}

/* Finds the intersection between a ray and `geometry`. If the ray does
   intersect the mesh and the index of the intersected triangle is not equal
   to `last_hit_triangle`, set `min_distance` to the distance from `origin` to
   the intersection and return the index of the triangle which the ray
   intersected, else return -1. */
__device__ int
intersect_mesh(const float3 &origin, const float3& direction,
	       Geometry *geometry, float &min_distance,
	       int last_hit_triangle = -1)
{
    int triangle_index = -1;

    float distance;
    min_distance = -1.0f;

    if (!intersect_node(origin, direction, geometry, geometry->start_node,
			min_distance))
	return -1;

    unsigned int stack[STACK_SIZE];

    unsigned int *head = &stack[0];
    unsigned int *node = &stack[1];
    unsigned int *tail = &stack[STACK_SIZE-1];
    *node = geometry->start_node;

    unsigned int i;

    do
    {
	unsigned int first_child = geometry->node_map[*node];
	unsigned int stop = geometry->node_map_end[*node];

	while (*node >= geometry->first_node && stop == first_child+1) {
	    *node = first_child;
	    first_child = geometry->node_map[*node];
	    stop = geometry->node_map_end[*node];
	}
		
	if (*node >= geometry->first_node) {
	    for (i=first_child; i < stop; i++) {
		if (intersect_node(origin, direction, geometry, i,
				   min_distance)) {
		    *node = i;
		    node++;
		}
	    }
	    
	    node--;
	}
	else {
	    // node is a leaf
	    for (i=first_child; i < stop; i++) {
		if (last_hit_triangle == i)
		    continue;

		Triangle triangle = get_triangle(geometry, i);

		if (intersect_triangle(origin, direction, triangle,
				       distance)) {
		    if (triangle_index == -1) {
			triangle_index = i;
			min_distance = distance;
			continue;
		    }

		    if (distance < min_distance) {
			triangle_index = i;
			min_distance = distance;
		    }
		}
	    } // triangle loop

	    node--;

	} // node is a leaf

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

extern "C"
{

__global__ void
distance_to_mesh(int nthreads, float3 *positions, float3 *directions,
		 Geometry *g, float *distances)
{
    int id = blockIdx.x*blockDim.x + threadIdx.x;

    if (id >= nthreads)
	return;

    float3 position = positions[id];
    float3 direction = directions[id];
    direction /= norm(direction);

    float distance;

    int triangle_index = intersect_mesh(position, direction, g, distance);

    if (triangle_index == -1)
	distances[id] = 1e9;
    else
	distances[id] = distance;
}

__global__ void
color_solids(int first_triangle, int nthreads, int *solid_id_map,
	     bool *solid_hit, unsigned int *solid_colors, Geometry *geometry)
{
    int id = blockIdx.x*blockDim.x + threadIdx.x;

    if (id >= nthreads)
	return;

    int triangle_id = first_triangle + id;
    int solid_id = solid_id_map[triangle_id];
    if (solid_hit[solid_id])
	geometry->colors[triangle_id] = solid_colors[solid_id];
}

} // extern "C"

#endif