float alpha_I(in vec3 wi, in vec2 alpha);

#include "heitz/slope/ggx.glsl"
#include "heitz/height/gaussian.glsl"
#include "heitz/slope.glsl"

vec3 generateMicrofacetNormal(in vec3 normal, in vec3 view, in float roughness) {
    vec3 facet = sampleD_wi(-Rotate(view, normal, vec3(0, 0, 1)), fract(randNext2F()), vec2(roughness));
    return Rotate(facet, vec3(0, 0, 1), normal);
}

#include "heitz/microsurface.glsl"

vec3 singleScatteringDiffuse(in vec3 albedo, in vec3 outgoingDirection, in vec3 incomingDirection, in vec3 normal, in vec3 wm, in float roughness) {
    incomingDirection = normalize(incomingDirection);
    outgoingDirection = normalize(outgoingDirection);
    vec2 alpha = vec2(roughness);

    if(dot(normal, -incomingDirection) < 0.0) {
        incomingDirection = reflect(incomingDirection, normal);
    }

    float lambdaI = smith_Lambda(-Rotate(incomingDirection, normal, vec3(0, 0, 1)), alpha);
    if(dot(normal, -incomingDirection) < 0.0) {
        lambdaI = 0.0;
    }
    float lambdaO = smith_Lambda( Rotate(outgoingDirection, normal, vec3(0, 0, 1)), alpha);
    float G2_given_G1 = (1.0 + lambdaI) / (1.0 + lambdaI + lambdaO);

    float value = 1.0 / pi * saturate(dot(wm, outgoingDirection)) * G2_given_G1;

    return albedo * value;
}

vec3 computeHeitzDiffuse(in vec3 albedo, in vec3 outgoingDirection, in vec3 incomingDirection, in vec3 normal, in vec3 n, in float roughness) {
    vec3 wm = generateMicrofacetNormal(normal, incomingDirection, roughness);

    vec3 f0 = vec3(IORTof0(n.r), IORTof0(n.g), IORTof0(n.b));

    vec3 fresnelV = 1.0 - fresnelFunction(saturate(dot(wm, -incomingDirection)), n, vec3(0.0), 0.0);
    vec3 fresnelL = 1.0 - fresnelFunction(saturate(dot(wm, outgoingDirection)), n, vec3(0.0), 0.0);

    #ifndef HEMISPHERE_INTEGRATION
        vec3 energyConservationFactor = 1 - (4 * sqrt(f0) + 5 * f0 * f0) / 9;

        fresnelL /= energyConservationFactor;
        fresnelV /= energyConservationFactor;
    #else
        vec3 energyConservationFactor = vec3(1.0);

        uint N = 32u;
        for(uint i = 0u; i < N; ++i) {
            vec2 hammersley = hammersley2d(i, N);
            float nDotV = 1.0 - hammersley.x; //This is equivalent to nDotV here.
            energyConservationFactor -= fresnelFunction(nDotV, n, vec3(0.0), 0.0) * (2.0 * nDotV) / N;
        }

        fresnelL /= energyConservationFactor;
        fresnelV /= energyConservationFactor;
    #endif

    vec3 single  = singleScatteringDiffuse(albedo, outgoingDirection, incomingDirection, normal, wm, roughness);
         single *= fresnelL;
         single *= fresnelV;

    vec3 lambertian  = albedo * saturate(dot(normal, outgoingDirection)) / pi;
         lambertian *= fresnelL;
         lambertian *= fresnelV;

    //This is most definitely not the correct fix, but it seems to work.
    if(dot(normal, outgoingDirection) <= 0.0) {
        return vec3(0.0);
    }

    return single;
}