

float ComputeOcclusion(vec3 N, vec3 P)
{
    float fovScale = gbufferProjection[1][1] * 0.36397;
    float radiusPix = fovScale * 0.1;
    float maxDist2 = P.z * P.z * fovScale * 0.22592;
    float invMaxDist2 = 1.0 / maxDist2;
    float falloffFactor = -2.88539 * invMaxDist2;

    float noise = blueNoise(gl_FragCoord.xy);
    float baseAngle = noise * 6.4;
    float alpha = noise;

    vec2 viewResF = vec2(viewResolution);
    vec2 kernelScale = alpha * radiusPix * vec2(viewResolution.x, viewResolution.y * aspectRatio);
    vec2 invRes = 1.0 / viewResF;

    vec2 accum = vec2(0.0);

    // only compute sin/cos ONCE
    vec2 dir = vec2(cos(baseAngle), sin(baseAngle));
    ivec2 resMinusOne = ivec2(viewResolution) - ivec2(1);

    for (int i = 0; i < 4; ++i)
    {
        vec2 sampleCoord = gl_FragCoord.xy + dir * kernelScale;
        ivec2 tc = ivec2(sampleCoord);
        tc = clamp(tc, ivec2(0), resMinusOne);

        float dSample = texelFetch(depthtex2, tc, 0).r;
        vec3 posS = toScreenSpace(vec3((vec2(tc) + 0.5) * invRes, dSample));

        vec3 d = posS - P;
        float d2 = dot(d, d);

        float mask = step(1e-5, d2) * step(d2, maxDist2);
        float nd = dot(d, N) * inversesqrt(max(d2, 1e-5));
        float atten = exp2(d2 * falloffFactor);

        accum += vec2(mask * nd * atten, mask);

        // rotate dir by 90 degrees CCW
        dir = vec2(dir.y, -dir.x);
    }

    return smoothstep(0.5, 1.0, 1.0 - (accum.x / max(accum.y, 1e-5)));
}

float rtao(vec3 normal, vec3 fragpos)
{

    //   return ComputeOcclusion(worldToView(normal), fragpos);

    const int samples = 4;  // Increasing this may further smooth corners at the cost of performance
    const int samples2 = 7; // Increasing this may further smooth corners at the cost of performance
    float occlusion = 0.0;

    // Precompute transformation matrices and screen-space constants
    mat3 modelView3 = mat3(gbufferModelView);
    vec3 clipPos = toClipSpace3(fragpos);
    vec2 invTexelSize = 1.0 / texelSize;
    const float TWO_PI = 6.2831853;

    // Generate random factors for dithering and spiral pattern
    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 = 128.0 * dither + 4.0;

    // AO sampling parameters: base rotation and spiral offset
    float spinAngle = rand.x;
    float r0 = rand.y;
    float spinFactor = spinAngle * TWO_PI;

    for (int i = 0; i < samples; i++)
    {
        // Generate a spiral pattern sample coordinate
        float alpha = (float(i) + r0) / samples2;
        float angle = float(i) * 2.39996323 + 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 sampleDir = vec3(sin(theta) * xy, cos(theta) * xy, y);

        // Bias the sampling direction toward the surface normal
        vec3 L = normalize(normal + sampleDir);
        vec3 dir = modelView3 * L;

        // Compute intersection with the near plane.
        float tNear = (-near - fragpos.z) / dir.z;

        // If tNear is negative or too small, the sample might not be valid.
        tNear = max(tNear, 0.0);

        // Choose a ray length tRay (you can clamp it to the range you want)
        float tRay = clamp(tNear, 0.01, far); // or choose based on scene scale

        vec3 samplePos = fragpos + dir * tRay;
        vec3 end = toClipSpace3(samplePos);
        vec3 delta = end - clipPos;
        float len = max(abs(delta.x) * invTexelSize.x, abs(delta.y) * invTexelSize.y);
        vec3 spos = clipPos + (delta / len) * jitterFactor;

        // Fetch depth at the sample location and compute a difference from the center depth
        float sampleDepth = ld(texelFetch(depthtex2, ivec2(spos.xy * invTexelSize), 0).r);
        float centerDepth = ld(spos.z);
        float diff = abs(sampleDepth - centerDepth) / centerDepth;

        // Instead of a binary occlusion test, smoothly blend the occlusion contribution.
        // If the sample depth is in front (occluding), use smoothstep to create a soft threshold.

        float depthDerivative = length(fwidth(spos));
        float adaptiveThreshold = depthDerivative * 2.0; // The factor 0.5 can be tuned as needed.

        float occlusionContribution = (sampleDepth < centerDepth) ? smoothstep(0.0, adaptiveThreshold, diff) : 1.0;

        // Reduce occlusion for grazing angles by blending with the dot-product factor.
        // Compute a weight based on the cosine, then use a smoothstep function for non-linear blending.
        float blendFactor = clamp(dot(normal, L), 0.0, 1.0); // adjust the exponent for less/more aggressive blending
        occlusionContribution = occlusionContribution * blendFactor;

        occlusion += occlusionContribution;
    }

    float occl = occlusion;
    occl = occl * occl;
    return clamp(occl, 0.0, 1.0);
}
