     return mat3(tangent, bitangent, normal);
 }
 
-// NDF
-float D_GGX(float NdotH, float roughness)
-{
+struct MicrofacetDistributionSample
+{
+    float probability;
+    float cosTheta;
+    float sinTheta;
+    float phi;
+};
+
+float D_GGX(float NdotH, float roughness) {
+    float a = NdotH * roughness;
+    float k = roughness / (1.0 - NdotH * NdotH + a * a);
+    return k * k * (1.0 / MATH_PI);
+}
+
+// GGX microfacet distribution
+// https://www.cs.cornell.edu/~srm/publications/EGSR07-btdf.html
+MicrofacetDistributionSample GGX(vec2 xi, float roughness)
+{
+    MicrofacetDistributionSample ggx;
+
+    // evaluate sampling equations
     float alpha = roughness * roughness;
-
-    float alpha2 = alpha * alpha;
-
-    float divisor = NdotH * NdotH * (alpha2 - 1.0) + 1.0;
-
-    return alpha2 / (MATH_PI * divisor * divisor);
+    ggx.cosTheta = sqrt(1.0 - alpha * alpha * xi.x / (1.0 - xi.x));
+    ggx.sinTheta = alpha * sqrt(xi.x) / (sqrt(1.0 - xi.x) * ggx.cosTheta);
+    ggx.phi = 2.0 * MATH_PI * xi.y;
+
+    // evaluate GGX pdf (for half vector)
+    ggx.probability = D_GGX(ggx.cosTheta, roughness) * ggx.cosTheta; // drop sinTheta as we are evaluating for solid angle
+
+    return ggx;
 }
 
 // NDF

    mat3 TBN = generateTBN(N);
    vec3 direction = TBN * localSpaceDirection;

    if(u_distribution == cGGX || u_distribution == cCharlie)
    {
        pdf = PDF(direction, N, roughness);
    }

float PDF(vec3 H, vec3 N, float roughness)
{
    float NdotH = dot(N, H);

    if(u_distribution == cCharlie)
    {
        float D = D_Charlie(roughness, NdotH);
        return max(D, 0.0);
