vec3 KuwaharaFilter(sampler2D tex, vec2 uv, vec2 texel)
{
    vec3 mean[4];
    vec3 sqrMean[4];

    for(int i = 0; i < 4; i++)
    {
        mean[i] = vec3(0.0);
        sqrMean[i] = vec3(0.0);
    }

    float n = float((KUWAHARA_RADIUS + 1) * (KUWAHARA_RADIUS + 1));

    // Cuadrante superior izquierdo
    for(int j = -KUWAHARA_RADIUS; j <= 0; j++)
    for(int i = -KUWAHARA_RADIUS; i <= 0; i++)
    {
        vec3 c = texture2D(tex, uv + vec2(i, j) * texel).rgb;
        mean[0] += c;
        sqrMean[0] += c * c;
    }

    // Cuadrante superior derecho
    for(int j = -KUWAHARA_RADIUS; j <= 0; j++)
    for(int i = 0; i <= KUWAHARA_RADIUS; i++)
    {
        vec3 c = texture2D(tex, uv + vec2(i, j) * texel).rgb;
        mean[1] += c;
        sqrMean[1] += c * c;
    }

    // Cuadrante inferior izquierdo
    for(int j = 0; j <= KUWAHARA_RADIUS; j++)
    for(int i = -KUWAHARA_RADIUS; i <= 0; i++)
    {
        vec3 c = texture2D(tex, uv + vec2(i, j) * texel).rgb;
        mean[2] += c;
        sqrMean[2] += c * c;
    }

    // Cuadrante inferior derecho
    for(int j = 0; j <= KUWAHARA_RADIUS; j++)
    for(int i = 0; i <= KUWAHARA_RADIUS; i++)
    {
        vec3 c = texture2D(tex, uv + vec2(i, j) * texel).rgb;
        mean[3] += c;
        sqrMean[3] += c * c;
    }

    float minVar = 1e20;
    vec3 result = vec3(0.0);

    for(int i = 0; i < 4; i++)
    {
        mean[i] /= n;
        sqrMean[i] /= n;

        vec3 variance = sqrMean[i] - mean[i] * mean[i];
        float sigma = variance.r + variance.g + variance.b;

        if(sigma < minVar)
        {
            minVar = sigma;
            result = mean[i];
        }
    }

    return result;
}
