#version 400

#define rnd21(p)    fract(sin(dot(p.xy ,vec2(12.9898,78.233))) * 43758.5453)
#define rnd22(p)    fract(sin((p) * mat2(127.1,311.7,269.5,183.3) )*43758.5453)
#define SHOW_TILING ( false )
#define SHOW_SINGULARITIES ( false )
#define PI 3.14159265359

in vec2 uv;
out vec4 color;

// Mouse control
uniform vec4 _ScaleTranslate;

// Textures
uniform sampler2D Tinput; // Gaussian input T(I)
uniform sampler2D invT; // Inverse histogram transformation T^{-1}

// Decorrelated color space vectors and origin
uniform vec3 _colorSpaceVector1;
uniform vec3 _colorSpaceVector2;
uniform vec3 _colorSpaceVector3;
uniform vec3 _colorSpaceOrigin;
uniform vec3 _DXTScalers; // only with DXT compression (Section 1.6)

vec3 ReturnToOriginalColorSpace(vec3 color)
{
	vec3 result = 
		_colorSpaceOrigin +
		_colorSpaceVector1 * color.r +
		_colorSpaceVector2 * color.g +
		_colorSpaceVector3 * color.b;
	return result;
}

float B0cos(vec2 uv)
{
    float cosu = sin(uv.x*PI);
    float cosv = sin(uv.y*PI);
    return pow(cosu*cosv*cosu*cosv, 0.5);
}

float B1cos(vec2 uv)
{
    uv+=vec2(0.5);
    float cosu = sin(uv.x*PI);
    float cosv = sin(uv.y*PI);
    return pow(cosu*cosv*cosu*cosv, 0.5);
}

float BSingularity(vec2 uv)
{
    uv -= vec2(0.5);
    uv*=1.41421356237;
    const float a = PI/4.0;
    float cosA = cos(a);
    float sinA = sin(a);
    mat2 rotateA = mat2(cosA, sinA, -sinA, cosA);
    vec2 V = rotateA*uv;
    float cosu = sin(V.x*PI);
    float cosv = sin(V.y*PI);
    return 0.02f*cosu*cosv*cosu*cosv;
}

vec2 T0(vec2 uv)
{
    vec2 seed = floor(uv)*2;
    return seed;
}

vec2 T1(vec2 uv)
{
    vec2 seed;
    seed = floor(uv + vec2(0.5))*2 + 1;
    return seed;
}

vec2 rotateRandom(vec2 T, vec2 UV, float startAngle, float stopAngle)
{
    vec2 r = rnd22(T);
    float a = fract(r.x+r.y)*(stopAngle-startAngle) + startAngle;
    float cosA = cos(a);
    float sinA = sin(a);
    mat2 rotateA = mat2(cosA, sinA, -sinA, cosA);
    return rotateA*UV;
}

vec2 reflectRandom(vec2 T, vec2 UV)
{
    float r = (round(rnd21(T))-0.5)*2.0;
    return mat2(1, 0, 0, r)*UV;
}

vec3 PrimalAndDualBlending(vec2 uv)
{
	float startAngle = 0.0, stopAngle = 2.0*PI;
	vec3 O;
	vec3 b;
    float bSum;
	vec3 m = vec3(0.5); //vec3 m = texture(Tinput, uv, 99.).rgb;     // mean texture color
	vec2 duvdx = dFdx(uv);
	vec2 duvdy = dFdy(uv);
    if(SHOW_SINGULARITIES)
    {
        b = vec3(B0cos(uv), 
                 B1cos(uv), 
                 0);
        bSum = b.x+b.y;
    }
    else
    {
        b = vec3(B0cos(uv), 
                 B1cos(uv), 
                 BSingularity(uv));
        bSum = b.x+b.y+b.z;
    }
    b/=bSum;
    vec2 t0 = T0(uv);
    vec2 t1 = T1(uv);
    if(SHOW_TILING)
    {
        O = vec3(cos(t0.x*PI), cos(t0.y*PI), 1.0); 
    }
    else
    {
//Cancelled improvement consisting in adding symetry and rotations to translations (incompatible with the importance map).
//        O = 	(textureGrad(Tinput, reflectRandom(t0, rotateRandom(t0, uv-rnd22(t0), startAngle, stopAngle)), duvdx, duvdy).rgb - m)*b.x 
//			+ 	(textureGrad(Tinput, reflectRandom(t1, rotateRandom(t1, uv-rnd22(t1), startAngle, stopAngle)), duvdx, duvdy).rgb - m)*b.y;

        O = 	(textureGrad(Tinput, uv-rnd22(t0), duvdx, duvdy).rgb - m)*b.x 
			+ 	(textureGrad(Tinput, uv-rnd22(t1), duvdx, duvdy).rgb - m)*b.y;

        O = m + O/length(b);  // contrast preserving interp. cf https://www.shadertoy.com/view/4dcSDr
    }
	
	// Compute LOD level to fetch the prefiltered look-up table invT
	float LOD = textureQueryLod(Tinput, uv).y / float(textureSize(invT, 0).y);

	// Fetch prefiltered LUT (T^{-1})
	vec3 color;
	color.r	= texture(invT, vec2(O.r, LOD)).r;
	color.g	= texture(invT, vec2(O.g, LOD)).g;
	color.b	= texture(invT, vec2(O.b, LOD)).b;
	
	// Original color space
	color = ReturnToOriginalColorSpace(color);

	return color;
}

void main()
{
	color = vec4(PrimalAndDualBlending(uv * _ScaleTranslate.xy + _ScaleTranslate.zw), 1.0);
}
