float smithGGXMasking(float ndotv, float a2) {
    float denomC = sqrt(a2 + (1.0 - a2) * ndotv * ndotv) + ndotv;

    return saturate(2.0f * ndotv / denomC);
}

float smithGGXMaskingShadowing(float ndotv, float ndotl, float a2) {
    float a = 2.0f * ndotl * ndotv;
    float denomA = ndotv * sqrt(a2 + (ndotv - ndotv * a2) * ndotv);
    float denomB = ndotl * sqrt(a2 + (ndotl - ndotl * a2) * ndotl);

    return saturate(a / (denomA + denomB));
}

float D_GGX(float ndoth, float roughnessSquared) {
    roughnessSquared = roughnessSquared < 1e-5 ? 0.0 : roughnessSquared;
	float p = (ndoth * roughnessSquared - ndoth) * ndoth + 1.0;
	return roughnessSquared / (pi * p * p);
}

vec3 ggxVndf(vec3 wo, vec2 roughness, vec2 U) {
    vec3 v = normalize(vec3(wo.x * roughness.x,
                            wo.y,
                            wo.z * roughness.y));

    vec3 t1 = (v.y < 0.999f) ? normalize(cross(v, vec3(0, 1, 0))) : vec3(1, 0, 0);
    vec3 t2 = cross(t1, v);

    float a = 1.0f / (1.0f + v.y);
    float r = sqrt(U.x);
    float phi = (U.y < a) ? (U.y / a) * pi 
                         : pi + (U.y - a) / (1.0f - a) * pi;
    float p1 = r * cos(phi);
    float p2 = r * sin(phi) * ((U.y < a) ? 1.0f : v.y);

    vec3 n = p1 * t1 + p2 * t2
             + sqrt(max(0.0f, 1.0f - p1 * p1 - p2 * p2)) * v;

    return normalize(vec3(roughness.x * n.x,
                          max(0.0f, n.y),
                          roughness.y * n.z));
}

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

vec3 specularBRDF(in float nDotL, in float nDotV, in float nDotH, in float vDotH) {
    nDotV = abs(nDotV);

    float G = smithGGXMaskingShadowing(nDotV, nDotL, sd.roughnessSquared);
    float D = D_GGX(nDotH, sd.roughnessSquared);
    vec3 F = fresnelFunction(vDotH, sd.n, sd.k);

    vec3 numerator = G * D * F;
    float denominator = 4.0 * nDotL * nDotV;

    vec3 specular = max(numerator / denominator, 0.0);

    return specular;
}