

vec3 rayTraceDiffuse(vec3 dir, vec3 position, float dither)
{
    int quality = int(8);
    float stepSizeAdjust = clamp(dither * 1.5, 0.25, 1.5);

    vec3 clipPosition = toClipSpace3(position);
    vec3 direction = normalize(toClipSpace3(position + dir) - clipPosition);
    vec3 maxLengths = ((step(0.0, direction) - clipPosition) / direction) * dither;
    float minLength = min(min(maxLengths.x, maxLengths.y), maxLengths.z);
    vec3 stepv = (direction * (minLength)*stepSizeAdjust) / quality;

    float tolerance = max(abs(stepv.z) * 10.0, 0.025 / (position.z * position.z));
    vec3 spos = clipPosition + stepv * (dither);

    float depthSample = texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r;

    if (depthSample < spos.z)
    {
        if (abs(spos.z - depthSample) < tolerance)
        {
            return vec3(spos.xy, depthSample);
        }
    }
    for (int i = 0; i <= quality; ++i)
    {
        vec3 spos = clipPosition + stepv * (float(i) + dither);
        float depthSample = texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r;

        if (depthSample < spos.z)
        {
            if (abs(spos.z - depthSample) < tolerance)
            {
                return vec3(spos.xy, depthSample);
            }
        }
    }

    return vec3(1.1);
}

#define SSGI

#define SSPTMIX1 1.0       //[0.0 0.001 0.01 0.015 0.025 0.05 0.1 0.125 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.9 1.0]
#define SSPTBIAS 0.035     //[0.0 0.001 0.01 0.015 0.025 0.05 0.1 0.125 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.9 1.0]
#define SSPTambient 0.25   // [0.01 0.025 0.05 0.075 0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.20 0.22 0.23 0.25 0.27 0.29 0.32 0.34 0.37 0.40 0.43 0.46 0.50 0.54 0.58 0.63 0.68 0.74 0.79 0.86 0.93 1.00 1.08 1.17 1.26 1.36 1.47 1.59 1.71 1.85 2.00 2.16 2.33 2.51 2.72 2.93 3.17 3.42 3.69 3.99 4.30 4.65 5.02 5.42 5.85 6.32 6.82 7.37 7.95 8.59 9.27 10.01 10.81 11.68 12.61 13.61 14.70 15.87 17.14 18.51 19.99 21.58 23.30 25.16 27.17 29.34 31.68 34.21 36.94 39.89 43.07 46.50 50.22 54.22 58.55 63.22 68.27 73.72 79.60 85.95 92.81 100.22 108.21 116.85 126.17 136.24 147.11 158.85 171.53 185.22 200.00]
#define SSPTweightmix 0.25 // [0.10 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.20 0.22 0.23 0.25 0.27 0.29 0.32 0.34 0.37 0.40 0.43 0.46 0.50 0.54 0.58 0.63 0.68 0.74 0.79 0.86 0.93 1.00 1.08 1.17 1.26 1.36 1.47 1.59 1.71 1.85 2.00 2.16 2.33 2.51 2.72 2.93 3.17 3.42 3.69 3.99 4.30 4.65 5.02 5.42 5.85 6.32 6.82 7.37 7.95 8.59 9.27 10.01 10.81 11.68 12.61 13.61 14.70 15.87 17.14 18.51 19.99 21.58 23.30 25.16 27.17 29.34 31.68 34.21 36.94 39.89 43.07 46.50 50.22 54.22 58.55 63.22 68.27 73.72 79.60 85.95 92.81 100.22 108.21 116.85 126.17 136.24 147.11 158.85 171.53 185.22 200.00]

#define STEPS 8         // [Unlimited 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99]
#define STEP_LENGTH 32. // [ 4.  5.  6.  7.  8.  9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 60. 90.]
#define SubSurfaceScattering

float linZ(float depth)
{
    return (2.0 * near) / (far + near - depth * (far - near));
    // l = (2*n)/(f+n-d(f-n))
    // f+n-d(f-n) = 2n/l
    // -d(f-n) = ((2n/l)-f-n)
    // d = -((2n/l)-f-n)/(f-n)
}

float invLinZ(float lindepth)
{
    return -((2.0 * near / lindepth) - far - near) / (far - near);
}

vec3 RT(vec3 dir, vec3 position, float noise)
{

    float ssptbias = SSPTBIAS;

    float stepSize = clamp(STEP_LENGTH * noise, 0.1, STEP_LENGTH);
    // float stepSize = 6;

    int maxSteps = int(STEPS);

    vec3 clipPosition = toClipSpace3(position);

    float rayLength = ((position.z + dir.z * sqrt(3.0) * far) > -sqrt(3.0) * near) ? (-sqrt(3.0) * near - position.z) / dir.z : sqrt(3.0) * far;

    vec3 end = toClipSpace3(position + dir * rayLength);
    vec3 direction = end - clipPosition; // convert to clip space

    float len = max(abs(direction.x) / texelSize.x, abs(direction.y) / texelSize.y) / stepSize;

    // get at which length the ray intersects with the edge of the screen
    vec3 maxLengths = (step(0., direction) - clipPosition) / direction;
    float mult = min(min(maxLengths.x, maxLengths.y), maxLengths.z);
    vec3 stepv = direction / len;
    int iterations = min(int(min(len, mult * len) - 2), maxSteps);

    // Do one iteration for closest texel (good contact shadows)
    vec3 spos = clipPosition + stepv / stepSize * 6;
    spos.xy += offsets[framemod8] * texelSize * 0.5;

    float sp = linZ(texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r);

    float currZ = linZ(spos.z);

    if (sp < currZ)
    {
        float dist = abs(sp - currZ) / currZ;

        if (dist <= 0.035)
            return vec3(spos.xy, invLinZ(sp));
    }

    spos += stepv;

    for (int i = 0; i < iterations; i++)
    {
        if (clamp(clipPosition.xy, 0, 1) != clipPosition.xy)
            break;

        float sp = linZ(texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r);

        float currZ = linZ(spos.z);
        if (sp < currZ)
        {
            if (spos.x < 0.0 || spos.y < 0.0 || spos.z < 0.0 || spos.x > 1.0 || spos.y > 1.0 || spos.z > 1.0)
                return vec3(1.1);
            float dist = abs(sp - currZ) / currZ;

            if (dist <= ssptbias)
                return vec3(spos.xy, invLinZ(sp));
        }

        spos += stepv;
    }

    return vec3(1.1);
}

// Adapted from "Sampling the GGX Distribution of Visible Normals",
// http://jcgt.org/published/0007/04/01/
vec3 SampleGGXVNDF(vec3 Ve, float alpha_x, float alpha_y, float U1, float U2)
{
    // Transforming the view direction to the hemisphere configuration
    vec3 Vh = normalize(vec3(alpha_x * Ve.x, alpha_y * Ve.y, Ve.z));

    // Orthonormal basis
    vec3 T1 = (Vh.z < 1.0f) ? normalize(cross(vec3(0, 0, 1), Vh)) : vec3(1, 0, 0);
    vec3 T2 = cross(Vh, T1);

    // Parameterization of the projected area
    float r = sqrt(U1);
    float phi = 2.0f * PI * U2;
    float t1 = r * cos(phi);
    float t2 = r * sin(phi);
    float s = 0.5f * (1.0f + Vh.z);
    t2 = (1.0f - s) * sqrt(1.0f - t1 * t1) + s * t2;

    // Reprojection onto hemisphere
    vec3 Nh = t1 * T1 + t2 * T2 + sqrt(max(0.0f, 1.0f - t1 * t1 - t2 * t2)) * Vh;

    // Transform the normal back to the ellipsoid configuration
    vec3 Ne = normalize(vec3(alpha_x * Nh.x, alpha_y * Nh.y, max(0.0f, Nh.z)));

    return Ne;
}

float vdc(uint bits)
{
    bits = (bits << 16u) | (bits >> 16u);
    bits = ((bits & 0x55555555u) << 1u) | ((bits & 0xAAAAAAAAu) >> 1u);
    bits = ((bits & 0x33333333u) << 2u) | ((bits & 0xCCCCCCCCu) >> 2u);
    bits = ((bits & 0x0F0F0F0Fu) << 4u) | ((bits & 0xF0F0F0F0u) >> 4u);
    bits = ((bits & 0x00FF00FFu) << 8u) | ((bits & 0xFF00FF00u) >> 8u);
    return float(bits) * 2.3283064365386963e-10;
}

// Hammersley point set
vec2 hs(uint i, uint N)
{
    return vec2(float(i) / float(N), vdc(i));
}

vec2 StratifiedSample(int i, int N)
{
    int sqrtN = int(sqrt(float(N)));
    int x = i % sqrtN;
    int y = i / sqrtN;
    return (vec2(x, y) + vec2(fract(sin(float(i)) * 43758.5453123), fract(cos(float(i)) * 43758.5453123))) /
           float(sqrtN);
}
vec3 importanceSampleGGX(vec2 Xi, vec3 N, float roughness, float rot)
{
    float a = roughness * roughness;

    float phi = 2.0 * PI * fract(Xi.x + rot);
    float cosTheta = sqrt((1.0 - Xi.y) / (1.0 + (a * a - 1.0) * Xi.y));
    float sinTheta = sqrt(1.0 - cosTheta * cosTheta);

    // From spherical coordinates to cartesian coordinates
    vec3 H;
    H.x = cos(phi) * sinTheta;
    H.y = sin(phi) * sinTheta;
    H.z = cosTheta;

    // From tangent-space vector to world-space sample vector
    vec3 up = abs(N.z) < 0.5 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
    vec3 tangent = normalize(cross(up, N));
    vec3 bitangent = cross(N, tangent);

    vec3 sampleVec = tangent * H.x + bitangent * H.y + N * H.z;

    return normalize(sampleVec);
}

void frisvad(in vec3 n, out vec3 f, out vec3 r)
{
    if (n.z < -0.999999)
    {
        f = vec3(0., -1, 0);
        r = vec3(-1, 0, 0);
    }
    else
    {
        float a = 1. / (1. + n.z);
        float b = -n.x * n.y * a;
        f = vec3(1. - n.x * n.x * a, b, -n.x);
        r = vec3(b, 1. - n.y * n.y * a, -n.y);
    }
}
mat3 CoordBase(vec3 n)
{
    vec3 x, y;
    frisvad(n, x, y);
    return mat3(x, y, n);
}
vec2 R2_samples(int n)
{
    vec2 alpha = vec2(0.75487765, 0.56984026);
    return fract(alpha * n);
}
float square(float x)
{
    return x * x;
}

// GGX distribution function G1
float g(float NdotL, float roughness)
{
    float k = (roughness * roughness) / 2.0;
    return NdotL / (NdotL * (1.0 - k) + k);
}
// const bool colortex4MipmapEnabled = true;

vec3 sampleHemisphere(float u1, float u2)
{
    float r = sqrt(u1);
    float theta = 2.0 * PI * u2;
    float x = r * cos(theta);
    float y = r * sin(theta);
    float z = sqrt(1.0 - u1);
    return vec3(x, y, z);
}
float G2OverG1SmithGGX(float NdotV, float NdotL, float alpha2)
{
    float a = sqrt(alpha2 + (NdotV - NdotV * alpha2) * NdotV);
    float b = sqrt(alpha2 + (NdotL - NdotL * alpha2) * NdotL);
    return clamp((NdotV + a) * NdotL / (NdotL * a + NdotV * b), 0, 1);
}
vec3 GenerateUnitVector(vec2 xy)
{
    const float tau = radians(360.0);

    xy.x *= tau;
    xy.y = xy.y * 2.0 - 1.0;
    return vec3(vec2(cos(xy.x), sin(xy.x)) * sqrt(1.0 - xy.y * xy.y), xy.y);
}

vec3 GenerateCosineVector(vec3 vector, vec2 xy)
{
    // Apparently this is actually this simple.
    // http://www.amietia.com/lambertnotangent.html
    // (cosine lobe around vector = lambertian BRDF)
    return normalize(vector + GenerateUnitVector(xy));
}
vec3 GenerateCosineVectorSafe(vec3 vector, vec2 xy)
{
    // Apparently this is actually this simple.
    // http://www.amietia.com/lambertnotangent.html
    // (cosine lobe around vector = lambertian BRDF)
    // This one deals with ther rare case where cosineVector == (0,0,0)
    // Can just normalize it instead if you are sure that won't be a problem
    vec3 cosineVector = vector + GenerateUnitVector(xy);
    float lenSq = dot(cosineVector, cosineVector);
    return lenSq > 0.0 ? cosineVector * inversesqrt(lenSq) : vector;
}

vec3 GetFacetGGX(
    vec3 viewDirection,
    vec2 roughness, // along x and y axis of input space used for the view direction
    vec2 xy         // uniform random numbers from 0 to 1 - x can be limited to < 1 to clamp tail
)
{
    const float tau = radians(360.0);
    // GGX VNDF sampling
    // http://www.jcgt.org/published/0007/04/01/

    // transform view direction to hemisphere (section 3.2)
    viewDirection = normalize(vec3(roughness * viewDirection.xy, viewDirection.z));

    // orthonrmal basis (section 4.1)
    float clsq = dot(viewDirection.yx, viewDirection.yx);
    vec3 T1 = vec3(clsq > 0.0 ? vec2(-viewDirection.y, viewDirection.x) * inversesqrt(clsq) : vec2(1.0, 0.0), 0.0);
    vec3 T2 = vec3(-T1.y * viewDirection.z, viewDirection.z * T1.x, viewDirection.x * T1.y - T1.x * viewDirection.y);

    // parameterization of the projected area (section 4.2)
    float r = sqrt(xy.x);
    float phi = tau * xy.y;

    float t1 = r * cos(phi);
    float tmp = clamp(1.0 - t1 * t1, 0.0, 1.0);
    float t2 = mix(sqrt(tmp), r * sin(phi), 0.5 + 0.5 * viewDirection.z);

    // reprojection onto hemisphere (section 4.3)
    vec3 normalH = t1 * T1 + t2 * T2 + sqrt(clamp(tmp - t2 * t2, 0.0, 1.0)) * viewDirection;

    // transform normal back to ellipsoid (section 3.4)
    return normalize(vec3(roughness * normalH.xy, normalH.z));
}

vec3 diffuseFunc(float skyLight, vec3 normal, vec3 fragpos, vec3 light)
{
    vec3 tempDiffuse = vec3(0.0);
    mat3 basis = CoordBase(normal);

    float occlusion = 0.0;
    for (int i = 0; i < SSPT_SAMPLES; i++)
    {
        int seed = frameCounter * int(SSPT_SAMPLES) + i;
        vec2 ij = fract(R2_samples(seed) + blueNoise4(gl_FragCoord.xy).rg);
        vec3 L = GenerateCosineVector(normal, ij);

        float NdotL = max(dot(normal, L), 0.0);

        if (NdotL > 0.0)
        {
            vec3 rtPos = RT(mat3(gbufferModelView) * L, fragpos.xyz, noise_standard(gl_FragCoord.xy));

            vec3 previousPosition = reproject(mat3(gbufferModelViewInverse) * toScreenSpace(rtPos) + gbufferModelViewInverse[3].xyz, false);
            if (all(greaterThanEqual(previousPosition.xy, vec2(0.0))) && rtPos.z < 1.0)
            {
                tempDiffuse += texture(colortex4, previousPosition.xy).rgb * NdotL;
            }
            else
            {
                occlusion += 1.0;
                float bounceAmount = float(L.y > 0.0) + clamp(-L.y * 0.1 + 0.1, 0.0, 1.0);
                tempDiffuse.rgb += skyFromTex2(L, colortex6).rgb * bounceAmount * pow(skyLight, 4.0) * NdotL * PI;
            }
        }
    }

    tempDiffuse /= float(SSPT_SAMPLES);
    float occlusionFactor = (occlusion / float(SSPT_SAMPLES));
    tempDiffuse *= occlusionFactor;
    return tempDiffuse;
}
