

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);
}

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 rayTraceRTAO(vec3 dir, vec3 position, float dither)
{
    const int quality = RTAO_TRACE_SAMPLES;
    float stepSizeAdjust = clamp(dither, 0.25, 1.0);

    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)) / float(32 * dither);

    float tolerance = max(abs(stepv.z) * 10.0, 0.02 / (position.z * position.z));

    for (int i = 0; i <= quality; ++i)
    {
        vec3 spos = clipPosition + stepv * (float(i) + blueNoise4(gl_FragCoord.xy).rgr);
        float depthSample = texelFetch(depthtex0, ivec2(spos.xy / texelSize), 0).r;

        if (depthSample < spos.z)
        {
            if (abs(spos.z - depthSample) < tolerance)
            {
                // Intersection found, return 0.0 immediately
                return 0.0;
            }
        }
    }

    // No intersection found after all steps
    return 1.0;
}

float rtao(float skyLight, vec3 normal, vec3 fragpos)
{

    mat3 basis = CoordBase(normal);
    vec4 bnoise = blueNoise4(gl_FragCoord.xy);
    float occlusion = 0.0;

    for (int i = 0; i < RTAO_SAMPLES; i++)
    {
        int seed = frameCounter * int(RTAO_SAMPLES) + i;
        vec2 ij = fract(R2_samples(seed) + bnoise.rg);

        // Sample directions over the hemisphere oriented around the normal
        vec3 H = sampleHemisphere(ij.x, ij.y);
        vec3 L = basis * H;

        float NdotL = max(dot(normal, L), 0.0);
        if (NdotL > 0.1)
        {
            // Compute occlusion and use early exit if no occlusion is found
            float rayResult = rayTraceRTAO(mat3(gbufferModelView) * L, fragpos.xyz, ij.x);
            occlusion += rayResult;
        }
    }

    return occlusion / float(RTAO_SAMPLES);
}
