uniform sampler3D volume;
uniform sampler1D pal;

uniform bool alph_on;

uniform float width;
uniform float iso_depth;

uniform vec3 cam_pos;
uniform vec3 cam_dir;
uniform vec3 cam_up;
uniform vec3 cam_right;

const float jump_size = 4.0;

float small_step = 1.0 / width;
float step = small_step * jump_size;

float far_d = 0.0;
float near_d = 0.0;

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

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

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

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

	// Slab Z
	if (dir.z >= 0.0)
	{
		tNearZ = (back - origin.z) / dir.z;
		tFarZ = (front - origin.z) / dir.z;
	}
	else if (dir.z < 0.0)
	{
		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;

	far_d = tFarMin;
	near_d = tNearMax;

	// Final test
	return (tNearMax < tFarMin);
}


vec4 back_step_sample(vec3 pos, vec3 dir, int num, int nsteps)
{
	vec4 vbackhit;
	for (int i=0; i<num; i++)
	{
		pos = pos + dir;
		vbackhit = texture3D(volume, pos);
		vbackhit = texture1D(pal, vbackhit.r);

		if (vbackhit.r <= iso_depth)
			break;
	}

	vec4 vbc;
	if (alph_on)
	{
		float l = 0.1;
		float sum_a = 0.0;
		vbackhit = vec4(0.0);

		do
		{
			vbc = texture3D(volume, pos);
			vbc = texture1D(pal, vbc.r);

			vbackhit += l * vbc.a * vbc;
			sum_a = sum_a + l * vbc.a;
			pos = pos - dir;
		} while(sum_a < 1.0 && nsteps-- > 0);
	}

	return vbackhit;
}

void main()
{
	float xt = gl_TexCoord[0].x - 0.5;
	float yt = gl_TexCoord[0].y - 0.5;

	vec3 pscreen = xt * cam_right + yt * cam_up + cam_dir;

	vec3 dir = normalize(pscreen);
	vec3 dirstep = dir * step;
	vec3 small_backstep = -dir * small_step;

	vec3 pos = cam_pos + vec3(0.5, 0.5, 0.5);


	bool bchk = bbox_check(pos, dirstep);
	if (!bchk)
	{
	 	gl_FragColor = vec4(0.0, 0.0, 1.0, 1.0);
		return;
	}


	if (near_d > 0.0)
	{
		vec3 mdir = dirstep * near_d;
		pos = pos + mdir;
	}
	else near_d = 0;

	int f_num_steps = int(far_d - near_d);


	float c;
	vec4 vgfx = vec4(0.0);

	int ns = 0;
	do
	{
		vec4 vb = texture3D(volume, pos);

		if (vb.r > iso_depth)
		{
			vb = back_step_sample(pos, small_backstep, int(jump_size), f_num_steps - ns);
			c = vb.r;
			break;
		}

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

	vgfx = texture1D(pal, c);
	gl_FragColor = vgfx;

}
