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 + (1.0 - a2) * ndotl * ndotl);
    float denomB = ndotl * sqrt(a2 + (1.0 - a2) * ndotv * ndotv);

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

float D_GGX(float ndoth, float alpha2) {
    //alpha2 = alpha2 < 1e-5 ? 0.0 : alpha2;
	float p = (ndoth * alpha2 - ndoth) * ndoth + 1.0;
	return alpha2 / (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 is_ggx(vec3 normal, vec4 noise, float alpha2) {
    noise.xyz = normalize(cross(normal, noise.xyz * 2.0 - 1.0));
    noise.w *= 0.9;
    return normalize(noise.xyz * sqrt(alpha2 * noise.w / (1.0 - noise.w)) + normal);
}

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