float G_1(in vec3 wi, in vec2 alpha, in float h0) {
    if(wi.z > 0.9999) {
        return 1.0;
    }
    if(wi.z <= 0.0) {
        return 0.0;
    }

    float C1_h0 = C1(h0);
    float lambda = smith_Lambda(wi, alpha);
    float value = pow(C1_h0, lambda);

    return value;
}

float G_1(in vec3 wi, in vec2 alpha) {
    if(wi.z > 0.9999) {
        return 1.0;
    }
    if(wi.z <= 0.0) {
        return 0.0;
    }

    float lambda = smith_Lambda(wi, alpha);
    float value = 1.0 / (1.0 + lambda);

    return value;
}

float sampleHeight(in vec3 wr, in vec2 alpha, in float hr, in float U) {
    if(wr.z > 0.9999) {
        return 1e35;
    }
    if(wr.z < -0.9999) {
        float value = invC1(U * C1(hr));
        return value;
    }
    if(abs(wr.z) < 0.0001) {
        return hr;
    }

    float G1 = G_1(wr, alpha, hr);

    if(U > 1.0 - G1) {
        return 1e35;
    }

    float h = invC1(
        C1(hr) / pow(1.0 - U, 1.0 / smith_Lambda(wr, alpha))
    );

    return h;
}

vec3 samplePhaseFunction(in vec3 direction, in vec3 normal, inout vec3 importance, in Bounce pathMaterial) {
    if(pathMaterial.isMetal > 0.9) {
        vec3 fresnel  = fresnelFunction(saturate(dot(normal, direction)), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);

        importance *= fresnel;
        #ifdef ALBEDO_TINT
            importance *= mix(pathMaterial.albedo, vec3(1.0), pathMaterial.albedoMetal);
        #endif

        return reflect(-direction, normal);
    } else {
        vec3 fresnel = fresnelFunction(dot(normal, direction), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);

        float fresnelLum = dot(fresnel, lumacoeff_rec709);
        float albedoLum = dot(pathMaterial.albedo, lumacoeff_rec709);
        float totalLum = albedoLum * (1.0 - fresnelLum) + fresnelLum; //I think this is named correctly? Not entirely sure.
        float specBounceProbability = fresnelLum / totalLum;

        bool specularBounce = specBounceProbability > randNextF();

        if(specularBounce) {
            importance /= specBounceProbability;
            
            importance *= fresnel;

            return reflect(-direction, normal);
        } else {
            importance /= 1.0 - specBounceProbability;

            #ifndef HEMISPHERE_INTEGRATION
                vec3 f03 = vec3(IORTof0(pathMaterial.n.r), IORTof0(pathMaterial.n.g), IORTof0(pathMaterial.n.b));

                vec3 energyConservationFactor = 1 - (4 * sqrt(f03) + 5 * f03 * f03) / 9;
            #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;
                }
            #endif

            importance *= 1.0 - fresnel;
            importance /= energyConservationFactor;

            vec3 newDirection = generateCosineVector(normal, randNext2F());
            //vec3 newFresnel = fresnelFunction(dot(normal, direction), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);

            importance *= pathMaterial.albedo;
            importance *= 1.0 - fresnelFunction(dot(normal, newDirection), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);

            return newDirection;
        }
    }
}

vec3 sampleMicrosurfaceBSDF(in vec3 direction, inout vec3 importance, in Bounce pathMaterial) {
    int scatteringOrder = 0;
    float hr = 1.0 + invC1(0.999);
    while(scatteringOrder <= MAX_BSDF_SCATTERING_EVENTS) {
        float U = randNextF();
        hr = sampleHeight(direction, vec2(pathMaterial.roughness), hr, U);

        if(hr == 1e35) {
            break;
        }

        scatteringOrder++;

        vec3 wm = sampleD_wi(-direction, fract(randNext2F()), vec2(pathMaterial.roughness));

        direction = samplePhaseFunction(-direction, wm, importance, pathMaterial);

        if(hr != hr || direction.z != direction.z) {
            return vec3(0, 0, 1);
        }
    }

    if(scatteringOrder > MAX_BSDF_SCATTERING_EVENTS) {
        importance = vec3(0.0);
    }

    return direction;
}

vec3 estimateBSDFPhase(in vec3 direction1, in vec3 direction2, in vec2 alpha, in Bounce pathMaterial) {
    vec3 halfway = normalize(direction1 + direction2);

    if(pathMaterial.isMetal > 0.9) {
        float vDotH = dot(direction1, halfway);
        vec3 specular  = fresnelFunction(vDotH, pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);
             specular *= D_wi(direction1, halfway, alpha);
             specular /= 4.0 * vDotH;

        #ifdef ALBEDO_TINT
            specular *= mix(pathMaterial.albedo, vec3(1.0), pathMaterial.albedoMetal);
        #endif

        return specular;
    } else {
        #ifndef HEMISPHERE_INTEGRATION
            vec3 f03 = vec3(IORTof0(pathMaterial.n.r), IORTof0(pathMaterial.n.g), IORTof0(pathMaterial.n.b));

            vec3 energyConservationFactor = 1 - (4 * sqrt(f03) + 5 * f03 * f03) / 9;
        #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;
            }
        #endif

        vec3 wm = sampleD_wi(direction1, fract(randNext2F()), alpha);

        vec3 diffuse  = pathMaterial.albedo / pi;
             diffuse *= saturate(dot(direction2, wm));
             diffuse *= 1.0 - fresnelFunction(          dot(wm, direction1), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);
             diffuse *= 1.0 - fresnelFunction(saturate(dot(wm, direction2)), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);
             diffuse /= energyConservationFactor;


        float vDotH = dot(direction1, halfway);
        vec3 specular  = fresnelFunction(saturate(dot(wm, direction1)), pathMaterial.n, pathMaterial.k, pathMaterial.albedoMetal);
             specular *= D_wi(direction1, halfway, alpha);
             specular /= 4.0 * vDotH;

        return specular + diffuse;
    }
}

vec3 estimateMicrosurfaceBSDF(in vec3 direction1, in vec3 direction2, in int bounceID, in Bounce pathMaterial) {
    float hr = 1.0 + invC1(0.999);

    vec2 alpha = vec2(pathMaterial.roughness);

    #ifdef CELESTIAL_PATH_ROUGHNESS_BOUNCE_INCREASE
        if(bounceID > 1) {
            alpha *= 10.0;
            alpha = saturate(alpha);
        }
    #endif

    int scatteringOrder = 0;
    vec3 importance = vec3(1.0);
    vec3 estimate = vec3(0.0);
    while(scatteringOrder <= MAX_BSDF_SCATTERING_EVENTS) {
        float U = randNextF();
        hr = sampleHeight(direction1, alpha, hr, U);

        if(hr == 1e35) {
            break;
        }

        scatteringOrder++;

        vec3 phase = estimateBSDFPhase(-direction1, direction2, alpha, pathMaterial);
        estimate += max(importance, 0.0) * max(phase, 0.0) * G_1(direction2, alpha, hr);

        vec3 wm = sampleD_wi(-direction1, fract(randNext2F()), alpha);

        direction1 = samplePhaseFunction(-direction1, wm, importance, pathMaterial);

        if(hr != hr || direction1.z != direction1.z) {
            return vec3(0, 0, 1);
        }
    }

    return estimate;
}