float unpackRoughness(float x)
{
    float r = 1.0 - x;
    return r * r;
}
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);
}
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 SampleBeckmannVNDF(vec3 Ve, float alpha, vec2 U)
{
    // Derive U1 and U2 from a single noise source
    float U1 = U.x;
    float U2 = U.y;

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

    // Beckmann sampling
    float tanTheta2 = -alpha * alpha * log(U1); // CDF inversion for Beckmann
    float theta = atan(sqrt(tanTheta2));
    float phi = 6.283185 * U2; // Azimuthal 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 * T1 + t2 * T2 + cosTheta * Vh);

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

float g(float NdotL, float roughness)
{
    float k = (roughness * roughness) / 2.0;
    float denom = NdotL * (1.0 - k) + k;
    denom = max(denom, 1e-5); // prevent division by zero
    return NdotL / denom;
}


vec3 SampleGGXVNDF2(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);

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

    float cosTheta = sqrt(max(0.0, 1.0 - t1 * t1 - t2 * t2));

    // 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 SampleVNDF2(vec3 Ve, float alpha)
{

    // Scale Ve's x and y components by alpha, keep z unchanged
    vec3 Vh = Ve * vec3(alpha, alpha, 1.0);

    // Compute the normalized tangent vector T10
    vec3 T10 = normalize(vec3(-Vh.y, Vh.x, 0.0));

    // Precompute the dither factor
    float dither = cos(6.2831 * R2_dither(gl_FragCoord.xy));

    // Compute the normal vector Nh with optimized coefficients
    vec3 Nh = (0.5 * Vh) + (T10 * 0.6667 * dither) + (0.1 * cross(Vh, T10));

    // Construct the final vector, applying alpha to x and y, and ensuring z is non-negative
    vec3 finalVec = vec3(alpha * Nh.xy, max(0.0, Nh.z));

    // Normalize the final vector before returning
    return normalize(finalVec);
}

vec3 SampleVNDF(vec3 np3, vec3 normal)
{
    // Ensure normal.z is within a stable range
    float nz = clamp(normal.z, -0.999999, 1.0);
    float a = 1.0 / (1.0 + nz);
    float b = -normal.x * normal.y * a;

    // Construct the basis matrix
    vec3 x = vec3(1.0 - normal.x * normal.x * a, b, -normal.x);
    vec3 y = vec3(b, 1.0 - normal.y * normal.y * a, -normal.y);
    mat3 basis = mat3(x, y, normal);

    // Transform and normalize the input vector
    vec3 Ve = normalize(-np3 * basis);
    vec3 Vh = normalize(vec3(0.0001 * Ve.xy, Ve.z));

    // Compute tangent vectors
    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);

    // Weighted combination for sampling
    vec3 weighted = Vh + T1 * -0.000043 + T2 * -0.25;
    vec3 H = normalize(vec3(0.0001 * weighted.xy, max(0.0, weighted.z)));

    // Compute reflected vector
    return basis * reflect(-Ve, H);
}

vec3 CalculateF0(vec3 f0, bool isMetal, int metalidx, vec3 albedo)
{
#ifdef HARDCODED_METALS
    if (isMetal)
    {

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

        vec3 A = (eta * eta) + (eta_k * eta_k) + vec3(1.0);
        vec3 N = A - (2.0 * eta);
        vec3 D = A + (2.0 * eta);
        f0 = N / D;

        if (metalidx == 255)
        {
            f0 = albedo;
        }
    }
#else
    if (metalidx == 255)
    {
        f0 = albedo;
    }
#endif

    return f0;
}
