

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 rtao6(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 = samples + 3; // 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;
    }

    // Normalize and apply an exponent to adjust the falloff
    return clamp(pow(occlusion / float(samples), 1.25) * 1.25, 0.0, 1.0);
}

vec3 toClipSpace32(vec3 viewSpacePosition)
{
    // extract the diagonal of the projection matrix
    vec3 scale = vec3(gbufferProjection[0].x, gbufferProjection[1].y, gbufferProjection[2].z);

    // perform scale-and-add (MAD)
    vec3 projected = scale * viewSpacePosition + gbufferProjection[3].xyz;

    // perspective divide, remap from [-1,1] → [0,1]
    return projected / -viewSpacePosition.z * 0.5 + 0.5;
}

vec3 RT(vec3 dir, vec3 position, vec3 clipStart, float randomOffset)
{

    // Configuration
    float STEP_SIZE = 4;
    const int MAX_STEPS = 6;
    const float DEPTH_THRESH = 0.035; // relative tolerance for early hit

    // Compute maximum ray length in world‐space before clipping at near/far planes
    float zDir = dir.z * sqrt(3.0);
    float worldMaxLen = ((position.z + zDir * far) > -sqrt(3.0) * near) ? (-sqrt(3.0) * near - position.z) / dir.z : (sqrt(3.0) * far - position.z) / dir.z;

    // Clip‐space end point
    vec3 clipEnd = toClipSpace32(position + dir * worldMaxLen);
    vec3 clipDir = clipEnd - clipStart;
    if (dot(clipDir.xy, clipDir.xy) < 2.0)
        return vec3(1.1); // Ray not crossing screen, skip

    // Number of steps based on screen‐space distance
    float lenX = abs(clipDir.x) / texelSize.x;
    float lenY = abs(clipDir.y) / texelSize.y;
    float numSteps = max(lenX, lenY) / STEP_SIZE;
    int stepCount = min(int(numSteps), MAX_STEPS);

    // Limit ray so it stays inside the clip cube
    vec3 tMin = (step(vec3(0.0), clipDir) - clipStart) / clipDir;
    float t_clip = min(min(tMin.x, tMin.y), tMin.z);
    stepCount = min(stepCount, int(t_clip * numSteps) - 2);

    // Precompute the per‐step increment in clip‐space
    vec3 clipStep = clipDir / numSteps;

    // Initial sample position, nudged forward to avoid self‐intersection
    vec3 samplePos = clipStart + clipStep * (6.0 / STEP_SIZE);

    // First‐contact test
    float sceneDepth = ld(texelFetch(depthtex2, ivec2(samplePos.xy / texelSize), 0).r);
    float rayDepth = ld(samplePos.z);
    if (sceneDepth < rayDepth)
    {
        float relDist = abs(sceneDepth - rayDepth) / rayDepth;
        if (relDist <= DEPTH_THRESH)
            return vec3(samplePos.xy, sceneDepth);
    }

    // Randomized offset to reduce banding
    samplePos += clipStep * randomOffset;

    // March along the ray
    for (int i = 0; i < stepCount; ++i)
    {

        sceneDepth = ld(texelFetch(depthtex2, ivec2(samplePos.xy / texelSize), 0).r);
        rayDepth = ld(samplePos.z);

        if (sceneDepth < rayDepth)
        {
            float relDist = abs(sceneDepth - rayDepth) / rayDepth;
            if (relDist <= DEPTH_THRESH)
                return vec3(samplePos.xy, sceneDepth);
        }
        samplePos += clipStep;
    }

    // No hit: return off‐screen marker
    return vec3(1.1);
}
vec3 RT2(vec3 dir, vec3 position, vec3 clipStart, float randomOffset)
{
    vec3 rtPos = vec3(1.1);

    float quality = 4;

    vec3 clipPosition = toClipSpace3(position);
    vec3 direction = normalize(toClipSpace3(position + dir) - clipPosition);
    vec3 maxLengths = (step(0.0, direction) - clipPosition) / direction;
    float mult = min(min(maxLengths.x, maxLengths.y), maxLengths.z);

    float biasPower = 2; // Optional: adaptive front-loading

    vec3 dirScaled = direction * mult;
    float invQuality = 1.0 / quality;

    for (int i = 0; i <= int(quality); ++i)
    {
        float stepFrac = pow((float(i) + randomOffset) * invQuality, biasPower);
        vec3 spos = clipPosition + dirScaled * stepFrac;

        ivec2 sampleCoord = ivec2(spos.xy / texelSize);
        float depthSample = texelFetch2D(depthtex2, sampleCoord, 0).r;

        if (depthSample < spos.z)
        {
            float dist = abs(depthSample - spos.z) / spos.z;
            if (dist <= 0.005)
            {
                rtPos = vec3(spos.xy, depthSample);
                break;
            }
        }
    }

    return rtPos;
}

float rtao(vec3 N, vec3 fragpos)
{

    // Constants
    const float PHI = 1.618033988749895; // golden ratio
    const int NRAYS = 4;                 // number of hemisphere taps
    const float PI2 = 6.28318530718;     // 2π
    vec3 N_view = worldToView(N);

    // Choose an arbitrary “up” vector depending on N to build the tangent frame
    vec3 UpVector = abs(N_view.z) < 0.999 ? vec3(0.0, 0.0, 1.0) : vec3(1.0, 0.0, 0.0);
    vec3 T = normalize(cross(UpVector, N_view));
    vec3 B = cross(N_view, T);

    // Jitter using blue noise
    float dither = blueNoise(gl_FragCoord.xy);
    vec3 clipStart = toClipSpace32(fragpos);

    float occlusion = 0.0;
    for (int i = 0; i < NRAYS; i++)
    {
        // Vogel disk / golden ratio indexing
        float index = float(i) + dither;
        float u = index / float(NRAYS);
        float v = fract(index * PHI);

        // Convert to spherical coordinates
        float theta = acos(u);
        float phi = PI2 * v;

        // Sample direction in tangent space
        vec3 H;
        H.x = sin(theta) * cos(phi);
        H.y = sin(theta) * sin(phi);
        H.z = cos(theta);

        // Transform to world space: H.x * T + H.y * B + H.z * N
        vec3 worldDir = normalize(T * H.x + B * H.y + N_view * H.z);

        // Ray-trace from fragpos along worldDir
        vec3 hit = RT(worldDir, fragpos, clipStart, fract(index / PHI + dither));

        // If we hit something before reaching “infinite” (z < 1.0), count occlusion
        if (hit.z < 1.0)
            occlusion += 1.0;
    }

    // Return 1 − (occluded rays / total rays)
    return clamp(1.0 - (occlusion / float(NRAYS)), 0.0, 1.0);
}
