EEVEE: GGX: Use distribution of visible normal for sampling

This changes the sampling routine to use the method described in
"A Simpler and Exact Sampling Routine for the GGXDistribution of Visible Normals"
by Eric Heitz.
http://jcgt.org/published/0007/04/01/slides.pdf

This avoids generating bad rays and thus improve noise level in screen-
space reflections / refraction.

8 files changed, 104 insertions, 60 deletions

diff --git a/source/blender/draw/engines/eevee/eevee_lightprobes.c b/source/blender/draw/engines/eevee/eevee_lightprobes.c
index 0ce94b8f1b1..3e1d4a8aaa6 100644
--- a/source/blender/draw/engines/eevee/eevee_lightprobes.c
+++ b/source/blender/draw/engines/eevee/eevee_lightprobes.c
@@ -245,7 +245,7 @@ void EEVEE_lightbake_cache_init(EEVEE_ViewLayerData *sldata,
     DRW_shgroup_uniform_float(grp, "intensityFac", &pinfo->intensity_fac, 1);
     DRW_shgroup_uniform_float(grp, "sampleCount", &pinfo->samples_len, 1);
     DRW_shgroup_uniform_float(grp, "invSampleCount", &pinfo->samples_len_inv, 1);
-    DRW_shgroup_uniform_float(grp, "roughnessSquared", &pinfo->roughness, 1);
+    DRW_shgroup_uniform_float(grp, "roughness", &pinfo->roughness, 1);
     DRW_shgroup_uniform_float(grp, "lodFactor", &pinfo->lodfactor, 1);
     DRW_shgroup_uniform_float(grp, "lodMax", &pinfo->lod_rt_max, 1);
     DRW_shgroup_uniform_float(grp, "texelSize", &pinfo->texel_size, 1);
@@ -1045,8 +1045,8 @@ void EEVEE_lightbake_filter_glossy(EEVEE_ViewLayerData *sldata,
     /* Disney Roughness */
     pinfo->roughness = square_f(pinfo->roughness);
     /* Distribute Roughness across lod more evenly */
-    pinfo->roughness = square_f(square_f(pinfo->roughness));
-    CLAMP(pinfo->roughness, 1e-8f, 0.99999f); /* Avoid artifacts */
+    pinfo->roughness = square_f(pinfo->roughness);
+    CLAMP(pinfo->roughness, 1e-4f, 0.9999f); /* Avoid artifacts */
 
 #if 1 /* Variable Sample count and bias (fast) */
     switch (i) {
diff --git a/source/blender/draw/engines/eevee/shaders/bsdf_common_lib.glsl b/source/blender/draw/engines/eevee/shaders/bsdf_common_lib.glsl
index 5f0fde134d1..c8eaa06094e 100644
--- a/source/blender/draw/engines/eevee/shaders/bsdf_common_lib.glsl
+++ b/source/blender/draw/engines/eevee/shaders/bsdf_common_lib.glsl
@@ -102,7 +102,7 @@ float D_ggx_opti(float NH, float a2)
   return M_PI * tmp * tmp; /* Doing RCP and mul a2 at the end */
 }
 
-float G1_Smith_GGX(float NX, float a2)
+float G1_Smith_GGX_opti(float NX, float a2)
 {
   /* Using Brian Karis approach and refactoring by NX/NX
    * this way the (2*NL)*(2*NV) in G = G1(V) * G1(L) gets canceled by the brdf denominator 4*NL*NV
@@ -122,7 +122,7 @@ float bsdf_ggx(vec3 N, vec3 L, vec3 V, float roughness)
   float NL = max(dot(N, L), 1e-8);
   float NV = max(dot(N, V), 1e-8);
 
-  float G = G1_Smith_GGX(NV, a2) * G1_Smith_GGX(NL, a2); /* Doing RCP at the end */
+  float G = G1_Smith_GGX_opti(NV, a2) * G1_Smith_GGX_opti(NL, a2); /* Doing RCP at the end */
   float D = D_ggx_opti(NH, a2);
 
   /* Denominator is canceled by G1_Smith */
diff --git a/source/blender/draw/engines/eevee/shaders/bsdf_lut_frag.glsl b/source/blender/draw/engines/eevee/shaders/bsdf_lut_frag.glsl
index 46ea8b747c8..4c1544654c1 100644
--- a/source/blender/draw/engines/eevee/shaders/bsdf_lut_frag.glsl
+++ b/source/blender/draw/engines/eevee/shaders/bsdf_lut_frag.glsl
@@ -25,7 +25,8 @@ void main()
       vec3 Xi = (vec3(i, j, 0.0) + 0.5) / sampleCount;
       Xi.yz = vec2(cos(Xi.y * M_2PI), sin(Xi.y * M_2PI));
 
-      vec3 H = sample_ggx(Xi, a2); /* Microfacet normal */
+      /* Microfacet normal */
+      vec3 H = sample_ggx(Xi, a, V);
       vec3 L = -reflect(V, H);
       float NL = L.z;
 
@@ -33,9 +34,10 @@ void main()
         float NH = max(H.z, 0.0);
         float VH = max(dot(V, H), 0.0);
 
-        float G1_v = G1_Smith_GGX(NV, a2);
-        float G1_l = G1_Smith_GGX(NL, a2);
-        float G_smith = 4.0 * NV * NL / (G1_v * G1_l); /* See G1_Smith_GGX for explanations. */
+        float G1_v = G1_Smith_GGX_opti(NV, a2);
+        float G1_l = G1_Smith_GGX_opti(NL, a2);
+        /* See G1_Smith_GGX_opti for explanations. */
+        float G_smith = 4.0 * NV * NL / (G1_v * G1_l);
 
         float brdf = (G_smith * VH) / (NH * NV);
 
diff --git a/source/blender/draw/engines/eevee/shaders/bsdf_sampling_lib.glsl b/source/blender/draw/engines/eevee/shaders/bsdf_sampling_lib.glsl
index bbc79a2d05b..004d884dc75 100644
--- a/source/blender/draw/engines/eevee/shaders/bsdf_sampling_lib.glsl
+++ b/source/blender/draw/engines/eevee/shaders/bsdf_sampling_lib.glsl
@@ -1,4 +1,5 @@
 
+#pragma BLENDER_REQUIRE(common_math_geom_lib.glsl)
 #pragma BLENDER_REQUIRE(bsdf_common_lib.glsl)
 
 uniform sampler1D texHammersley;
@@ -8,6 +9,11 @@ vec3 tangent_to_world(vec3 vector, vec3 N, vec3 T, vec3 B)
   return T * vector.x + B * vector.y + N * vector.z;
 }
 
+vec3 world_to_tangent(vec3 vector, vec3 N, vec3 T, vec3 B)
+{
+  return vec3(dot(T, vector), dot(B, vector), dot(N, vector));
+}
+
 #ifdef HAMMERSLEY_SIZE
 vec3 hammersley_3d(float i, float invsamplenbr)
 {
@@ -22,9 +28,18 @@ vec3 hammersley_3d(float i, float invsamplenbr)
 
 /* -------------- BSDFS -------------- */
 
-float pdf_ggx_reflect(float NH, float a2)
+#define USE_VISIBLE_NORMAL 1
+
+float pdf_ggx_reflect(float NH, float NV, float VH, float alpha)
 {
+  float a2 = sqr(alpha);
+#if USE_VISIBLE_NORMAL
+  float D = a2 / D_ggx_opti(NH, a2);
+  float G1 = NV * 2.0 / G1_Smith_GGX_opti(NV, a2);
+  return G1 * VH * D / NV;
+#else
   return NH * a2 / D_ggx_opti(NH, a2);
+#endif
 }
 
 float pdf_hemisphere()
@@ -32,22 +47,50 @@ float pdf_hemisphere()
   return 0.5 * M_1_PI;
 }
 
-vec3 sample_ggx(vec3 rand, float a2)
+vec3 sample_ggx(vec3 rand, float alpha, vec3 Vt)
 {
+#if USE_VISIBLE_NORMAL
+  /* From:
+   * "A Simpler and Exact Sampling Routine for the GGXDistribution of Visible Normals"
+   * by Eric Heitz.
+   * http://jcgt.org/published/0007/04/01/slides.pdf
+   * View vector is expected to be in tangent space. */
+
+  /* Stretch view. */
+  vec3 Th, Bh, Vh = normalize(vec3(alpha * Vt.xy, Vt.z));
+  make_orthonormal_basis(Vh, Th, Bh);
+  /* Sample point with polar coordinates (r, phi). */
+  float r = sqrt(rand.x);
+  float x = r * rand.y;
+  float y = r * rand.z;
+  float s = 0.5 * (1.0 + Vh.z);
+  y = (1.0 - s) * sqrt(1.0 - x * x) + s * y;
+  float z = sqrt(saturate(1.0 - x * x - y * y));
+  /* Compute normal. */
+  vec3 Hh = x * Th + y * Bh + z * Vh;
+  /* Unstretch. */
+  vec3 Ht = normalize(vec3(alpha * Hh.xy, saturate(Hh.z)));
+  /* Microfacet Normal. */
+  return Ht;
+#else
   /* Theta is the cone angle. */
-  float z = sqrt((1.0 - rand.x) / (1.0 + a2 * rand.x - rand.x)); /* cos theta */
-  float r = sqrt(max(0.0, 1.0 - z * z));                         /* sin theta */
+  float z = sqrt((1.0 - rand.x) / (1.0 + sqr(alpha) * rand.x - rand.x)); /* cos theta */
+  float r = sqrt(max(0.0, 1.0 - z * z));                                 /* sin theta */
   float x = r * rand.y;
   float y = r * rand.z;
-
   /* Microfacet Normal */
   return vec3(x, y, z);
+#endif
 }
 
-vec3 sample_ggx(vec3 rand, float a2, vec3 N, vec3 T, vec3 B, out float NH)
+vec3 sample_ggx(vec3 rand, float alpha, vec3 V, vec3 N, vec3 T, vec3 B, out float pdf)
 {
-  vec3 Ht = sample_ggx(rand, a2);
-  NH = Ht.z;
+  vec3 Vt = world_to_tangent(V, N, T, B);
+  vec3 Ht = sample_ggx(rand, alpha, Vt);
+  float NH = saturate(Ht.z);
+  float NV = saturate(Vt.z);
+  float VH = saturate(dot(Vt, Ht));
+  pdf = pdf_ggx_reflect(NH, NV, VH, alpha);
   return tangent_to_world(Ht, N, T, B);
 }
 
@@ -69,18 +112,23 @@ vec3 sample_hemisphere(vec3 rand, vec3 N, vec3 T, vec3 B)
 }
 
 #ifdef HAMMERSLEY_SIZE
-vec3 sample_ggx(float nsample, float inv_sample_count, float a2, vec3 N, vec3 T, vec3 B)
+vec3 sample_ggx(float nsample,
+                float inv_sample_count,
+                float alpha,
+                vec3 V,
+                vec3 N,
+                vec3 T,
+                vec3 B,
+                out float pdf)
 {
   vec3 Xi = hammersley_3d(nsample, inv_sample_count);
-  vec3 Ht = sample_ggx(Xi, a2);
-  return tangent_to_world(Ht, N, T, B);
+  return sample_ggx(Xi, alpha, V, N, T, B, pdf);
 }
 
 vec3 sample_hemisphere(float nsample, float inv_sample_count, vec3 N, vec3 T, vec3 B)
 {
   vec3 Xi = hammersley_3d(nsample, inv_sample_count);
-  vec3 Ht = sample_hemisphere(Xi);
-  return tangent_to_world(Ht, N, T, B);
+  return sample_hemisphere(Xi, N, T, B);
 }
 
 vec3 sample_cone(float nsample, float inv_sample_count, float angle, vec3 N, vec3 T, vec3 B)
diff --git a/source/blender/draw/engines/eevee/shaders/btdf_lut_frag.glsl b/source/blender/draw/engines/eevee/shaders/btdf_lut_frag.glsl
index 2ffe23a9197..2f1298e2707 100644
--- a/source/blender/draw/engines/eevee/shaders/btdf_lut_frag.glsl
+++ b/source/blender/draw/engines/eevee/shaders/btdf_lut_frag.glsl
@@ -37,7 +37,7 @@ void main()
       Xi.yz = vec2(cos(Xi.y * M_2PI), sin(Xi.y * M_2PI));
 
       /* Microfacet normal. */
-      vec3 H = sample_ggx(Xi, a2);
+      vec3 H = sample_ggx(Xi, a2, V);
 
       float VH = dot(V, H);
 
@@ -59,7 +59,7 @@ void main()
         float LH = dot(L, H);
 
         /* Balancing the adjustments made in G1_Smith. */
-        float G1_l = NL * 2.0 / G1_Smith_GGX(NL, a2);
+        float G1_l = NL * 2.0 / G1_Smith_GGX_opti(NL, a2);
 
         /* btdf = abs(VH*LH) * (ior*ior) * D * G(V) * G(L) / (Ht2 * NV)
          * pdf = (VH * abs(LH)) * (ior*ior) * D * G(V) / (Ht2 * NV) */
diff --git a/source/blender/draw/engines/eevee/shaders/effect_ssr_frag.glsl b/source/blender/draw/engines/eevee/shaders/effect_ssr_frag.glsl
index a9139b7a146..11aa85f2f2d 100644
--- a/source/blender/draw/engines/eevee/shaders/effect_ssr_frag.glsl
+++ b/source/blender/draw/engines/eevee/shaders/effect_ssr_frag.glsl
@@ -64,27 +64,22 @@ uniform sampler2D specroughBuffer;
 
 layout(location = 0) out vec4 hitData;
 
-void do_planar_ssr(
-    int index, vec3 V, vec3 N, vec3 T, vec3 B, vec3 viewPlaneNormal, vec3 vP, float a2, vec4 rand)
+void do_planar_ssr(int index,
+                   vec3 V,
+                   vec3 N,
+                   vec3 T,
+                   vec3 B,
+                   vec3 viewPlaneNormal,
+                   vec3 vP,
+                   float alpha,
+                   vec4 rand)
 {
-  float NH;
+  float pdf;
   /* Microfacet normal */
-  vec3 H = sample_ggx(rand.xzw, a2, N, T, B, NH);
-  float pdf = pdf_ggx_reflect(NH, a2);
-
+  vec3 H = sample_ggx(rand.xzw, alpha, V, N, T, B, pdf);
   vec3 R = reflect(-V, H);
   R = reflect(R, viewPlaneNormal);
 
-  /* If ray is bad (i.e. going below the plane) regenerate. */
-  if (dot(R, viewPlaneNormal) > 0.0) {
-    /* Microfacet normal */
-    vec3 H = sample_ggx(rand.xzw * vec3(1.0, -1.0, -1.0), a2, N, T, B, NH);
-    pdf = pdf_ggx_reflect(NH, a2);
-
-    R = reflect(-V, H);
-    R = reflect(R, viewPlaneNormal);
-  }
-
   Ray ray;
   ray.origin = vP;
   ray.direction = R * 1e16;
@@ -92,7 +87,7 @@ void do_planar_ssr(
   RayTraceParameters params;
   params.jitter = rand.y;
   params.trace_quality = ssrQuality;
-  params.roughness = a2;
+  params.roughness = alpha * alpha;
 
   HitData data;
   data.is_planar = true;
@@ -104,26 +99,26 @@ void do_planar_ssr(
   hitData = encode_hit_data(data);
 }
 
-void do_ssr(vec3 V, vec3 N, vec3 T, vec3 B, vec3 vP, float a2, vec4 rand)
+void do_ssr(vec3 V, vec3 N, vec3 T, vec3 B, vec3 vP, float alpha, vec4 rand)
 {
-  float NH;
+  float pdf;
   /* Microfacet normal */
-  vec3 H = sample_ggx(rand.xzw, a2, N, T, B, NH);
+  vec3 H = sample_ggx(rand.xzw, alpha, V, N, T, B, pdf);
   vec3 R = reflect(-V, H);
 
   Ray ray;
-  ray.origin = vP;
+  ray.origin = vP + N * 1e-4;
   ray.direction = R * 1e16;
 
   RayTraceParameters params;
   params.thickness = ssrThickness;
   params.jitter = rand.y;
   params.trace_quality = ssrQuality;
-  params.roughness = a2;
+  params.roughness = alpha * alpha;
 
   HitData data;
   data.is_planar = true;
-  data.ray_pdf_inv = safe_rcp(pdf_ggx_reflect(NH, a2));
+  data.ray_pdf_inv = safe_rcp(pdf);
   data.is_hit = raytrace(ray, params, true, data.hit_dir);
   data.hit_dir = get_view_space_from_depth(data.hit_dir.xy, data.hit_dir.z);
   data.hit_dir -= ray.origin;
@@ -170,8 +165,7 @@ void main()
   }
 
   float roughness = speccol_roughness.a;
-  float roughnessSquared = max(1e-3, roughness * roughness);
-  float a2 = roughnessSquared * roughnessSquared;
+  float alpha = max(1e-3, roughness * roughness);
 
   /* Early out */
   if (roughness > ssrMaxRoughness + 0.2) {
@@ -204,12 +198,12 @@ void main()
       tracePosition = transform_point(ViewMatrix, tracePosition);
       vec3 viewPlaneNormal = transform_direction(ViewMatrix, pd.pl_normal);
 
-      do_planar_ssr(i, vV, vN, vT, vB, viewPlaneNormal, tracePosition, a2, rand);
+      do_planar_ssr(i, vV, vN, vT, vB, viewPlaneNormal, tracePosition, alpha, rand);
       return;
     }
   }
 
-  do_ssr(vV, vN, vT, vB, vP, a2, rand);
+  do_ssr(vV, vN, vT, vB, vP, alpha, rand);
 }
 
 #else /* STEP_RESOLVE */
diff --git a/source/blender/draw/engines/eevee/shaders/lightprobe_filter_glossy_frag.glsl b/source/blender/draw/engines/eevee/shaders/lightprobe_filter_glossy_frag.glsl
index 35fdbcb715f..99cbf2839ad 100644
--- a/source/blender/draw/engines/eevee/shaders/lightprobe_filter_glossy_frag.glsl
+++ b/source/blender/draw/engines/eevee/shaders/lightprobe_filter_glossy_frag.glsl
@@ -3,7 +3,7 @@
 #pragma BLENDER_REQUIRE(common_math_geom_lib.glsl)
 
 uniform samplerCube probeHdr;
-uniform float roughnessSquared;
+uniform float roughness;
 uniform float texelSize;
 uniform float lodFactor;
 uniform float lodMax;
@@ -52,7 +52,9 @@ void main()
   float weight = 0.0;
   vec3 out_radiance = vec3(0.0);
   for (float i = 0; i < sampleCount; i++) {
-    vec3 H = sample_ggx(i, invSampleCount, roughnessSquared, N, T, B); /* Microfacet normal */
+    float pdf;
+    /* Microfacet normal */
+    vec3 H = sample_ggx(i, invSampleCount, roughness, V, N, T, B, pdf);
     vec3 L = -reflect(V, H);
     float NL = dot(N, L);
 
@@ -62,7 +64,6 @@ void main()
       /* Coarse Approximation of the mapping distortion
        * Unit Sphere -> Cubemap Face */
       const float dist = 4.0 * M_PI / 6.0;
-      float pdf = pdf_ggx_reflect(NH, roughnessSquared);
       /* http://http.developer.nvidia.com/GPUGems3/gpugems3_ch20.html : Equation 13 */
       float lod = clamp(lodFactor - 0.5 * log2(pdf * dist), 0.0, lodMax);
 
diff --git a/source/blender/draw/engines/eevee/shaders/ssr_lib.glsl b/source/blender/draw/engines/eevee/shaders/ssr_lib.glsl
index 37624fc31ea..612e95832e4 100644
--- a/source/blender/draw/engines/eevee/shaders/ssr_lib.glsl
+++ b/source/blender/draw/engines/eevee/shaders/ssr_lib.glsl
@@ -15,21 +15,20 @@ uniform float refractionDepth;
 
 vec4 screen_space_refraction(vec3 vP, vec3 N, vec3 V, float ior, float roughnessSquared, vec4 rand)
 {
-  float a2 = max(5e-6, roughnessSquared * roughnessSquared);
+  float alpha = max(0.002, roughnessSquared);
 
   /* Importance sampling bias */
   rand.x = mix(rand.x, 0.0, BTDF_BIAS);
 
   vec3 T, B;
-  float NH;
   make_orthonormal_basis(N, T, B);
-  vec3 H = sample_ggx(rand.xzw, a2, N, T, B, NH); /* Microfacet normal */
-  float pdf = pdf_ggx_reflect(NH, a2);
+  float pdf;
+  /* Microfacet normal */
+  vec3 H = sample_ggx(rand.xzw, alpha, V, N, T, B, pdf);
 
   /* If ray is bad (i.e. going below the plane) regenerate. */
   if (F_eta(ior, dot(H, V)) < 1.0) {
-    H = sample_ggx(rand.xzw * vec3(1.0, -1.0, -1.0), a2, N, T, B, NH); /* Microfacet normal */
-    pdf = pdf_ggx_reflect(NH, a2);
+    H = sample_ggx(rand.xzw * vec3(1.0, -1.0, -1.0), alpha, V, N, T, B, pdf);
   }
 
   vec3 vV = viewCameraVec(vP);