#ifdef GL_ES
precision highp float;
precision highp int;
precision highp sampler2D;
precision highp sampler2DArray;
#endif

in vec4 fragTexCoord;

layout(location=0) out vec4 fragOut0;

layout (binding = 0) uniform sampler2D composite;
layout (binding = 1) uniform sampler2D emissive;
layout (binding = 2) uniform sampler2D depth;
layout (binding = 3) uniform sampler3D volume_tex;
#ifdef NOISE
layout (binding = 4) uniform sampler3D noise_volume_tex;
#endif

layout (std140) uniform genericData {
	mat4 p_inv;
	mat4 v_inv;
	vec3 camera;
	float zNear;
	vec3 globalLightDirection;
	float zFar;
	vec3 globalLightDiffuse;
	float stepsize;
	vec3 nebPos;
	float opacitydistance;
	vec3 nebSize;
	float alphalim;
	vec3 nebulaColor;
	float udfScale;
	float emissiveSpreadFactor;
	float emissiveIntensity;
	float emissiveFalloff;
	float henyeyGreensteinCoeff;
	vec3 noiseColor;
	int directionalLightSampleSteps;
	float directionalLightStepSize;
	float noiseColorScale1;
	float noiseColorScale2;
	float noiseIntensity;
	float aspect;
	float fov;
};

const float sqrt4pi_inv = inversesqrt(4.0 * 3.14159);
const float beer_powder_norm = 3.0/2.0 * sqrt(3.0);

//This HG implementation assumes that the scatter vectors between which is theta both point away from the point of scatter
float henyey_greenstein(float cosTheta) {
	float radicant = 1.0 + henyeyGreensteinCoeff * henyeyGreensteinCoeff + 2.0 * henyeyGreensteinCoeff * cosTheta;
	return sqrt4pi_inv * (1.0 - henyeyGreensteinCoeff * henyeyGreensteinCoeff) / pow(radicant, 3.0 / 2.0);
}

/*
Base idea of sampling volumetrics from multiple LOD 3D-Textures with cloud models / noise as well as the shading for the diffuse
part of the volumetrics is adapted from Andrew Schneider's 2015 SIGGRAPH talk "The real-time volumetric cloudscapes of Horizon: Zero Dawn".
*/

void main()
{
	vec4 eyeDirection = p_inv * vec4(fragTexCoord.xy * 2.0 - 1.0, -1, 1);
    eyeDirection.w = 0;
    vec3 rayDirection = normalize((v_inv * eyeDirection).xyz);

	vec4 color_in = texture(composite, fragTexCoord.xy);

	vec3 lCorner = nebPos - nebSize * 0.5;
	vec3 rCorner = nebPos + nebSize * 0.5;
	
	vec3 t1 = (lCorner - camera) / rayDirection;
	vec3 t2 = (rCorner - camera) / rayDirection;

	vec3 tMin = min(t1, t2);
	vec3 tMax = max(t1, t2);

	vec2 fragcoordAngle = (fragTexCoord.xy - 0.5) * fov;
	fragcoordAngle.x *= aspect;
	float depth = (2.0 * zNear * zFar / (zFar + zNear - (2 * texture(depth, fragTexCoord.xy).x - 1) * (zFar - zNear))) * sqrt(1.0 + tan(fragcoordAngle.x) * tan(fragcoordAngle.x) + tan(fragcoordAngle.y) * tan(fragcoordAngle.y));

	//The t at which the ray enters the nebula cube
	float maxtMin = max(0,max(tMin.x, max(tMin.y, tMin.z)));
	//The t at which the ray leaves the nebula cube
	float mintMax = min(depth, min(tMax.x, min(tMax.y, tMax.z)));

	//The cumulative one minus alpha value of the nebula. Corresponds to the multiplier of the current step when additively drawing to a fragment for which the steps closer to the camera are already drawn.
	float cumOMAlpha = 1;
	//The cumulative distance already travelled through the nebula.
	float cumnebdist = 0;
	//The color of the current fragment
	vec3 cumcolor = vec3(0, 0, 0);

	//While not perfectly accurate, it's sufficient to assume all gradients for the 3d texture lookups are identical
	vec3 initialPos = (camera + rayDirection * maxtMin) / nebSize + 0.5;
	vec3 gradX = dFdx(initialPos);
	vec3 gradY = dFdy(initialPos);

	vec3 sidestep = 1.0 / vec3(textureSize(volume_tex, 0));

	for(float stept = maxtMin; stept < mintMax;) {
		//Step setup
		vec3 position = camera + rayDirection * stept - nebPos;
		vec3 sampleposition = position / nebSize + 0.5;
		vec4 volume_sample = textureGrad(volume_tex, sampleposition, gradX, gradY);

#ifdef DO_EDGE_SMOOTHING
		//Average 3D texel with texels on corner, in an attempt to reduce jagged edges.
		float stepcolor_alpha = volume_sample.a;
		if(cumOMAlpha > 0.8) {
			//Draw additional samples only if we currently encounter an edge
			stepcolor_alpha = stepcolor_alpha / 2.0 + (
				textureGrad(volume_tex, sampleposition + vec3(sidestep.x, sidestep.y, sidestep.z), gradX, gradY).a +
				textureGrad(volume_tex, sampleposition + vec3(sidestep.x, sidestep.y, -sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(sidestep.x, -sidestep.y, sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(sidestep.x, -sidestep.y, -sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(-sidestep.x, sidestep.y, sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(-sidestep.x, sidestep.y, -sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(-sidestep.x, -sidestep.y, sidestep.z), gradX, gradY).a + 
				textureGrad(volume_tex, sampleposition + vec3(-sidestep.x, -sidestep.y, -sidestep.z), gradX, gradY).a) / 16.0;
		}
#else
		float stepcolor_alpha = volume_sample.a;
#endif

		float stepsize_current = min(max(stepsize, step(stepcolor_alpha, 0.01) * volume_sample.x * udfScale), mintMax - stept);

		float stepalpha = -(pow(alphalim, 1.0 / (opacitydistance / stepsize_current)) - 1.0f) * stepcolor_alpha;
		//All the following computations are just required if we have a stepcoloralpha that is non-zero.
		if(stepcolor_alpha > 0.01)
		{
			//Noise - Technically, we need to do edge noise here but I'm not very happy with how it turned out. The cloud needs to be "eroded" by the noise as per the SIGGRAPH talk, but it's unclear as to how that can be achieved.
			//So for now, don't do edge noise. noise_volume_tex.ba are still reserved for noise data for edge erosion.

			//Diffuse
#ifdef NOISE
			//Mix the actual nebula color with the noise color depending on the product of the noise channels multiplied with intensity, smooth-clamped from 0 to 1
			vec3 stepcolor_neb = mix(nebulaColor, noiseColor,
				smoothstep(0, 1, (textureGrad(noise_volume_tex, position / noiseColorScale1, gradX, gradY).r + textureGrad(noise_volume_tex, position / noiseColorScale2, gradX, gradY).g) / 2.0 * noiseIntensity));
#else
			vec3 stepcolor_neb = nebulaColor;
#endif
			vec3 stepcolor_diffuse = stepcolor_neb * henyey_greenstein(dot(rayDirection, globalLightDirection));
			float directionalLightStep = 4.0 / float(directionalLightSampleSteps);
			float directionalLightDepth = 0.1;
			//From the current position in the volumetric cloud, sample towards the sun to test whether or not the current part of the nebula is bright or not
			for(int dlstep = 1; dlstep <= directionalLightSampleSteps; dlstep++){
				//TODO Sample these in a cone towards the sun, not just a line. Also add the high-distance sample for concave nebulae
				vec3 dlsteppos = (position - globalLightDirection * (dlstep * directionalLightStepSize)) / nebSize + 0.5;
				float dlstepalpha = textureGrad(volume_tex, dlsteppos, gradX, gradY).a * step(0, dlsteppos.x) * step(dlsteppos.x, 1) * step(0, dlsteppos.y) * step(dlsteppos.y, 1) * step(0, dlsteppos.z) * step(dlsteppos.z, 1);
				directionalLightDepth += dlstepalpha * directionalLightStep;
			}
			stepcolor_diffuse *= beer_powder_norm * (1 - exp(-directionalLightDepth * 2.0)) * exp(-directionalLightDepth);

			//Emissive
			cumnebdist += stepcolor_alpha * stepsize_current;
			vec3 emissive_lod = textureLod(emissive, fragTexCoord.xy, clamp(cumnebdist * emissiveSpreadFactor, 0, float(textureQueryLevels(emissive) - 1))).rgb;
			vec3 stepcolor_emissive = clamp(emissive_lod.rgb * pow(alphalim, 1.0 / (opacitydistance / ((depth - stept) * emissiveFalloff + 0.01))) * emissiveIntensity, 0, 1);

			//Step finish
			vec3 stepcolor = clamp(stepcolor_diffuse + stepcolor_emissive, 0, 1);
			cumcolor += stepalpha * cumOMAlpha * stepcolor;
		}

		cumOMAlpha *= 1.0 - stepalpha;
		stept += stepsize_current;

		if(cumOMAlpha < alphalim)
			break;
	}

	fragOut0 = vec4(cumOMAlpha * color_in.rgb + ((1.0 - cumOMAlpha) * cumcolor), 1);
}
