//--// Settings //------------------------------------------------------------//

#include "/settings.glsl"
#include "/internalSettings.glsl"

//--// Uniforms //------------------------------------------------------------//

uniform sampler2D colortex0; // Accumulation buffer

uniform sampler3D colortex1; // 3D noise texture

uniform sampler3D colortex2; // Atmosphere Scattering LUT
uniform sampler2D colortex3; // Atmosphere Transmittance LUT

#define atlas   colortex4 // Atlas, base color
#define atlas_n colortex5 // Atlas, _n
#define atlas_s colortex6 // Atlas, _s
uniform sampler2D atlas;
uniform sampler2D atlas_n;
uniform sampler2D atlas_s;

uniform sampler2D noisetex;

uniform int frameCounter;

uniform float sunAngle;

uniform int worldDay;
uniform int worldTime;

//--// Camera uniforms

uniform float aspectRatio;

uniform float eyeAltitude;

uniform vec3 cameraPosition;
uniform vec3 previousCameraPosition;

uniform mat4 gbufferModelView;
uniform mat4 gbufferModelViewInverse;
uniform mat4 gbufferPreviousModelView;
uniform mat4 gbufferProjection;
uniform mat4 gbufferProjectionInverse;
uniform mat4 gbufferPreviousProjection;

//--// Shadow uniforms

#ifdef SHADOW_MAPPING
	uniform mat4 shadowModelView;
	uniform mat4 shadowProjection;
	uniform mat4 shadowProjectionInverse;
#endif

uniform sampler2D shadowtex0;
uniform sampler2D shadowcolor0;
uniform sampler2D shadowcolor1;

//--// Custom uniforms

uniform vec2 viewResolution;
uniform vec2 viewPixelSize;

uniform vec2 frameR2;

uniform vec3 sunVector;
uniform vec3 moonVector;
uniform vec3 shadowLightVector;

uniform float reset;

#if defined STAGE_VERTEX
	//--// Vertex Functions //------------------------------------------------//

	void main() {
		gl_Position = vec4(gl_Vertex.xy * 2.0 - 1.0, 1.0, 1.0);
	}
#elif defined STAGE_FRAGMENT
	//--// Fragment Outputs //------------------------------------------------//

	/* DRAWBUFFERS:0 */

	out vec4 fragdata;
	#define accumulated fragdata.rgb
	#define historyFrames fragdata.a
	//layout (location = 0, component = 0) out vec3  accumulated;
	//layout (location = 0, component = 3) out float historyFrames; // total accumulated frames so far for each pixel

	//--// Fragment Includes //-----------------------------------------------//

	#include "/include/utility.glsl"
	#include "/include/utility/color.glsl"
	#include "/include/utility/math.glsl"
	#include "/include/utility/rotation.glsl"
	#include "/include/utility/sampling.glsl"
	#include "/include/utility/voxelization.glsl"

	#ifdef SHADOW_MAPPING
	#include "/include/shadowDistortion.glsl"
	#endif

	#include "/include/atmosphere/constants.glsl"
	#include "/include/atmosphere/lookup.glsl"
	#include "/include/atmosphere/phase.glsl"
	#include "/include/atmosphere/transmittance.glsl"
	#include "/include/atmosphere/scattering.glsl"

	#include "/include/celestialConstants.glsl"

	uint lowbias32(uint x) {
		// https://nullprogram.com/blog/2018/07/31/
		x ^= x >> 16u;
		x *= 0x7feb352du;
		x ^= x >> 15u;
		x *= 0x846ca68bu;
		x ^= x >> 16u;
		return x;
	}
	uint randState;
	void InitRand(uint seed) { randState = lowbias32(seed); }
	uint RandNext() { return randState = lowbias32(randState); }
	#define RandNext2() uvec2(RandNext(), RandNext())
	#define RandNext3() uvec3(RandNext2(), RandNext())
	#define RandNext4() uvec4(RandNext3(), RandNext())
	#define RandNextF() (float(RandNext()) / float(0xffffffffu))
	#define RandNext2F() (vec2(RandNext2()) / float(0xffffffffu))
	#define RandNext3F() (vec3(RandNext3()) / float(0xffffffffu))
	#define RandNext4F() (vec4(RandNext4()) / float(0xffffffffu))

	#include "/include/clouds3D.fsh"

	#include "/include/voxelRaytracer.glsl"

	//--// Fragment Functions //----------------------------------------------//

	float FresnelDielectric(float VdotH, float eta) { // H = N when no roughness
		// Assumes non-polarized
		float p = 1.0 - (eta * eta * (1.0 - VdotH * VdotH));
		if (p <= 0.0) { return 1.0; } p = sqrt(p);

		vec2 r = vec2(VdotH, p);
		r = (eta * r - r.yx) / (eta * r + r.yx);
		return dot(r, r) * 0.5;
	}
	vec3 FresnelDielectric(float VdotH, vec3 eta) { // H = N when no roughness
		// Assumes non-polarized
		vec3 p = sqrt(Max0(1.0 - (eta * eta * (1.0 - VdotH * VdotH))));

		vec3 rs = (eta * VdotH - p) / (eta * VdotH + p);
		vec3 rp = (eta * p - VdotH) / (eta * p + VdotH);

		return Clamp01((rs * rs + rp * rp) * 0.5);
	}
	vec3 FresnelConductor(float VdotH, vec3 n2, vec3 k2) {
		// Found this while looking into complex fresnel online.
		// No idea about the original source.
		// Seems to assume that n1 and k1 are 1 and 0, respectively.
		vec3 ksq = k2 * k2;
		vec3 rs = (Pow2(n2 -        VdotH ) + ksq) / (Pow2(n2 +        VdotH ) + ksq);
		vec3 rp = (Pow2(n2 - (1.0 / VdotH)) + ksq) / (Pow2(n2 + (1.0 / VdotH)) + ksq);

		return Clamp01((rs + rp) * 0.5);
	}

	float DistributionGGX(float NdotH, float alpha2) {
		float p = (NdotH * alpha2 - NdotH) * NdotH + 1.0;
		return alpha2 / (pi * p * p);
	}

	float G1SmithGGX(float NdotV, float alpha2) {
		return Clamp01(2.0 * NdotV / (sqrt(alpha2 + (NdotV - NdotV * alpha2) * NdotV) + NdotV));
	}
	float G2SmithGGX(float NdotV, float NdotL, float alpha2) {
		float a = 2.0 * NdotV * NdotL;
		float b = NdotL * sqrt(alpha2 + (NdotV - NdotV * alpha2) * NdotV);
		float c = NdotV * sqrt(alpha2 + (NdotL - NdotL * alpha2) * NdotL);
		return Clamp01(a / (b + c));
	}
	vec2 G1G2SmithGGX(float NdotV, float NdotL, float alpha2) {
		vec2 delta = vec2(NdotV, NdotL);
		     delta = abs(sqrt((delta - delta * alpha2) * delta + alpha2) / delta);
		return Clamp01(2.0 / (delta.x + vec2(1.0, delta.y)));
	}
	float G2OverG1SmithGGX(float NdotV, float NdotL, float alpha2) {
		float a = sqrt(alpha2 + (NdotV - NdotV * alpha2) * NdotV);
		float b = sqrt(alpha2 + (NdotL - NdotL * alpha2) * NdotL);
		return Clamp01((NdotV + a) * NdotL / (NdotL * a + NdotV * b));
	}

	bool IsEmissiveId(int id) {
		return
		   id == 18 // lava
		|| id == 20 // torches on ground
		|| (id >= 147 && id <= 150) // torches on walls
		|| id == 240; // emissive full blocks
	}

	void ShadowLight(
		vec3 pos, ivec3 idx, vec3 viewDir,
		vec3 normal, vec3 mappedNormal,
		vec3 albedo, float roughness, vec3 eta, // TODO: Remaining material properties here
		inout vec3 diffuse, inout vec3 specular
	) {
		const int samples = 1;

		float NdotV = dot(mappedNormal, -viewDir);

		for (int i = 0; i < samples; ++i) {
			vec3 rayDir = GenerateConeVector(shadowLightVector, sunAngularRadius, RandNext2F());

			float NdotL = dot(mappedNormal, rayDir);
			if (NdotL <= 0.0 || dot(normal, rayDir) <= 0.0) { continue; }

			bool hit = RaytraceVoxelsAlphatest(pos, idx, rayDir);

			if (!hit) {
				#ifdef CLOUDS
				float cloudsTransmittance = EstimateTransmittance(SceneSpaceFromVoxelSpace(pos) + cameraPosition, rayDir);
				if (cloudsTransmittance <= 0.0) { continue; }
				#endif

				vec3 shadowLightIlluminance  = sunAngle < 0.5 ? sunIlluminance : moonIlluminance;
				     shadowLightIlluminance *= AtmosphereTransmittance(colortex3, atmosphere_planetRadius, rayDir.y);
				#ifdef CLOUDS
				     shadowLightIlluminance *= cloudsTransmittance;
				#endif

				float NdotH = dot( mappedNormal, normalize(rayDir - viewDir));
				float VdotH = dot(-viewDir,      normalize(rayDir - viewDir));

				float g1 = G1SmithGGX(NdotV, roughness * roughness);
				float g2 = G2SmithGGX(NdotV, NdotL, roughness * roughness);
				float d  = DistributionGGX(NdotH, roughness * roughness);
				vec3  f  = FresnelDielectric(VdotH, eta);

				mat3 rot = GetRotationMatrix(vec3(0.0, 0.0, 1.0), mappedNormal);
				vec3 facet = rot * GetFacetGGX(-viewDir * rot, vec2(roughness), RandNext2F());
				float MdotL = dot(facet, rayDir);

				vec3 sampleDiffuse  = shadowLightIlluminance * (g2 / g1) * Clamp01(MdotL) / pi;
				vec3 sampleSpecular = shadowLightIlluminance * g2 * d * f / (4.0 * NdotV);
				diffuse  += sampleDiffuse  / samples;
				specular += sampleSpecular / samples;
			}
		}
	}

	float mbnm_itp(vec3 rayDir, vec3 surfNormal, vec3 mapNormal) {
		vec3 tangent = -normalize(mapNormal - dot(mapNormal, surfNormal) * surfNormal);
		return Clamp01(dot(rayDir, tangent) / (dot(mapNormal, surfNormal) * dot(rayDir, surfNormal)));
	}
	void EnvironmentLight(vec3 viewPos, vec3 viewDir, out vec3 light) {
		const int samples = 1;
		const int bounces = MAX_BOUNCES + 1;

		light = vec3(0.0);
		for (int i = 0; i < samples; ++i) {
			vec3 contribution = vec3(1.0);

			vec3  hitEta;
			float hitRoughness;
			vec3  hitNormal;
			vec3  hitMappedNormal;
			mat3  hitTbnMatrix;
			vec3  hitAlbedo;

			vec3 rayPos = viewPos;
			ivec3 rayIdx = ivec3(rayPos);
			vec3 rayDir = viewDir;

			vec3 sampleResult = vec3(0.0);
			for (int bounce = 0; bounce < bounces; ++bounce) {
				if (bounce > 0) {
					#if defined WORLD_OVERWORLD
					#ifdef RANDOM_SUNLIGHT
					contribution *= 2.0;
					if (RandNextF() < 0.5)
					#endif
					{
						vec3 bounceDiffuse  = vec3(0.0);
						vec3 bounceSpecular = vec3(0.0);
						ShadowLight(
							rayPos, rayIdx, rayDir,
							hitNormal, hitMappedNormal,
							hitAlbedo, hitRoughness, hitEta, // TODO: Remaining material properties here
							bounceDiffuse, bounceSpecular
						);

						sampleResult += contribution * (bounceDiffuse * hitAlbedo + bounceSpecular);

						#ifdef RANDOM_SUNLIGHT
						break;
						#endif
					}
					#endif

					{ // calc new ray dir
						float NdotV = dot(hitMappedNormal, -rayDir);

						mat3 rot = GetRotationMatrix(vec3(0.0, 0.0, 1.0), hitMappedNormal);
						vec3 facetNormal = rot * GetFacetGGX(-rayDir * rot, vec2(hitRoughness), RandNext2F());

						vec3 fresnel = vec3(FresnelDielectric(dot(-rayDir, facetNormal), hitEta));

						float fresnelLum = dot(fresnel, RgbToXyz[1]);
						float albedoLum = dot(hitAlbedo, RgbToXyz[1]);
						float totalLum = fresnelLum + (1.0 - fresnelLum) * albedoLum;

						#ifdef PROBABILISTIC_PATH_TERMINATION
						float continuationProb = totalLum;//pow(totalLum, 0.5);
						if (RandNextF() < (1.0 - continuationProb)) { break; }
						contribution /= continuationProb;
						#endif

						float specularProb = fresnelLum / totalLum;
						bool isSpecularBounce = RandNextF() < specularProb;

						if (isSpecularBounce) {
							contribution *= fresnel / specularProb;
							rayDir = reflect(rayDir, facetNormal);
						} else {
							contribution *= (1.0 - fresnel) / (1.0 - specularProb);
							rayDir = GenerateCosineVector(facetNormal, RandNext2F());
							contribution *= hitAlbedo;
						}
					}

					float NdotL = dot(hitNormal, rayDir);
					if (NdotL <= 0.0) { break; }
				}

				vec4[2] hitVoxel; vec2 hitUv;
				bool hit = RaytraceVoxelsAlphatest(
					rayPos, rayIdx, rayDir,
					rayPos, rayIdx, hitNormal,
					hitVoxel, hitUv, hitTbnMatrix
				);

				if (hit) {
					vec4 base = textureLod(atlas, hitUv, 0.0);
					hitAlbedo = SrgbToLinear(base.rgb * hitVoxel[0].rgb);

					vec4 spec = textureLod(atlas_s, hitUv, 0.0);

					// Emissive light
					float emissive = float(IsEmissiveId(int((255.5 - OCTREE_LEVELS) - hitVoxel[0].a * 255.0)));
					emissive = spec.a < 254.5 / 255.0 ? max(emissive, spec.a * (255.0 / 254.0)) : emissive;
					sampleResult += contribution * hitAlbedo * EMISSIVE_MAX_LUMINANCE * emissive;

					vec4 norm = textureLod(atlas_n, hitUv, 0.0);

					hitMappedNormal = normalize(hitTbnMatrix * (norm.xyz * 2.0 - 1.0));

					hitRoughness = 1.0 - spec.r; hitRoughness *= hitRoughness;
					vec3 f0 = spec.g > (229.5 / 255.0) ? hitAlbedo : vec3(spec.g * spec.g);
					hitEta = (1.0 - sqrt(f0)) / (1.0 + sqrt(f0));
					if (spec.g > (229.5 / 255.0)) { hitAlbedo = vec3(0.0); }
				} else {
					// SKy
					#if defined WORLD_OVERWORLD
					vec3 sky = vec3(0.0);

					#ifdef CLOUDS
					float cloudsDistance = -1.0;
					vec4 clouds = Render3DClouds(rayDir, cloudsDistance);
					if (clouds.a < 1.0) {
						sky += sunIlluminance  * AtmosphereScattering(colortex2, vec3(0.0, rayPos.y + atmosphere_planetRadius, 0.0) + rayDir * cloudsDistance, rayDir, sunVector );
						sky += moonIlluminance * AtmosphereScattering(colortex2, vec3(0.0, rayPos.y + atmosphere_planetRadius, 0.0) + rayDir * cloudsDistance, rayDir, moonVector);
						sky *= clouds.a;
						sky += clouds.rgb;
						sky *= AtmosphereTransmittance(colortex3, rayPos.y + atmosphere_planetRadius, rayDir.y, cloudsDistance);
						sky += sunIlluminance  * AtmosphereScattering(colortex2, colortex3, vec3(0.0, rayPos.y + atmosphere_planetRadius, 0.0), rayDir, sunVector,  cloudsDistance);
						sky += moonIlluminance * AtmosphereScattering(colortex2, colortex3, vec3(0.0, rayPos.y + atmosphere_planetRadius, 0.0), rayDir, moonVector, cloudsDistance);
					} else
					#endif
					{
						sky  = sunIlluminance  * AtmosphereScattering(colortex2, atmosphere_planetRadius, rayDir.y, sunVector.y,  dot(rayDir, sunVector));
						sky += moonIlluminance * AtmosphereScattering(colortex2, atmosphere_planetRadius, rayDir.y, moonVector.y, dot(rayDir, moonVector));
					}

					sampleResult += contribution * sky;
					#endif

					break;
				}
			}

			light += Max0(sampleResult) / samples;
		}
	}

	void main() {
		InitRand(uint(gl_FragCoord.x) + uint(viewResolution.x) * (uint(gl_FragCoord.y) + uint(viewResolution.y) * uint(frameCounter)));

		vec2 uv = viewPixelSize * gl_FragCoord.xy;

		//--------------------------------------------------------------------//

		#ifdef DOF
			// calc dof offset here, so it can be used to tweak focal length for coma
			vec2 dofxy = RandNext2F();

			#if CAMERA_IRIS_BLADE_COUNT > 0
				const float ang = tau / CAMERA_IRIS_BLADE_COUNT;

				dofxy.x *= CAMERA_IRIS_BLADE_COUNT;
				float blade = floor(dofxy.x);
				dofxy.x = fract(dofxy.x);

				mat2 rot = GetRotationMatrix((blade + CAMERA_IRIS_BLADE_ROTATION) * ang);

				vec2 offsPoly = vec2(cos(ang / 2.0), sin(ang / 2.0) * (dofxy.x * 2.0 - 1.0));
				vec2 offsCirc = vec2(cos(ang * (dofxy.x - 0.5)), sin(ang * (dofxy.x - 0.5)));
				vec2 dofOffset = rot * mix(offsPoly, offsCirc, CAMERA_IRIS_BLADE_ROUNDING);
			#else // circular aperture
				vec2 dofOffset = vec2(cos(tau * dofxy.x), sin(tau * dofxy.x));
			#endif
			dofOffset *= sqrt(dofxy.y);
			float r = length(dofOffset);
		#endif

		const vec2 sensorSize = CAMERA_SCALE * 1e-3 * vec2(CAMERA_SENSOR_WIDTH_MM, CAMERA_SENSOR_HEIGHT_MM);
		#if CAMERA_FOCAL_LENGTH_MM > 0
		vec2 focalLength = CAMERA_SCALE * 1e-3 * vec2(CAMERA_FOCAL_LENGTH_MM * sensorSize.x / (aspectRatio * sensorSize.y), CAMERA_FOCAL_LENGTH_MM);
		#else
		vec2 focalLength = 0.5 * sensorSize * Diagonal(gbufferProjection).xy;
		#endif

		#ifdef DOF
		focalLength *= 1.0 + r * CAMERA_ABBERATION_COMA;
		#endif

		#if CAMERA_FOCAL_LENGTH_MM == 0
		vec2 fieldOfView = 2.0 * atan(sensorSize / (2.0 * focalLength));
		#endif

		//--//

		vec2 xy  = (uv + vec2(CAMERA_SHIFT_X, CAMERA_SHIFT_Y)) * 2.0 - 1.0;
		     xy += viewPixelSize * (RandNext2F() * 2.0 - 1.0);

		#if PROJECTION == PROJECTION_PANINI_COMPRESSED || PROJECTION == PROJECTION_PANINI
			#if CAMERA_FOCAL_LENGTH_MM > 0
			xy *= sensorSize / (2.0 * focalLength);
			#else
			xy.x *= aspectRatio;
			xy *= tan(fieldOfView.y / 2.0);
			#endif

			const float d = 1.0; // > 0
			float h = xy.x;
			float v = xy.y;

			float k = h*h / pow(d + 1.0, 2.0);

			float delta = k * k * d * d - (k + 1) * (k * d * d - 1);

			float cosPhi = (-k * d + sqrt(delta)) / (k + 1);
			float S = (d + 1.0) / (d + cosPhi);
			float phi = atan(h, S * cosPhi);
			float theta = atan(v, S);

			#if PROJECTION == PROJECTION_PANINI_COMPRESSED
			vec3 viewDir = vec3(sin(phi) * cos(theta), S * sin(theta) * cos(phi), -cos(phi) * cos(theta));
			#else
			vec3 viewDir = vec3(sin(phi) * cos(theta), sin(theta), -cos(phi) * cos(theta));
			#endif
		#elif PROJECTION == PROJECTION_EQUIRECTANGULAR
			xy  = xy * 0.5 + 0.5;
			xy *= vec2(tau, pi);

			vec3 viewDir = vec3(-sin(xy.x) * sin(xy.y), -cos(xy.y), cos(xy.x) * sin(xy.y));
		#elif PROJECTION == PROJECTION_EQUISOLID_ANGLE
			#if CAMERA_FOCAL_LENGTH_MM > 0
			xy *= sensorSize / (2.0 * focalLength);
			#else
			xy.x *= aspectRatio;
			float y =  sin(fieldOfView.y / 2.0);
			float z = -cos(fieldOfView.y / 2.0);
			xy *= sqrt(2.0 / (1.0 - z)) * y;
			#endif

			vec3 viewDir = vec3(xy * sqrt(1.0 - dot(xy, xy) / 4.0), dot(xy, xy) / 2.0 - 1.0);
		#elif PROJECTION == PROJECTION_EQUIDISTANT
			#if CAMERA_FOCAL_LENGTH_MM > 0
			xy *= sensorSize / (2.0 * focalLength);
			#else
			xy.x *= aspectRatio;
			xy *= fieldOfView.y / 2.0;
			#endif
			float theta = length(xy);
			xy = xy / theta;
			vec3 viewDir = vec3(xy * sin(theta), -cos(theta));
		#elif PROJECTION == PROJECTION_STEREOGRAPHIC
			#if CAMERA_FOCAL_LENGTH_MM > 0
			xy *= sensorSize / (4.0 * focalLength);
			#else
			xy.x *= aspectRatio;
			xy /= sin(fieldOfView.y / 2.0) / (1.0 - cos(fieldOfView.y / 2.0));
			#endif

			vec3 viewDir = vec3(2.0 * xy, dot(xy, xy) - 1.0) / (1.0 + dot(xy, xy));
		#else // PROJECTION == PROJECTION_RECTILINEAR
			xy *= sensorSize / (2.0 * focalLength);
			vec3 viewDir = normalize(vec3(xy, -1.0));
		#endif

		vec3 viewPos = vec3(0.0);

		#ifdef DOF
			vec2 apertureDiameter = focalLength / CAMERA_FSTOP;
			vec2 apertureRadius = apertureDiameter / 2.0;

			viewPos.xy += apertureRadius * dofOffset;

			// spherical abberation is (approximately) proportional to square of r
			float focus = CAMERA_FOCUS * (1.0 - r*r * CAMERA_ABBERATION_SPHERICAL);

			const float tiltTheta = radians(CAMERA_TILT_THETA);
			const float tiltPhi   = radians(CAMERA_TILT_PHI);
			const vec3 n = vec3(vec2(-sin(tiltPhi), cos(tiltPhi)) * sin(2.0 * tiltTheta), cos(2.0 * tiltTheta));
			vec3 p = vec3(0.0, 0.0, focus);

			float focusDist = dot(p, n) / abs(dot(viewDir, n));

			vec3 focusPoint = viewDir * focusDist;
			viewDir = normalize(focusPoint - viewPos);
		#endif

		//--//

		viewPos = VoxelSpaceFromSceneSpace(mat3(gbufferModelViewInverse) * viewPos + gbufferModelViewInverse[3].xyz);
		viewDir = mat3(gbufferModelViewInverse) * viewDir;

		//--// Raytrace lighting for current frame

		vec3 current = vec3(0.0);
		//if (viewDir.y > 0.0) {
			EnvironmentLight(viewPos, viewDir, current);
		//}

		current = Max0(current);

		//--// Accumulate & estimate variance

		bool    moved  = gbufferProjection != gbufferPreviousProjection;
		moved = moved || gbufferModelView  != gbufferPreviousModelView;
		moved = moved || cameraPosition    != previousCameraPosition;
		historyFrames = !moved ? texture(colortex0, uv).a + 1 : 1.0;

		vec3 history = texture(colortex0, uv).rgb;

		float wCurrent = 1.0 / max(historyFrames, 1);
		accumulated = mix(history, current, wCurrent);

		accumulated = clamp(accumulated, 0.0, 65e3);
	}
#endif
