const mat3 tmp_r709 = mat3(0.64 / 0.33, 0.3 / 0.6, 0.15 / 0.06, 1.0, 1.0, 1.0, 0.03 / 0.33, 0.1 / 0.6, 0.79 / 0.06);
const vec3 lc_r709 = vec3(0.3127 / 0.329, 1.0, 0.3583 / 0.329) * inverse(tmp_r709);

const mat3 R709ToXyz = mat3(lc_r709 * tmp_r709[0], lc_r709, lc_r709 * tmp_r709[2]);

const mat3 RgbToXyz = R709ToXyz;

vec4 atrous4(sampler2D tex1) {

    ivec2 fragCoord = ivec2(gl_FragCoord.xy);
    vec2 uv = texelSize * gl_FragCoord.xy;
    float depth = texelFetch(depthtex0, fragCoord, 0).r;
    vec4 diffuse;
	#ifdef FILTER_REVERSE_ORDER
    const int stepSize = 1 << (FILTER_ITERATION_FINAL - FILTER_ITERATION);
	#else
    const int stepSize = 1 << FILTER_ITERATION;
	#endif

			// Init for filter - start weight and result with center sample
    float diffuseWeightAccum = 1.0;
    #if FILTER_STAGEDEF == 5    
    diffuse = texelFetch(tex1, ivec2(gl_FragCoord.xy * SSPTSCALE), 0);
    //diffuse.a = computeVariance(tex1, ivec2(gl_FragCoord.xy * SSPTSCALE));
    diffuse.a = texelFetch(colortex10, ivec2(gl_FragCoord.xy * SSPTSCALE), 0).y;
    #else
    diffuse = texelFetch(tex1, fragCoord, 0);
    #endif
    diffuse.a *= diffuse.a;
    if (diffuse.a > 0.0000001) {

			// Misc prep
        vec3 posCenter = ViewSpaceFromScreenSpace(vec3(uv, depth), gbufferProjectionInverse);
        vec3 normalCenter = (decodeNormal(texelFetch(colortex1, fragCoord, 0).zw));

        vec3 normalCenterVS = mat3(gbufferModelView) * normalCenter;
        float rcpPhiPos = 1.0 / ((FILTER_SIGMAP * stepSize) * texelSize.x * gbufferProjectionInverse[0].x * -posCenter.z);

        float diffuseLumCenter = dot(diffuse.rgb, RgbToXyz[1]);
        float rcpDiffusePhiLum = max(1.0 / (FILTER_SIGMAL * sqrt(diffuse.a)), 0.0);

			// Perform spatial filter
			// 3x3 seems to give similar quality as 5x5

        for (int x = -1; x <= 1; ++x) {
            for (int y = -1; y <= 1; ++y) {
                if (x == 0 && y == 0) {
                    continue;
                }

                ivec2 offset = ivec2(x, y) * stepSize;
                ivec2 samplePos = fragCoord + offset;

                vec3 posSample = vec3(texelSize * (samplePos + 0.5), texelFetch(depthtex0, samplePos, 0).r);
                posSample = ViewSpaceFromScreenSpace(posSample, gbufferProjectionInverse);
                float wPos = dot(normalCenterVS, posSample - posCenter) * rcpPhiPos;

					#ifdef FILTER_USE_NORMAL_WEIGHT
                vec3 normalSample = (decodeNormal(texelFetch(colortex1, samplePos, 0).zw));

                float wNormal = pow(dot(normalCenter, normalSample) * 0.5 + 0.5, FILTER_SIGMAN);
					#endif

                float wShared = exp2(-x * x - y * y); // gaussian weight
                wShared *= exp(-wPos * wPos);
					#ifdef FILTER_USE_NORMAL_WEIGHT
                wShared *= wNormal;
					#endif

					// fetch diffuse
    #if FILTER_STAGEDEF == 5    
                vec4 diffuseSample = texelFetch(tex1, ivec2(samplePos * SSPTSCALE), 0);
                //diffuseSample.a = sqrt(max(computeVariance(tex1, ivec2(samplePos * SSPTSCALE)), 0));
                diffuseSample.a = texelFetch(colortex10, ivec2(samplePos * SSPTSCALE), 0).y;

    #else
                vec4 diffuseSample = texelFetch(tex1, samplePos, 0);

    #endif

                diffuseSample.a *= diffuseSample.a;
                float diffuseLumSample = dot(diffuseSample.rgb, RgbToXyz[1]);

					// two types of luminance weights:
					// symmetric & asymmetric
					// symmetric methods considers the variances of both the sampled & filtered pixels
					// asymmetric only considers the variance of either the filtered or sampled pixel (not both)
					// symmetric methods are less stable
					// asymmetric methods mess with the overall brightness
					#ifdef FILTER_SYMMETRIC
                float rcpDiffusePhiLumSample = max(1.0 / (FILTER_SIGMAL * sqrt(diffuseSample.a)), 0);
                float wDiffuse = wShared * max(exp(-abs(diffuseLumCenter - diffuseLumSample) * sqrt(rcpDiffusePhiLum * rcpDiffusePhiLumSample)), 0.0); // low stability
					#else
                float wDiffuse = wShared * max(exp(-abs(diffuseLumCenter - diffuseLumSample) * rcpDiffusePhiLum), 0.0); // asymmetric
					#endif

                diffuse += diffuseSample * vec4(vec3(wDiffuse), wDiffuse * wDiffuse);
                diffuseWeightAccum += wDiffuse;
            }
        }

        diffuse.a = sqrt(diffuse.a);
        diffuse /= diffuseWeightAccum;
    }
    return diffuse;
}
