#version 330 core

uniform sampler2D imageTexture;
uniform vec2 resolution;
uniform float filterRadius;
uniform uint frameId;

uniform float avgLuminance;
uniform float whitePoint;
uniform float keyValue;

// interpolated params
in vec2 v_texCoord;

layout(location=0) out vec3 outputColor;

// Reinhard's tone mapping.
// https://www.cs.utah.edu/~reinhard/cdrom/tonemap.pdf
// https://www.gamedev.net/topic/407348-reinhards-tone-mapping-operator/
vec3 ToneMapReinhard(in vec3 color)
{
	// RGB -> XYZ conversion
	//const mat3 RGB2XYZ = mat3(0.5141364, 0.3238786,  0.16036376,
	//						  0.265068,  0.67023428, 0.06409157,
	//						  0.0241188, 0.1228178,  0.84442666); 
	const mat3 RGB2XYZ = mat3(0.5141364,  0.265068,   0.0241188,
							  0.3238786,  0.67023428, 0.1228178,
							  0.16036376, 0.06409157, 0.84442666); 
	vec3 XYZ = RGB2XYZ * color;

	// XYZ -> Yxy conversion
	vec3 Yxy;
	Yxy.r = XYZ.g;                            // copy luminance Y
	Yxy.g = XYZ.r / (XYZ.r + XYZ.g + XYZ.b); // x = X / (X + Y + Z)
	Yxy.b = XYZ.g / (XYZ.r + XYZ.g + XYZ.b); // y = Y / (X + Y + Z)

	// (Lp) Map average luminance to the middlegrey zone by scaling pixel luminance
	float Lp = Yxy.r * keyValue / avgLuminance;                       
	// (Ld) Scale all luminance within a displayable range of 0 to 1
	Yxy.r = (Lp * (1.0f + Lp/(whitePoint * whitePoint)))/(1.0f + Lp);

	// Yxy -> XYZ conversion
	XYZ.r = Yxy.r * Yxy.g / Yxy. b;               // X = Y * x / y
	XYZ.g = Yxy.r;                                // copy luminance Y
	XYZ.b = Yxy.r * (1 - Yxy.g - Yxy.b) / Yxy.b;  // Z = Y * (1-x-y) / y

	// XYZ -> RGB conversion
	//const mat3 XYZ2RGB  = mat3(2.5651,-1.1665,-0.3986,
	//						  -1.0217, 1.9777, 0.0439, 
	//						   0.0753, -0.2543, 1.1892);
	const mat3 XYZ2RGB  = mat3(2.5651, -1.0217,  0.0753, 
							  -1.1665,  1.9777, -0.2543,
							  -0.3986,  0.0439,  1.1892);
	color = XYZ2RGB * XYZ;
	return color;
}

//vec3 ToneMap(in vec3 color)
//{
//	//	return color;			
//	float luminance = 0.299 * color.x + 0.587 * color.y + 0.114 * color.z;
//	// sigmoid mapping - seems to have not very intuitive parameters
//	// float m = contrast/100.0;
//	// float f = brightness/100.0;
//	// return (luminance/(luminance + pow(f*avgLuminance, m)))*color;
//	// simple linear clamp with contrast / brightness (currently data independent)
//	float m = contrast / 100.0;
//	float f = brightness / 100.0 - 1.0;
//	return color*m + vec3(1, 1, 1)*f;
//}

//Pre-computed table with exp(-sqrt(locX * locX + locY * locY))
//For 17x17 filter surroundings
#define USE_TABLE

#ifdef USE_TABLE
const float expSqrt[289]=float[289]
(
    0.243117f, 0.264803f, 0.286505f, 0.307510f, 0.326922f, 0.343695f, 0.356730f, 0.365028f, 0.367879f, 0.365028f, 0.356730f, 0.343695f, 0.326922f, 0.307510f, 0.286505f, 0.264803f, 0.243117f,
    0.264803f, 0.290127f, 0.315864f, 0.341199f, 0.365028f, 0.385979f, 0.402519f, 0.413175f, 0.416862f, 0.413175f, 0.402519f, 0.385979f, 0.365028f, 0.341199f, 0.315864f, 0.290127f, 0.264803f,
    0.286505f, 0.315864f, 0.346227f, 0.376709f, 0.406006f, 0.432348f, 0.453586f, 0.467505f, 0.472367f, 0.467505f, 0.453586f, 0.432348f, 0.406006f, 0.376709f, 0.346227f, 0.315864f, 0.286505f,
    0.307510f, 0.341199f, 0.376709f, 0.413175f, 0.449154f, 0.482454f, 0.510101f, 0.528677f, 0.535261f, 0.528677f, 0.510101f, 0.482454f, 0.449154f, 0.413175f, 0.376709f, 0.341199f, 0.307510f,
    0.326922f, 0.365028f, 0.406006f, 0.449154f, 0.493069f, 0.535261f, 0.571771f, 0.597269f, 0.606531f, 0.597269f, 0.571771f, 0.535261f, 0.493069f, 0.449154f, 0.406006f, 0.365028f, 0.326922f,
    0.343695f, 0.385979f, 0.432348f, 0.482454f, 0.535261f, 0.588411f, 0.637186f, 0.673488f, 0.687289f, 0.673488f, 0.637186f, 0.588411f, 0.535261f, 0.482454f, 0.432348f, 0.385979f, 0.343695f,
    0.356730f, 0.402519f, 0.453586f, 0.510101f, 0.571771f, 0.637186f, 0.702188f, 0.756155f, 0.778801f, 0.756155f, 0.702188f, 0.637186f, 0.571771f, 0.510101f, 0.453586f, 0.402519f, 0.356730f,
    0.365028f, 0.413175f, 0.467505f, 0.528677f, 0.597269f, 0.673488f, 0.756155f, 0.837967f, 0.882497f, 0.837967f, 0.756155f, 0.673488f, 0.597269f, 0.528677f, 0.467505f, 0.413175f, 0.365028f,
    0.367879f, 0.416862f, 0.472367f, 0.535261f, 0.606531f, 0.687289f, 0.778801f, 0.882497f, 1.000000f, 0.882497f, 0.778801f, 0.687289f, 0.606531f, 0.535261f, 0.472367f, 0.416862f, 0.367879f,
    0.365028f, 0.413175f, 0.467505f, 0.528677f, 0.597269f, 0.673488f, 0.756155f, 0.837967f, 0.882497f, 0.837967f, 0.756155f, 0.673488f, 0.597269f, 0.528677f, 0.467505f, 0.413175f, 0.365028f,
    0.356730f, 0.402519f, 0.453586f, 0.510101f, 0.571771f, 0.637186f, 0.702188f, 0.756155f, 0.778801f, 0.756155f, 0.702188f, 0.637186f, 0.571771f, 0.510101f, 0.453586f, 0.402519f, 0.356730f,
    0.343695f, 0.385979f, 0.432348f, 0.482454f, 0.535261f, 0.588411f, 0.637186f, 0.673488f, 0.687289f, 0.673488f, 0.637186f, 0.588411f, 0.535261f, 0.482454f, 0.432348f, 0.385979f, 0.343695f,
    0.326922f, 0.365028f, 0.406006f, 0.449154f, 0.493069f, 0.535261f, 0.571771f, 0.597269f, 0.606531f, 0.597269f, 0.571771f, 0.535261f, 0.493069f, 0.449154f, 0.406006f, 0.365028f, 0.326922f,
    0.307510f, 0.341199f, 0.376709f, 0.413175f, 0.449154f, 0.482454f, 0.510101f, 0.528677f, 0.535261f, 0.528677f, 0.510101f, 0.482454f, 0.449154f, 0.413175f, 0.376709f, 0.341199f, 0.307510f,
    0.286505f, 0.315864f, 0.346227f, 0.376709f, 0.406006f, 0.432348f, 0.453586f, 0.467505f, 0.472367f, 0.467505f, 0.453586f, 0.432348f, 0.406006f, 0.376709f, 0.346227f, 0.315864f, 0.286505f,
    0.264803f, 0.290127f, 0.315864f, 0.341199f, 0.365028f, 0.385979f, 0.402519f, 0.413175f, 0.416862f, 0.413175f, 0.402519f, 0.385979f, 0.365028f, 0.341199f, 0.315864f, 0.290127f, 0.264803f,
    0.243117f, 0.264803f, 0.286505f, 0.307510f, 0.326922f, 0.343695f, 0.356730f, 0.365028f, 0.367879f, 0.365028f, 0.356730f, 0.343695f, 0.326922f, 0.307510f, 0.286505f, 0.264803f, 0.243117f    
);

#define getExpSqrt(x,y) expSqrt[17*y + x]
#endif

void main()
{
	// show red square so that we know which shader is active
	/*if (gl_FragCoord.x < 20 && gl_FragCoord.y < 20) {
		outputColor.rgb = vec3(1, 0, 0);
		return;
	}*/

    float x = v_texCoord.x;
    float y = v_texCoord.y;

	float stepx = 1.0 / resolution.x;
	float stepy = 1.0 / resolution.y;

	// MK
	float samplesCount = 0;
	float sumWeight = 0.0;
	vec4 sumColor = vec4(0.0, 0.0, 0.0, 0.0);
	vec4 color = texture(imageTexture, vec2(x, y)).rgba;
	uint tmpNum = floatBitsToUint(color.a);
	uint countNum = tmpNum & uint(0xFFFF);
	uint frameIdNum = (tmpNum >> uint(16)) & uint(0x3FFF);

	// count samples found in first round twice
	bool overweight = true;

	if ((countNum > uint(0)) && (frameIdNum==frameId))
    {
#ifdef USE_TABLE
    float weight = getExpSqrt(8, 8);//expSqrt[8][8];
#else
    float weight = exp(-0.0*0.0 - 0.0*0.0);
#endif
		sumColor += (weight / float(countNum))*color;
		sumWeight += weight;
		samplesCount += float(countNum);
		if (overweight)
        {
			sumColor += (weight / float(countNum))*color;
			sumWeight += weight;
			overweight = false;
		}
	}

	int steps = 1;
#ifdef USE_TABLE
    const int locD = 1;
#else
	const float locD = 0.125;
#endif

	float startSamplesCount;

	while ((samplesCount < filterRadius) && (steps<9))
	{
		vec2 v = vec2(-stepx*steps, -stepy*steps);

#ifdef USE_TABLE
		int locX, locY;
        locX = locY = -steps;
#else
        float locX, locY;
        locX = locY = -(steps * locD);
#endif

		startSamplesCount = samplesCount;

		if (startSamplesCount > 0.0)
			overweight = false;

		for (int i = 0; i < 2 * steps; ++i)
		{
			color = texture(imageTexture, vec2(x + v.x, y + v.y)).rgba;
			tmpNum = floatBitsToUint(color.a);
			countNum = tmpNum & uint(0xFFFF);
			frameIdNum = (tmpNum >> uint(16)) & uint(0x3FFF);
			if ((countNum > uint(0)) && (frameIdNum==frameId)) {
#ifdef USE_TABLE
                float weight = getExpSqrt(8+locX, 8+locY);
#else
                float weight = exp(-sqrt(locX * locX + locY * locY));
#endif
				sumColor += (weight / float(countNum))*color;
				sumWeight += weight;
				samplesCount += float(countNum);

				if (overweight) {
					sumColor += (weight / float(countNum))*color;
					sumWeight += weight;
					if (startSamplesCount > 0.0)
						overweight = false;
				}
			}

			v.x += stepx;
			locX += locD;
		}

		for (int i = 0; i < 2 * steps; ++i)
		{
			color = texture(imageTexture, vec2(x + v.x, y + v.y)).rgba;
			tmpNum = floatBitsToUint(color.a);
			countNum = tmpNum & uint(0xFFFF);
			frameIdNum = (tmpNum >> uint(16)) & uint(0x3FFF);
			if ((countNum > uint(0)) && (frameIdNum==frameId)) {
#ifdef USE_TABLE
                float weight = getExpSqrt(8+locX, 8+locY);
#else
                float weight = exp(-sqrt(locX * locX + locY * locY));
#endif
				sumColor += (weight / float(countNum))*color;
				sumWeight += weight;
				samplesCount += float(countNum);

				if (overweight) {
					sumColor += (weight / float(countNum))*color;
					sumWeight += weight;
					if (startSamplesCount > 0.0)
						overweight = false;
				}
			}

			v.y += stepy;
			locY += locD;
		}

		for (int i = 0; i < 2 * steps; ++i)
		{
			color = texture(imageTexture, vec2(x + v.x, y + v.y)).rgba;
			tmpNum = floatBitsToUint(color.a);
			countNum = tmpNum & uint(0xFFFF);
			frameIdNum = (tmpNum >> uint(16)) & uint(0x3FFF);
			if ((countNum > uint(0)) && (frameIdNum==frameId)) {
#ifdef USE_TABLE
                float weight = getExpSqrt(8+locX, 8+locY);
#else
                float weight = exp(-sqrt(locX * locX + locY * locY));
#endif
				sumColor += (weight / float(countNum))*color;
				sumWeight += weight;
				samplesCount += float(countNum);

				if (overweight) {
					sumColor += (weight / float(countNum))*color;
					sumWeight += weight;
					if (startSamplesCount > 0.0)
						overweight = false;
				}
			}

			v.x -= stepx;
			locX -= locD;
		}

		for (int i = 0; i < 2 * steps; ++i)
		{
			color = texture(imageTexture, vec2(x + v.x, y + v.y)).rgba;
			tmpNum = floatBitsToUint(color.a);
			countNum = tmpNum & uint(0xFFFF);
			frameIdNum = (tmpNum >> uint(16)) & uint(0x3FFF);
			if ((countNum > uint(0)) && (frameIdNum==frameId)) {
#ifdef USE_TABLE
                float weight = getExpSqrt(8+locX, 8+locY);
#else
                float weight = exp(-sqrt(locX * locX + locY * locY));
#endif
				sumColor += (weight / float(countNum))*color;
				sumWeight += weight;
				samplesCount += float(countNum);

				if (overweight) {
					sumColor += (weight / float(countNum))*color;
					sumWeight += weight;
					if (startSamplesCount > 0.0)
						overweight = false;
				}
			}

			v.y -= stepy;
			locY -= locD;
		}
		steps += 1;
	}

	// = texture(imageTexture, vec2(x,y)).rgb;
	//outputColor.rgb = ToneMap(sumColor.rgb / sumWeight);
	outputColor.rgb = ToneMapReinhard(sumColor.rgb / sumWeight);
}
