

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

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

    float biasdist = clamp(position.z * position.z, 8, 16); // shrink sample size as distance increases

    float ssptbias = 1 / biasdist;
    float stepSize = 64 * noise;

    int maxSteps = RTAO_TRACE_SAMPLES;

    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;
    float len = max(abs(direction.x) / texelSize.x, abs(direction.y) / texelSize.y) / stepSize;

    vec3 maxLengths = (step(0.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.0), maxSteps);

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

    stepv *= noise;
    spos += stepv;

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

        float sp = ld(texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r);
        float currZ = ld(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, invld2(sp));
        }

        spos += stepv;
    }

    return vec3(1.1);
}

vec2 R2_samples(int n)
{
    vec2 alpha = vec2(0.75487765, 0.56984026);
    return fract(alpha * n);
}
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);
}

float rtao(vec3 normal, vec3 fragpos)
{

    int samples = RTAO_SAMPLES;
    float occlusion = 0.0;

    float dither = blueNoise(gl_FragCoord.xy);
    vec2 rand = fract(vec2(dither, R2_dither(gl_FragCoord.xy)) + float(frameCounter % 10000) * vec2(0.75487765, 0.56984026));
    float jitterFactor = 64.0 * (rand.x + rand.y) + 1.0;

    // AO sampling parameters
    float spinAngle = rand.x;
    float r0 = rand.y;
    const float nbFactor = 552.920306;
    float spinFactor = spinAngle * TWO_PI;

    for (int i = 0; i < samples; i++)
    {
        // Compute sample coordinates in a spiral pattern
        float alpha = (float(i) + r0) / 7.0;
        float angle = alpha * nbFactor + spinFactor;
        vec2 tap = vec2(cos(angle), sin(angle)) * alpha;

        // Project tap to a hemisphere direction
        float theta = tap.x * TWO_PI;
        float y = tap.y;
        float xy = sqrt(1.0 - y * y);
        vec3 L = normalize(normal + vec3(sin(theta) * xy, cos(theta) * xy, y));
        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(rtPos);

            if (all(greaterThanEqual(previousPosition.xy, vec2(0.0))) && rtPos.z < 1.0)
            {
            }
            else
            {
                occlusion += 1.0;
            }
        }
    }

    return pow(occlusion / float(samples), 1.25);
}
