
float4 bbox_check(float4 origin, float4 dir)
{
	float left = 0.0f;
	float right = 1.0f;
	float down = 0.0f;
	float up = 1.0f;
	float back = 0.0f;
	float front = 1.0f;

	float tNearX, tFarX;
	float tNearY, tFarY;
	float tNearZ, tFarZ;

	// Slab X
	if (dir.x >= 0.0f)
	{
		tNearX = (left - origin.x) / dir.x;
		tFarX = (right - origin.x) / dir.x;
	}
	else if (dir.x < 0.0f)
	{
		tNearX = (right - origin.x) / dir.x;
		tFarX = (left - origin.x) / dir.x;
	}

	// Slab Y
	if (dir.y >= 0.0f)
	{
		tNearY = (down - origin.y) / dir.y;
		tFarY = (up - origin.y) / dir.y;
	}
	else if (dir.y < 0.0f)
	{
		tNearY = (up - origin.y) / dir.y;
		tFarY = (down - origin.y) / dir.y;
	}

	// Slab Z
	if (dir.z >= 0.0f)
	{
		tNearZ = (back - origin.z) / dir.z;
		tFarZ = (front - origin.z) / dir.z;
	}
	else if (dir.z < 0.0f)
	{
		tNearZ = (front - origin.z) / dir.z;
		tFarZ = (back - origin.z) / dir.z;
	}

	// We have to find now the the largest of the Near and the the smallest of the Far

	float tNearMax, tFarMin;

	if (tNearX > tNearY)
		tNearMax = tNearX;
	else
		tNearMax = tNearY;
	if (tNearZ > tNearMax)
		tNearMax = tNearZ;

	if (tFarX < tFarY)
		tFarMin = tFarX;
	else
		tFarMin = tFarY;
	if (tFarZ < tFarMin)
		tFarMin = tFarZ;

	float4 results;
	results.x = tNearMax - tFarMin;
	results.y = tFarMin;
	results.z = tNearMax;

	return results;
}

unsigned char getVoxel(__global unsigned char* vol, float4 pos, int vol_width)
{
	int xp = pos.x * vol_width;
	int yp = pos.y * vol_width;
	int zp = pos.z * vol_width;

	int ii = zp * vol_width * vol_width + yp * vol_width + xp;
	if (ii < 0) ii = 0;
	if (ii >= vol_width * vol_width * vol_width) ii = vol_width * vol_width * vol_width - 1;

	return vol[ii];
}

__kernel void volume0(__global unsigned int* screen, __global unsigned char* volume, __global float* cam_pos, __global float* cam_dir, __global float* cam_up, __global float* cam_right, int scr_width, int scr_height, int vol_width, int iso_depth, int ptime)
{
	int c_read = 0;

	int x = get_global_id(0);
	int y = get_global_id(1);
//	int t = ptime >> 4;

	float4 cm_pos = (float4)(cam_pos[0], cam_pos[1], cam_pos[2], cam_pos[3]);
	float4 cm_dir = (float4)(cam_dir[0], cam_dir[1], cam_dir[2], cam_dir[3]);
	float4 cm_up = (float4)(cam_up[0], cam_up[1], cam_up[2], cam_up[3]);
	float4 cm_right = (float4)(cam_right[0], cam_right[1], cam_right[2], cam_right[3]);

	int jump_size = 4;
	float small_step = 1.0f / vol_width;
	float step = small_step * (float)jump_size;

	float xt = ((float)x / (float)scr_width) - 0.5f;
	float yt = ((float)y / (float)scr_height) - 0.5f;

	float4 pscreen = xt * cm_right + yt * cm_up + cm_dir;

	float4 dir = pscreen / (pscreen.x * pscreen.x + pscreen.y * pscreen.y + pscreen.z * pscreen.z);
	float4 dirstep = dir * step;

	float4 pos = cm_pos + (float4)(0.5f, 0.5f, 0.5f, 0.0f);

	float4 res = bbox_check(pos, dirstep);
	if (res.x >= 0)
	{
		c_read = y >> 3;
		screen[y * scr_width + x] = ((c_read >> 2) << 16) | c_read;
		return;
	}

	if (res.z > 0.0f)
	{
		pos = pos + dirstep * res.z;
	}
	else res.z = 0.0f;

	float f_num_steps = res.y - res.z;


	int c_sum = 0;
	int a_sum = 0;

	float ns = 0;
	do
	{
		c_read = getVoxel(volume, pos, vol_width);
		int a_read = c_read >> 6;

		c_sum = (((255 - a_sum) * c_read) >> 8) + c_sum;
		a_sum = (((255 - a_sum) * a_read) >> 8) + a_sum;

		if (c_sum > 255) c_sum = 255;
		if (a_sum > 255) a_sum = 255;

		pos = pos + dirstep;
	}while(ns++ < f_num_steps);

	screen[y * scr_width + x] = (a_sum << 24) | (c_sum << 16) | (c_sum << 8) | c_sum;
}
