
float unpackRoughness(float x)
{
    float r = 1.0 - x;
    return r * r;
}

vec3 rayTrace(vec3 dir, vec3 position, float dither, float rayQuality)
{
    float quality = 32;

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

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

    for (int i = 0; i <= int(rayQuality); ++i)
    {

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

    return vec3(1.1); // Indicates no intersection found
}

vec3 RT5(vec3 dir, vec3 position, float dither, float steps)
{

    float quality = steps;

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

    vec3 stepv = (direction * minLength) / quality;
    float tolerance = max(abs(stepv.z) * 2.0, 0.01 / (position.z * position.z));

    for (int i = 0; i <= int(quality); ++i)
    {

        vec3 spos = clipPosition + stepv * (i + dither);

        //  spos.xy += offsets[framemod8] * texSize * 0.5;

        float depthSample = texelFetch2D(depthtex0, ivec2(spos.xy / texelSize), 0).r;
        if (depthSample < spos.z && abs(spos.z - depthSample) < tolerance)
        {
            return vec3(spos.xy, depthSample);
        }
    }

    return vec3(1.1); // Indicates no intersection found
}
float invld2(float lindepth)
{
    return -((2.0 * near / lindepth) - far - near) / (far - near);
}

vec3 RT_SIMPLE(vec3 dir, vec3 position, float dither, float roughness)
{
    // Constants
    const float DEFAULT_NOISE = 0.5;
    const int MAX_STEPS = 2;
    const vec3 DEFAULT_RT_POS = vec3(1.1); // Indicates no intersection found
    const float ROUGHNESS_THRESHOLD = 0.25;

    // Initialize return position
    vec3 rtPos = DEFAULT_RT_POS;

    // Compute the farthest point in world space and convert to clip space
    vec3 farthestWorldPos = position + dir * DEFAULT_NOISE * 64.0;
    vec3 farthestClipPos = toClipSpace3(farthestWorldPos);

    // Check if the farthest point is within normalized device coordinates (NDC) bounds
    if (farthestClipPos.x > 0.0 && farthestClipPos.y > 0.0 &&
        farthestClipPos.x < 1.0 && farthestClipPos.y < 1.0)
    {
        // Sample the depth texture at the farthest point's screen space coordinates
        float sampledDepth = texture(depthtex0, farthestClipPos.xy).r;

        // If the sampled depth is closer than the farthest point, update rtPos
        if (sampledDepth < 1.0)
        {
            rtPos = farthestClipPos;
        }
    }

    // Early exit for materials with high roughness
    if (roughness > ROUGHNESS_THRESHOLD)
    {
        return rtPos;
    }

    // Convert the current position to clip space
    vec3 clipPosition = toClipSpace3(position);

    // Calculate the normalized direction in clip space
    vec3 clipDirection = normalize(toClipSpace3(position + dir) - clipPosition);

    // Determine the maximum lengths along each axis before exiting the clip space
    vec3 maxLengths = (step(0.0, clipDirection) - clipPosition) / clipDirection;
    float minLength = min(min(maxLengths.x, maxLengths.y), maxLengths.z);

    // Compute the step vector based on the minimum length and number of steps
    vec3 stepVector = clipDirection * (minLength / float(MAX_STEPS));

    // Tolerance based on the z-component of the step vector
    float tolerance = abs(stepVector.z);

    // Precompute the scaling factor for texel coordinates
    vec2 texCoordScale = 1.0 / texelSize;

    // Ray marching loop to find the intersection
    for (int i = 0; i < MAX_STEPS; ++i)
    {
        // Calculate the sample position in clip space
        vec3 sampleClipPos = clipPosition + stepVector * (float(i) + 0.1);

        // Convert clip space coordinates to integer texel coordinates
        ivec2 texCoords = ivec2(sampleClipPos.xy * texCoordScale);

        // Fetch the depth value from the depth texture
        float depthSample = texelFetch(depthtex0, texCoords, 0).r;

        // Check for intersection within the tolerance
        if (depthSample < sampleClipPos.z && abs(sampleClipPos.z - depthSample) < tolerance)
        {
            rtPos = vec3(sampleClipPos.xy, depthSample);
            break; // Intersection found; exit the loop
        }
    }

    return rtPos; // Return the intersection position or default value if not found
}

// 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;
}
vec3 SampleGGXVNDF3(vec3 Ve, float alpha, float U1, float U2)
{
    // Scale the view direction for the stretched ellipsoid
    vec3 Vh = normalize(vec3(alpha * Ve.x, alpha * Ve.y, Ve.z));

    // Create an orthonormal basis around Vh
    vec3 T1 = (Vh.z < 1.0) ? normalize(cross(vec3(0.0, 0.0, 1.0), Vh)) : vec3(1.0, 0.0, 0.0);
    vec3 T2 = cross(Vh, T1);

    // Projected area sampling
    float r = sqrt(U1);
    float phi = 2.0 * PI * U2;
    float t1 = r * cos(phi);
    float t2 = r * sin(phi);
    t2 = mix(sqrt(1.0 - t1 * t1), t2, 0.5 * (1.0 + Vh.z));

    // Map back to the hemisphere
    vec3 Nh = normalize(t1 * T1 + t2 * T2 + sqrt(max(0.0, 1.0 - t1 * t1 - t2 * t2)) * Vh);

    // Scale the hemisphere normal back to the ellipsoid
    return normalize(vec3(alpha * Nh.x, alpha * Nh.y, max(0.0, Nh.z)));
}
vec3 SampleGGXVNDF2(vec3 Ve, float alpha, float U2)
{
    // Scale the view direction for the stretched ellipsoid
    vec3 Vh = normalize(vec3(alpha * Ve.x, alpha * Ve.y, Ve.z));

    // Create an orthonormal basis around Vh
    vec3 T10 = (Vh.z < 1.0) ? normalize(cross(vec3(0.0, 0.0, 1.0), Vh)) : vec3(1.0, 0.0, 0.0);
    vec3 T20 = cross(Vh, T10);

    // Projected area sampling
    float r = 0.6667;
    float phi = 6.2831 * U2;
    float t1 = r * cos(phi);
    float t2 = r * sin(phi);
    t2 = mix(sqrt(1.0 - t1 * t1), t2, 0.75);

    // Map back to the hemisphere
    vec3 Nh = normalize(t1 * T10 + t2 * T20 + sqrt(max(0.0, 1.0 - t1 * t1 - t2 * t2)) * Vh);

    // Scale the hemisphere normal back to the ellipsoid
    return normalize(vec3(alpha * Nh.x, alpha * Nh.y, max(0.0, Nh.z)));
}
vec3 SampleBeckmannVNDF(vec3 Ve, float alpha, float U1)
{
    // Scale the view direction for the stretched ellipsoid
    vec3 Vh = normalize(vec3(alpha * Ve.x, alpha * Ve.y, Ve.z));

    // Create an orthonormal basis around Vh
    vec3 T10 = (Vh.z < 1.0) ? normalize(cross(vec3(0.0, 0.0, 1.0), Vh)) : vec3(1.0, 0.0, 0.0);
    vec3 T20 = cross(Vh, T10);

    // Beckmann sampling
    float theta = atan(sqrt(-log(U1)) * alpha); // Sampled angle based on Beckmann distribution
    float phi = 2.0 * 3.14159265359 * U1;       // Uniform azimuth angle
    float sinTheta = sin(theta);
    float cosTheta = cos(theta);

    // Convert polar coordinates to Cartesian
    float t1 = sinTheta * cos(phi);
    float t2 = sinTheta * sin(phi);

    // Map back to the hemisphere
    vec3 Nh = normalize(t1 * T10 + t2 * T20 + cosTheta * Vh);

    // Scale the hemisphere normal back to the ellipsoid
    return normalize(vec3(alpha * Nh.x, alpha * Nh.y, max(0.0, Nh.z)));
}

vec3 SampleGGXVNDF_SIMPLE(vec3 Ve, float alpha_x, float alpha_y, 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 = 0.707107;
    float phi = 2.0 * 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;
}

mat3 CoordBase(vec3 n)
{
    vec3 x, y;
    float nz = clamp(n.z, -0.999999, 1.0);
    float a = 1.0 / (1.0 + nz);
    float b = -n.x * n.y * a;
    x = vec3(1.0 - n.x * n.x * a, b, -n.x);
    y = vec3(b, 1.0 - n.y * n.y * a, -n.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;
}

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);
}
const bool colortex0MipmapEnabled = true;
const bool colortex4MipmapEnabled = true;
const bool colortex6MipmapEnabled = true;

////

vec3 GetN(int idx)
{
    switch (idx)
    {
    case 230:
        return vec3(2.9114, 2.9497, 2.5845);
    case 231:
        return vec3(0.18299, 0.42108, 1.3734);
    case 232:
        return vec3(1.3456, 0.96521, 0.61722);
    case 233:
        return vec3(3.1071, 3.1812, 2.3230);
    case 234:
        return vec3(0.27105, 0.67693, 1.3164);
    case 235:
        return vec3(1.9100, 1.8300, 1.4400);
    case 236:
        return vec3(2.3757, 2.0847, 1.8453);
    case 237:
        return vec3(0.15943, 0.14512, 0.13547);
    case 255:
        return vec3(1.0);
    default:
        return vec3(0.0);
    }
}

vec3 GetK(int idx)
{
    switch (idx)
    {
    case 230:
        return vec3(3.0893, 2.9318, 2.7670);
    case 231:
        return vec3(3.4242, 2.3459, 1.7704);
    case 232:
        return vec3(7.4746, 6.3995, 5.3031);
    case 233:
        return vec3(3.3314, 3.3291, 3.1350);
    case 234:
        return vec3(3.6092, 2.6248, 2.2921);
    case 235:
        return vec3(3.5100, 3.4000, 3.1800);
    case 236:
        return vec3(4.2655, 3.7153, 3.1365);
    case 237:
        return vec3(3.9291, 3.1900, 2.3808);
    case 255:
        return vec3(1.0);
    default:
        return vec3(1.0);
    }
}

vec3 sampleGGXVNDF(vec3 V_, float alpha_x, float alpha_y, float U1, float U2)
{
    // stretch view
    vec3 V = normalize(vec3(alpha_x * V_.x, alpha_y * V_.y, V_.z));
    // orthonormal basis
    vec3 T1 = (V.z < 0.9999) ? normalize(cross(V, vec3(0, 0, 1))) : vec3(1, 0, 0);
    vec3 T2 = cross(T1, V);
    // sample point with polar coordinates (r, phi)
    float a = 1.0 / (1.0 + V.z);
    float r = sqrt(U1);
    float phi = (U2 < a) ? U2 / a * 3.141592653589793 : 3.141592653589793 + (U2 - a) / (1.0 - a) * 3.141592653589793;
    float P1 = r * cos(phi);
    float P2 = r * sin(phi) * ((U2 < a) ? 1.0 : V.z);
    // compute normal
    vec3 N = P1 * T1 + P2 * T2 + sqrt(max(0.0, 1.0 - P1 * P1 - P2 * P2)) * V;
    // unstretch
    N = normalize(vec3(alpha_x * N.x, alpha_y * N.y, max(0.0, N.z)));
    return N;
}

#include "/lib/specular.glsl"
float g(float NdotL, float roughness)
{
    float k = (roughness * roughness) / 2.0;
    return NdotL / (NdotL * (1.0 - k) + k);
}

vec3 calculateSpecularLighting(
    vec3 normal,
    vec3 np3,
    vec3 albedo,
    float shading,
    float skylight,
    float roughness,
    vec3 f0,
    float porosity,
    vec3 fragpos,
    vec3 scene, bool hand, inout vec3 a, inout vec3 b)
{
    float noise = blueNoise(gl_FragCoord.xy);

    // Compute metal index from f0.y
    int metalidx = int((f0.y * 255));
    bool isMetal = (metalidx > 229);
    bool isNotHardcoded = (metalidx == 255);
    // roughness = 0.0;
    if (isMetal)
    {

#ifdef HARDCODED_METALS

        vec3 eta_k = GetK(metalidx);
        vec3 eta = GetN(metalidx);

        f0 = ((eta - vec3(1.0)) * (eta - vec3(1.0)) + (eta_k * eta_k)) / ((eta + vec3(1.0)) * (eta + vec3(1.0)) + (eta_k * eta_k));
        if (isNotHardcoded)
#endif
            f0 = albedo;
    }
    vec3 p3 = mat3(gbufferModelViewInverse) * fragpos + cameraPosition;
    vec2 uv = (p3.xz) * 0.5;
    float waterAmount = 0.0;

    //  waterAmount += length(texture(noisetex, mod(uv * 0.05, 1.0)).rgb);

    waterAmount = waterAmount * sqrt(skylight);

    float rainMult = sqrt(skylight) * wetness * (1.0 - pow(porosity, 2.0));
    roughness = mix(roughness, 0.0, rainMult * waterAmount);
    // roughness = 0.2;

    vec3 specTerm = shading * GGX2(normal, -np3, lightPos, roughness, f0) * pow(skylight, 3.0);

    // Indirect specular calculation
    vec3 indirectSpecular = vec3(0.0);
    vec3 fresnelDiffuse = vec3(0.0);
    int nSpecularSamples = int(mix(float(REFLECTION_SAMPLES), float(REFLECTION_SAMPLES_MAX), roughness));

    // Coordinate basis for sampling
    mat3 basis = CoordBase(normal);
    vec3 normSpaceView = normalize(-np3 * basis);

    // Texture size and LOD calculation
    // ivec2 texSize = textureSize(colortex4, 0); // Get the size of the base mipmap level (LOD 0)
    // float maxLod = float(log2(float(max(texSize.x, texSize.y))));
    float lod = mix(1.5, 8.0, roughness);
    vec3 reflectionColorScene = scene;

    if (!hand)
    {
        for (int i = 0; i < nSpecularSamples; i++)
        {
            // Generate random samples
            int seed = int(frameCounter) * nSpecularSamples + i;
            vec2 ij = fract(R2_samples(seed) + blueNoise4(gl_FragCoord.xy).rg);

            // Sample GGX VNDF
            // vec3 H = SampleGGXVNDF(normSpaceView, roughness, roughness, ij.x, ij.y);
            vec3 H = SampleGGXVNDF2(normSpaceView, roughness, ij.y);
            vec3 Ln = reflect(-normSpaceView, H);
            vec3 L = basis * Ln;

            // Fresnel term
            vec3 F = f0 + (1.0 - f0) * pow(clamp(1.0 + dot(-Ln, H), 0.0, 1.0), 5.0);

            // Geometry term
            float NoL = clamp(dot(normal, L), 0.0, 1.0);
            float G1 = g(NoL, roughness);

            // Ray contribution
            vec3 rayContrib = F * G1;
            vec3 reflectionColor = vec3(0.0);

            // Skip calculations if ray does not contribute much to the lighting
            if (luma(rayContrib) > 0.02)
            {

                //   reflectionColor = vec3(1.0, 0, 0);
                float rayQuality = 35.0 * sqrt(luma(rayContrib));

                if (rayQuality > 5.0)
                {

                    vec3 rtPos = RT5(mat3(gbufferModelView) * L, fragpos, fract(ij.x + ij.y), rayQuality);

                    if (rtPos.z < 1.0)
                    {
                        vec3 previousPosition = mat3(gbufferModelViewInverse) * toScreenSpace(rtPos) + gbufferModelViewInverse[3].xyz + cameraPosition - previousCameraPosition;
                        previousPosition = mat3(gbufferPreviousModelView) * previousPosition + gbufferPreviousModelView[3].xyz;
                        previousPosition.xy = projMAD(gbufferPreviousProjection, previousPosition).xy / -previousPosition.z * 0.5 + 0.5;

                        if (distance(texcoord, rtPos.xy) * 100 > 1)
                        {
                            if (previousPosition.x > 0.0 && previousPosition.y > 0.0 && previousPosition.x < 1.0 && previousPosition.y < 1.0)
                            {
                                reflectionColor = textureLod(colortex4, previousPosition.xy, lod).rgb;
                            }
                        }
                    }
                    else
                    {
                        float bounceAmount = float(L.y > 0.0) + clamp(-L.y * 0.1 + 0.1, 0.0, 1.0);
                        vec4 sky = skyFromTexLod(L, colortex6, lod);
                        sky.rgb = mix(sky.rgb, vec3(sky.a * sunLight), clamp(sky.a, 0, 1));

                        reflectionColor = mix(reflectionColorScene, sky.rgb, pow(skylight, 3.0) * bounceAmount);
                    }
                }

                indirectSpecular += reflectionColor * rayContrib;

                // Accumulate Fresnel diffuse term
                fresnelDiffuse += rayContrib;
            }
        }
    }
    else
    {

        vec3 H = SampleGGXVNDF(normSpaceView, roughness, roughness, 0.5, 0.5);
        vec3 Ln = reflect(-normSpaceView, H);
        vec3 L = basis * Ln;

        // Fresnel term
        vec3 F = f0 + (1.0 - f0) * pow(clamp(1.0 + dot(-Ln, H), 0.0, 1.0), 5.0);

        // Geometry term
        float NoL = clamp(dot(normal, L), 0.0, 1.0);
        float G1 = g(NoL, roughness);

        // Ray contribution
        vec3 rayContrib = F * G1;

        // Skip calculations if ray does not contribute much to the lighting

        vec4 sky = skyFromTexLod(L, colortex6, lod);
        sky.rgb = mix(sky.rgb, vec3(sky.a * sunLight), clamp(sky.a, 0, 1));
        vec3 reflectionColor = mix(reflectionColorScene, sky.rgb, skylight);

        indirectSpecular += reflectionColor * rayContrib;

        // Accumulate Fresnel diffuse term
        fresnelDiffuse += rayContrib;
    }
    // Average indirect specular and fresnel diffuse terms
    a = (indirectSpecular / float(nSpecularSamples)) + (specTerm * sunLight * SPECULAR_BOOST);
    b = (vec3(1.0) - fresnelDiffuse / float(nSpecularSamples));
    b = (isMetal ? vec3(0.0) : b);
    vec3 effectiveAlbedo = (metalidx < 229) ? vec3(1.0) : albedo;
    vec3 finalout = vec3((a + b * scene) * effectiveAlbedo);

    return finalout;
}

vec3 computeReflections(vec3 albedo, float lightmapY, float smoothness, vec3 f0, vec3 normal, vec3 fragpos, float materialsX, vec3 scene, bool hand, inout vec3 a, inout vec3 b)
{

    vec3 p3 = toWorldSpace(fragpos);
    float porosity = (materialsX * 255 <= 64.5) ? step(0.0, materialsX) : 0.0;
    float roughness = unpackRoughness(smoothness);
    float shading = texture(colortex1, texcoord).a;
    vec3 reflections = calculateSpecularLighting(normal, normalize(p3), albedo, shading, lightmapY, roughness, f0, porosity, fragpos, scene, hand, a, b);

    return vec3(reflections);
}
