#version 330 core
#extension GL_ARB_compute_shader : require
#extension GL_ARB_shader_image_load_store : require

#define LOCAL_SIZE 8

uniform sampler3D texVolume; //input data, only intensity is set
uniform writeonly image3D imgTarget; //output data, xyz = gradient, alpha = intensity

uniform float fMaxValue; //for rescaling of alpha values (0,1]

ivec3 transformRighthanded(in ivec3 pos)
{
	ivec3 texSize = textureSize(texVolume,0);

	int newY = (texSize.y-1) - pos.y;
	//return ivec3(pos.xzy);
	return ivec3(pos.xz,newY);
}

float getVal(in ivec3 pos)
{
	//the voxel's value is rescaled to fit in [0,1] range
	return texelFetch(texVolume,pos,0).a / fMaxValue;
}

//gradient calculation compute shader
layout(local_size_x = LOCAL_SIZE, local_size_y = LOCAL_SIZE, local_size_z = LOCAL_SIZE) in;
void main(void)
{
	//invocation id equals to 3d coords here
	ivec3 pos = ivec3(gl_GlobalInvocationID);
	
	//original value	
	float value = getVal(pos);

	//gradient calculation using central differences
	//note the right handed transformation here (exchanging YZ, flipping sign of new Z)
	vec3 diff = vec3(
		getVal(pos + ivec3(1,0,0)) - getVal(pos - ivec3(1,0,0)),
		getVal(pos + ivec3(0,0,1)) - getVal(pos - ivec3(0,0,1)),
		getVal(pos - ivec3(0,1,0)) - getVal(pos + ivec3(0,1,0)));

	//vec3 diff = vec3(
	//	getVal(pos + ivec3(1,0,0)) - getVal(pos - ivec3(1,0,0)),
	//	getVal(pos + ivec3(0,1,0)) - getVal(pos - ivec3(0,1,0)),
	//	getVal(pos + ivec3(0,0,1)) - getVal(pos - ivec3(0,0,1)));

	//we can normalize here, this saves a normalize in the vertex shader per fragment
	float len = length(diff);
	if(len==0)
	{	
		diff=vec3(0,0,0);
	}
	else
		diff = diff / len;

	vec4 finalVal = vec4(diff,value);
	//store in the righthanded position
	//imageStore(imgTarget,pos,finalVal);
	imageStore(imgTarget,transformRighthanded(pos),finalVal);
}

