//-*-c-*-

__device__ Matrix inv(const Matrix&m, float& determinant)
{
	determinant = det(m);

	return make_matrix(m.a11*m.a22 - m.a12*m.a21,
			   m.a02*m.a21 - m.a01*m.a22,
			   m.a01*m.a12 - m.a02*m.a11,
			   m.a12*m.a20 - m.a10*m.a22,
			   m.a00*m.a22 - m.a02*m.a20,
			   m.a02*m.a10 - m.a00*m.a12,
			   m.a10*m.a21 - m.a11*m.a20,
			   m.a01*m.a20 - m.a00*m.a21,
			   m.a00*m.a11 - m.a01*m.a10)/determinant;
}

__device__ bool intersect_triangle(const float3 &x, const float3 &p, float3 *triangle, float3 &intersection)
{
	float3 v0 = triangle[0];
	float3 v1 = triangle[1];
	float3 v2 = triangle[2];

	float determinant;
	float3 u = inv(make_matrix(v1-v0,v2-v0,-p), determinant)*(x-v0);

	if (determinant == 0.0)
		return false;

	if (u.x < 0.0 || u.y < 0.0 || u.z < 0.0 || (1-u.x-u.y) < 0.0)
		return false;

	intersection = x + p*u.z;

	return true;
}

__device__ int get_color(const float3 &p, float3 *vertex)
{
	float3 v1 = vertex[1] - vertex[0];
	float3 v2 = vertex[2] - vertex[0];

	float3 normal = cross(v1,v2);

	float scale;
	scale = dot(normal,-p)/(norm(normal)*norm(p));
	if (scale < 0.0)
		scale = dot(-normal,-p)/(norm(normal)*norm(p));

	int rgb = floorf(255*scale);

	return rgb*65536 + rgb*256 + rgb;
}

extern "C"
{

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

	if (idx > max_idx)
		return;

	x[idx] += v;
}

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

	if (idx > max_idx)
		return;

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

__global__ void intersect_triangle_mesh(int max_idx, float3 *xarr, float3 *parr, int n, float3* mesh, int *pixelarr)
{
	int idx = blockIdx.x*blockDim.x + threadIdx.x;

	if (idx > max_idx)
		return;

	float3 x = xarr[idx];
	float3 p = parr[idx];
	int *pixel = pixelarr+idx;

	bool hit = false;

	float distance, min_distance;
	float3 intersection, min_intersection;

	int i;
	for (i=0; i < n; i++)
	{
		if (intersect_triangle(x, p, mesh+3*i, intersection))
		{
			if (!hit)
			{
				*pixel = get_color(p, mesh+3*i);

				min_distance = norm(intersection-x);
				min_intersection = intersection;
				hit = true;
				continue;
			}

			distance = norm(intersection-x);

			if (distance < min_distance)
			{
				*pixel = get_color(p, mesh+3*i);

				min_distance = distance;
				min_intersection = intersection;
			}
		}
	}
	
	if (!hit)
		*pixel = 0;
	
}

} // extern "c"