
#include <shaders/materials/commons_sampling.glsl>
#include <shaders/commons_hlsl.glsl>

float sampleShadowPCF(in sampler2DShadow smpl, in vec4 coords, out float in_frustum, int samples, float sampling_range)
{
	float fact = 0.0;
	in_frustum = 0.0;

	// no shading if fragment behind
	if (coords.w <= 0.0)
		return 1.0;

	// no shading if fragment outside of light frustum
	if (coords.x < -coords.w || coords.x > coords.w || coords.y < -coords.w || coords.y > coords.w || coords.z < 0.0)
		return 1.0;

	in_frustum = 1.0;

	// NOTE: we need to shift after projection. this works here but we dont do it in directional. most likely because projection matrix is simple so it works...
	//       but this shift in directional might actualy introduce tiny inaccuracies
	// NOTE: can we factor that into projection matrix? most likely we can!
	coords.xy = coords.xy * vec2(0.5) + vec2(0.5) * coords.w;

#ifdef SPIRV_VULKAN
	coords.y = coords.y /coords.w;
	coords.y = 1.0 - coords.y;
	coords.y = coords.y * coords.w;
#endif

	vec3 samp = coords.xyz / coords.w;

#if 1
	float sampling_r = 1.0;
	float y = -sampling_r;
	for(int iy = 0; iy < samples; iy++, y += (2.0 * sampling_r) / float(samples - 1))
	{
		float x = -sampling_r;
		for(int ix = 0; ix < samples; ix++, x += (2.0 * sampling_r) / float(samples - 1))
		{
			float shadow = textureProjLod(
				smpl,
				coords + vec4(x * sampling_range * coords.w, y * sampling_range * coords.w, 0.0, 0.0),
				0.0
			);
			fact += shadow;
		}
	}

	float f = fact * (1.0 / (samples * samples));
	return f;

#else

	int h_samples = samples * samples / 2;
	for(int n = 0; n < h_samples; n++)
	{
		vec2 xy = (HammersleySample2D(n, h_samples) - 0.5) * 2.0;
		float shadow = textureProj(smpl, coords + vec4(xy.x * sampling_range * coords.w, xy.y * sampling_range * coords.w, 0.0, 0.0));
		fact += shadow;
	}
	return fact * (1.0 / h_samples);
#endif
}

// Occluder search

float AvgBlockersDepthToPenumbra(float light_size, float z_shadowMapView, float avgBlockersDepth)
{
  float penumbra = light_size * (z_shadowMapView - avgBlockersDepth) / avgBlockersDepth;
  //penumbra *= penumbra;
  return penumbra;
}

float PenumbraFromOccluderSearch(in sampler2D smpl, float light_size, vec2 hash, vec4 coords, float sampling_range, int samplesCount, float shadowmap_bias)
{
	float avgBlockersDepth = 0.0f;
	float blockersCount = 0.0f;
	const float penumbraFilterMaxSize = sampling_range;

	if (coords.x < -coords.w || coords.x > coords.w || coords.y < -coords.w || coords.y > coords.w || coords.z < 0.0)
		return 0.0;

	vec2 shadowMapUV = coords.xy / coords.w * 0.5 + 0.5;
	shadowMapUV.y = 1.0 - shadowMapUV.y;

	float z_shadowMapView = coords.z / coords.w;

	{
		[[loop]]
		for(int i = 0; i < samplesCount; i++)
		{
			const float golden_ratio = 1.61803398875;
			//vec2 sampleUV = fract(hash.xy + float(i + ((globals.monotonic * samplesCount) & 127)) * golden_ratio) - 0.5;
			vec2 sampleUV = fract(hash.xy + float(i) * golden_ratio) - 0.5;
			//vec2 sampleUV = fract(hash.xy + float(i) * golden_ratio) - 0.5;
			//vec2 sampleUV = fract(hash.xy + (i + (globals.monotonic & 127)) * golden_ratio) - 0.5;
			//vec2 sampleUV = fract(HammersleySample2D(i, samplesCount) + hash.rg + float(globals.monotonic) * golden_ratio) - 0.5;
			//vec2 sampleUV = fract(HammersleySample2D(i, samplesCount) + hash.rg) - 0.5;

			sampleUV = shadowMapUV + sampleUV * penumbraFilterMaxSize;

			#if 0
			float sampleDepth = texture(smpl, sampleUV).r;
			if(sampleDepth < z_shadowMapView)
			{
				avgBlockersDepth += sampleDepth;
				blockersCount    += 1.0f;
			}
			#else
			vec4 sampleDepth = textureGather(smpl, sampleUV);
			if(sampleDepth.r < z_shadowMapView)
			{
				avgBlockersDepth += sampleDepth.r;
				blockersCount    += 1.0f;
			}
			if(sampleDepth.g < z_shadowMapView)
			{
				avgBlockersDepth += sampleDepth.g;
				blockersCount    += 1.0f;
			}
			if(sampleDepth.b < z_shadowMapView)
			{
				avgBlockersDepth += sampleDepth.b;
				blockersCount    += 1.0f;
			}
			if(sampleDepth.a < z_shadowMapView)
			{
				avgBlockersDepth += sampleDepth.a;
				blockersCount    += 1.0f;
			}
			#endif
		}
	}

	if(blockersCount > 0.0f)
	{
		avgBlockersDepth /= blockersCount;
		return AvgBlockersDepthToPenumbra(light_size, z_shadowMapView, avgBlockersDepth);
	}
	else
	{
		return 0.0f;
	}
}

float sampleShadowPCFNoiseOffset(in sampler2DShadow smpl, in vec4 coords, ivec2 noise_offset, out float in_frustum, int samples, float sampling_range, float shadowmap_bias)
{
	float fact = 0.0;
	in_frustum = 0.0;

	// no shading if fragment behind
	if (coords.w <= 0.0)
		return 1.0;

	// no shading if fragment outside of light frustum
	if (coords.x < -coords.w || coords.x > coords.w || coords.y < -coords.w || coords.y > coords.w || coords.z < 0.0)
		return 1.0;

	in_frustum = 1.0;

	// NOTE: we need to shift after projection. this works here but we dont do it in directional. most likely because projection matrix is simple so it works...
	//       but this shift in directional might actualy introduce tiny inaccuracies
	// NOTE: can we factor that into projection matrix? most likely we can!
	coords.xy = coords.xy * vec2(0.5) + vec2(0.5) * coords.w;

	//if (coords.w != 1.0)
	//	return 1.0;

#ifdef SPIRV_VULKAN
	coords.y = coords.y /coords.w;
	coords.y = 1.0 - coords.y;
	coords.y = coords.y * coords.w;
#endif

	//vec3 samp = coords.xyz / coords.w;
	float bias = shadowmap_bias; //-0.0001;
	float bias2 = 0.3;
	const float golden_ratio = 1.61803398875;

#if 0
	uint h_samples = 0;
	vec2 hash = texelFetch(s_BlueNoise, ivec3(noise_offset.xy & ivec2(127), 0), 0).rg;

	int samples_coarse = samples;//int(min(8, samples));	// with current scheme there is little point integrating > 8 elements
	int samples_fine   = samples - samples_coarse;

	// coarse
	float sampling_r = 1.0;
	for(int ii=0; ii < samples_coarse; ii ++)
	{
		//vec2 hash_for_sample = fract(hash.xy + float(h_samples + (globals.monotonic & 127)) * golden_ratio);
		vec2 hash_for_sample = fract(hash.xy + float(h_samples) * golden_ratio);
		hash_for_sample = hash_for_sample - 0.5;
		vec2 xy = vec2(hash_for_sample.x, hash_for_sample.y) * sampling_r;
		//vec2 xy = vec2(x, y);
		vec4 jittered_coords = coords + vec4(xy.x * sampling_range * coords.w, xy.y * sampling_range * coords.w, bias, 0.0);
		float shadow = textureProjLod(
			smpl,
			jittered_coords,
			0.0);
		fact += shadow;
		h_samples += 1;
	}

	float avg = fact * (1.0 / float(samples_coarse));
#if 0
	//return avg;
	if (abs(avg - 0.5) > 0.25 && abs(avg - 0.5) < 0.5) // add more in the in-out areas
	{
		for(int ii=0; ii < samples_fine; ii ++)
		{
			vec2 hash_for_sample = fract(hash.xy + float(h_samples) * golden_ratio);
			hash_for_sample = hash_for_sample - 0.5;
			vec2 xy = vec2(hash_for_sample.x, hash_for_sample.y) * sampling_r;
			//vec2 xy = vec2(x, y);
			vec4 jittered_coords = coords + vec4(xy.x * sampling_range * coords.w, xy.y * sampling_range * coords.w, bias, 0.0);
			float shadow = textureProjLod(
				smpl,
				jittered_coords,
				0.0);
			fact += shadow;
			h_samples += 1;
		}

		avg = fact * (1.0 / float(samples_coarse + samples_fine));
		return 100.0;
	}
#endif
	return avg;
#else

	int taken_samples = 0;
	int h_samples = samples;
	vec3 hash = texelFetch(s_BlueNoise, ivec3(noise_offset.xy & ivec2(127), globals.monotonic & 15), 0).rgb;

	[[loop]]
	for(int n = 0; n < h_samples; n++)
	{
		//vec2 xy = HammersleySample2D(n, h_samples);
		//vec2 xy = fract(HammersleySample2D(n, h_samples) + hash.rg + float(globals.monotonic & 127) * golden_ratio) - 0.5;
		//vec2 xy = fract(HammersleySample2D(n, h_samples) + hash.rg) - 0.5;
		vec2 xy = fract(hash.rg + float(n) * golden_ratio) - 0.5;

		vec4 jittered_coords = coords + vec4(xy.x * sampling_range * coords.w, xy.y * sampling_range * coords.w, bias, 0.0);
		float shadow = textureProjLod(
			smpl,
			jittered_coords,
			0.0);

		fact += shadow;
		taken_samples += 1;
	}

	return fact * (1.0 / float(h_samples));
#endif
}

float sampleShadow(in sampler2DShadow smpl, in vec4 coords, out float in_frustum, out vec4 projector_color)
{
	in_frustum = 0.0;

	// no shading if fragment behind
	if (coords.w <= 0.0)
		return 0.0;

	// no shading if fragment outside of light frustum
	if (coords.x < -coords.w || coords.x > coords.w || coords.y < -coords.w || coords.y > coords.w)
		return 1.0;

	in_frustum = 1.0;

	// NOTE: we need to shift after projection. this works here but we dont do it in directional. most likely because projection matrix is simple so it works...
	//       but this shift in directional might actualy introduce tiny inaccuracies
	// NOTE: can we factor that into projection matrix? most likely we can!
	coords.xy = coords.xy * vec2(0.5) + vec2(0.5) * coords.w;

#ifdef SPIRV_VULKAN
	coords.y = coords.y /coords.w;
	coords.y = 1.0 - coords.y;
	coords.y = coords.y * coords.w;
#endif

	vec3 samp = coords.xyz / coords.w;

	float shadow = textureProj(smpl, coords);
	// TODO: Reimplement
	projector_color = vec4(0.0);// texture(LightProjectorSamplers[light.projector_sampler], coords.xy / coords.w);// * 10.0;
	return shadow;
}

float light_calculate_spot_attenuation(in LightProperties light, in vec3 pos)
{
	float cutoff = light.cutoff;
	float light_distance = length(light.position.xyz - pos.xyz);
	float falloff = dot(light.direction.xyz, (pos.xyz - light.position.xyz) / light_distance);
	float attenuation = 0.0;

	if (falloff > cutoff && light_distance < light.range)
	{
		attenuation = 1.0 - (1.0 - falloff) / (1.0 - cutoff);
		attenuation *= (1.0 - clamp(light_distance / light.range, 0.0, 1.0));
		attenuation = pow(attenuation, 2.0);
	}
	
	return attenuation;
}

float light_calculate_point_attenuation(in LightProperties light, in vec3 pos)
{
	float light_distance = length(light.position.xyz - pos.xyz);
	float attenuation = 0.0;

	if (light_distance < light.range)
	{
		attenuation = 1.0;
		attenuation *= (1.0 - clamp(light_distance / light.range, 0.0, 1.0));
		attenuation = pow(attenuation, 2.0);
	}
	
	return attenuation;
}

float light_calculate_area_attenuation(in LightProperties light, in vec3 pos)
{
	float light_distance = length(light.position.xyz - pos.xyz);
	float attenuation = 0.0;

	if (light_distance < light.range)
	{
		attenuation = (1.0 - clamp(light_distance / light.range, 0.0, 1.0));
		attenuation = pow(attenuation, 2.0);
	}
	
	return attenuation;
}

const float c_MinReflectance = 0.04;

struct MaterialInfo
{
    float perceptualRoughness; // roughness value, as authored by the model creator (input to shader)
    vec3 reflectance0;         // full reflectance color (normal incidence angle)

    float alphaRoughness;      // roughness mapped to a more linear change in the roughness (proposed by [2])
    vec3 diffuseColor;         // color contribution from diffuse lighting

    vec3 reflectance90;        // reflectance color at grazing angle
    vec3 specularColor;        // color contribution from specular lighting
};

struct AngularInfo
{
    float NdotL;                  // cos angle between normal and light direction
    float NdotV;                  // cos angle between normal and view direction
    float NdotH;                  // cos angle between normal and half vector
    float LdotH;                  // cos angle between light direction and half vector
    float VdotH;                  // cos angle between view direction and half vector
};

AngularInfo getAngularInfo(vec3 normalized_point_to_light, vec3 normal, vec3 normalized_view)
{
    // Standard one-letter names
#if 1
    vec3 n = normal;                      // Outward direction of surface point
    vec3 v = normalized_view;             // Direction from surface point to view
    vec3 l = normalized_point_to_light;   // Direction from surface point to light
#else // sanity checking
	vec3 n = normalize(normal);
    vec3 v = normalize(normalized_view);
    vec3 l = normalize(normalized_point_to_light);
#endif

    vec3 h = normalize(l + v);            // Direction of the vector between l and v

    //float NdotL = clamp(dot(n, l), 0.0, 1.0);
    float NdotL = dot(n, l);
    //float NdotV = clamp(dot(n, v), 0.0, 1.0);
    float NdotV = dot(n, v);
    float NdotH = clamp(dot(n, h), 0.0, 1.0);
    float LdotH = clamp(dot(l, h), 0.0, 1.0);
    float VdotH = clamp(dot(v, h), 0.0, 1.0);

    AngularInfo angularInfo = {
        NdotL,
        NdotV,
        NdotH,
        LdotH,
        VdotH
	};

    return angularInfo;
}

float getPerceivedBrightness(vec3 vec)
{
    return sqrt(0.299 * vec.r * vec.r + 0.587 * vec.g * vec.g + 0.114 * vec.b * vec.b);
}

// https://github.com/KhronosGroup/glTF/blob/master/extensions/2.0/Khronos/KHR_materials_pbrSpecularGlossiness/examples/convert-between-workflows/js/three.pbrUtilities.js#L34
float solveMetallic(vec3 diffuse, vec3 specular, float oneMinusSpecularStrength)
{
    float specularBrightness = getPerceivedBrightness(specular);

    if (specularBrightness < c_MinReflectance)
	{
		return 0.0;
    }

    float diffuseBrightness = getPerceivedBrightness(diffuse);

    float a = c_MinReflectance;
    float b = diffuseBrightness * oneMinusSpecularStrength / (1.0 - c_MinReflectance) + specularBrightness - 2.0 * c_MinReflectance;
    float c = c_MinReflectance - specularBrightness;
    float D = b * b - 4.0 * a * c;

    return clamp((-b + sqrt(D)) / (2.0 * a), 0.0, 1.0);
}

// Lambert lighting
// see https://seblagarde.wordpress.com/2012/01/08/pi-or-not-to-pi-in-game-lighting-equation/
vec3 diffuse(MaterialInfo materialInfo)
{
    return materialInfo.diffuseColor / M_PI;
}

// The following equation models the Fresnel reflectance term of the spec equation (aka F())
// Implementation of fresnel from [4], Equation 15
vec3 specularReflection(MaterialInfo materialInfo, AngularInfo angularInfo)
{
    return materialInfo.reflectance0 + (materialInfo.reflectance90 - materialInfo.reflectance0) * pow(clamp(1.0 - angularInfo.VdotH, 0.0, 1.0), 5.0);
}

// Smith Joint GGX
// Note: Vis = G / (4 * NdotL * NdotV)
// see Eric Heitz. 2014. Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs. Journal of Computer Graphics Techniques, 3
// see Real-Time Rendering. Page 331 to 336.
// see https://google.github.io/filament/Filament.md.html#materialsystem/specularbrdf/geometricshadowing(specularg)
float visibilityOcclusion(MaterialInfo materialInfo, AngularInfo angularInfo)
{
    float NdotL = clamp(angularInfo.NdotL, 0.0, 1.0);
    float NdotV = clamp(angularInfo.NdotV, 0.0, 1.0);
    float alphaRoughnessSq = materialInfo.alphaRoughness * materialInfo.alphaRoughness;

    float GGXV = NdotL * sqrt(NdotV * NdotV * (1.0 - alphaRoughnessSq) + alphaRoughnessSq);
    float GGXL = NdotV * sqrt(NdotL * NdotL * (1.0 - alphaRoughnessSq) + alphaRoughnessSq);

    float GGX = GGXV + GGXL;
    if (GGX > 0.0)
    {
        return 0.5 / GGX;
    }
    return 0.0;
}

// The following equation(s) model the distribution of microfacet normals across the area being drawn (aka D())
// Implementation from "Average Irregularity Representation of a Roughened Surface for Ray Reflection" by T. S. Trowbridge, and K. P. Reitz
// Follows the distribution function recommended in the SIGGRAPH 2013 course notes from EPIC Games [1], Equation 3.
float microfacetDistribution(MaterialInfo materialInfo, AngularInfo angularInfo)
{
    float alphaRoughnessSq = materialInfo.alphaRoughness * materialInfo.alphaRoughness;
    float f = (angularInfo.NdotH * alphaRoughnessSq - angularInfo.NdotH) * angularInfo.NdotH + 1.0;
    return alphaRoughnessSq / (M_PI * f * f + 0.000001f);
}

vec3 fakeSkinBRDF(float NdotL)
{
	float curvature = 0.195;

	NdotL = mad(NdotL, 0.5, 0.5); // map to 0 to 1 range
	float curva = (1.0/mad(curvature, 0.5 - 0.0625, 0.0625) - 2.0) / (16.0 - 2.0); // curvature is within [0, 1] remap to normalized r from 2 to 16
	float oneMinusCurva = 1.0 - curva;
	float3 curve0;
	{
		float3 rangeMin = float3(0.0, 0.3, 0.3);
		float3 rangeMax = float3(1.0, 0.7, 0.7);
		float3 offset = float3(0.0, 0.06, 0.06);
		float3 t = saturate( mad(float3(NdotL), vec3(1.0) / (rangeMax - rangeMin), (offset + rangeMin) / (rangeMin - rangeMax)  ) );
		float3 lowerLine = (t * t) * float3(0.65, 0.5, 0.9);
		lowerLine.r += 0.045;
		lowerLine.b *= t.b;
		float3 m = float3(1.75, 2.0, 1.97);
		float3 upperLine = mad(float3(NdotL), m, float3(0.99, 0.99, 0.99) -m );
		upperLine = saturate(upperLine);
		float3 lerpMin = float3(0.0, 0.35, 0.35);
		float3 lerpMax = float3(1.0, 0.7 , 0.6 );
		float3 lerpT = saturate( mad(float3(NdotL), float3(1.0)/(lerpMax-lerpMin), lerpMin/ (lerpMin - lerpMax) ));
		curve0 = lerp(lowerLine, upperLine, lerpT * lerpT);
	}
	float3 curve1;
	{
		float3 m = float3(1.95, 2.0, 2.0);
		float3 upperLine = mad( float3(NdotL), m, float3(0.99, 0.99, 1.0) - m);
		curve1 = saturate(upperLine);
	}
	float oneMinusCurva2 = oneMinusCurva * oneMinusCurva;
	float3 brdf = lerp(curve0, curve1, mad(oneMinusCurva2, -1.0 * oneMinusCurva2, 1.0) );

	return brdf;
}

vec3 getPointShade(vec3 normalized_point_to_light, MaterialInfo materialInfo, vec3 normal, vec3 normalized_view)
{
    AngularInfo angularInfo = getAngularInfo(normalized_point_to_light, normal, normalized_view);

    if (angularInfo.NdotL > 0.0 || angularInfo.NdotV > 0.0)
    {
        // Calculate the shading terms for the microfacet specular shading model
        vec3 F    = specularReflection(materialInfo, angularInfo);
        float Vis = visibilityOcclusion(materialInfo, angularInfo);
        float D   = microfacetDistribution(materialInfo, angularInfo);

        // Calculation of analytical lighting contribution
        vec3 diffuseContrib = (vec3(1.0) - F) * diffuse(materialInfo);
        vec3 specContrib    = F * Vis * D;

        // Obtain final intensity as reflectance (BRDF) scaled by the energy of the light (cosine law)
        return clamp(angularInfo.NdotL, 0.0, 1.0) * (diffuseContrib + specContrib);
        //return fakeSkinBRDF(angularInfo.NdotL) * (diffuseContrib + specContrib);
    }

    return vec3(0.0, 0.0, 0.0);
}
