Cycles: improve SSS Fresnel and retro-reflection in Principled BSDF

For details see the "Extending the Disney BRDF to a BSDF with Integrated
Subsurface Scattering" paper.

We split the diffuse BSDF into a lambertian and retro-reflection component.
The retro-reflection component is always handled as a BSDF, while the
lambertian component can be replaced by a BSSRDF.

For the BSSRDF case, we compute Fresnel separately at the entry and exit
points, which may have different normals. As the scattering radius decreases
this converges to the BSDF case.

A downside is that this increases noise for subsurface scattering in the
Principled BSDF, due to some samples going to the retro-reflection component.
However the previous logic (also in 2.93) was simple wrong, using a
non-sensical view direction vector at the exit point. We use an importance
sampling weight estimate for the retro-reflection to try to better balance
samples between the BSDF and BSSRDF.

Differential Revision: https://developer.blender.org/D12801

3 files changed, 84 insertions, 30 deletions

diff --git a/intern/cycles/kernel/closure/bsdf.h b/intern/cycles/kernel/closure/bsdf.h
index 4eb8bcae997..87aa6339f80 100644
--- a/intern/cycles/kernel/closure/bsdf.h
+++ b/intern/cycles/kernel/closure/bsdf.h
@@ -128,7 +128,6 @@ ccl_device_inline int bsdf_sample(const KernelGlobals *kg,
 
   switch (sc->type) {
     case CLOSURE_BSDF_DIFFUSE_ID:
-    case CLOSURE_BSDF_BSSRDF_ID:
       label = bsdf_diffuse_sample(sc,
                                   Ng,
                                   sd->I,
@@ -401,7 +400,6 @@ ccl_device_inline int bsdf_sample(const KernelGlobals *kg,
       break;
 #  ifdef __PRINCIPLED__
     case CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID:
-    case CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID:
       label = bsdf_principled_diffuse_sample(sc,
                                              Ng,
                                              sd->I,
@@ -481,7 +479,6 @@ ccl_device_inline
   if (!is_transmission) {
     switch (sc->type) {
       case CLOSURE_BSDF_DIFFUSE_ID:
-      case CLOSURE_BSDF_BSSRDF_ID:
         eval = bsdf_diffuse_eval_reflect(sc, sd->I, omega_in, pdf);
         break;
 #ifdef __SVM__
@@ -550,7 +547,6 @@ ccl_device_inline
         break;
 #  ifdef __PRINCIPLED__
       case CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID:
-      case CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID:
         eval = bsdf_principled_diffuse_eval_reflect(sc, sd->I, omega_in, pdf);
         break;
       case CLOSURE_BSDF_PRINCIPLED_SHEEN_ID:
@@ -576,7 +572,6 @@ ccl_device_inline
   else {
     switch (sc->type) {
       case CLOSURE_BSDF_DIFFUSE_ID:
-      case CLOSURE_BSDF_BSSRDF_ID:
         eval = bsdf_diffuse_eval_transmit(sc, sd->I, omega_in, pdf);
         break;
 #ifdef __SVM__
@@ -637,7 +632,6 @@ ccl_device_inline
         break;
 #  ifdef __PRINCIPLED__
       case CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID:
-      case CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID:
         eval = bsdf_principled_diffuse_eval_transmit(sc, sd->I, omega_in, pdf);
         break;
       case CLOSURE_BSDF_PRINCIPLED_SHEEN_ID:
diff --git a/intern/cycles/kernel/closure/bsdf_principled_diffuse.h b/intern/cycles/kernel/closure/bsdf_principled_diffuse.h
index 0d611f40096..04963ca1dc5 100644
--- a/intern/cycles/kernel/closure/bsdf_principled_diffuse.h
+++ b/intern/cycles/kernel/closure/bsdf_principled_diffuse.h
@@ -19,50 +19,98 @@
 /* DISNEY PRINCIPLED DIFFUSE BRDF
  *
  * Shading model by Brent Burley (Disney): "Physically Based Shading at Disney" (2012)
+ *
+ * "Extending the Disney BRDF to a BSDF with Integrated Subsurface Scattering" (2015)
+ * For the separation of retro-reflection, "2.3 Dielectric BRDF with integrated
+ * subsurface scattering"
  */
 
 #include "kernel/closure/bsdf_util.h"
 
 CCL_NAMESPACE_BEGIN
 
+enum PrincipledDiffuseBsdfComponents {
+  PRINCIPLED_DIFFUSE_FULL = 1,
+  PRINCIPLED_DIFFUSE_LAMBERT = 2,
+  PRINCIPLED_DIFFUSE_LAMBERT_EXIT = 4,
+  PRINCIPLED_DIFFUSE_RETRO_REFLECTION = 8,
+};
+
 typedef ccl_addr_space struct PrincipledDiffuseBsdf {
   SHADER_CLOSURE_BASE;
 
   float roughness;
+  int components;
 } PrincipledDiffuseBsdf;
 
 static_assert(sizeof(ShaderClosure) >= sizeof(PrincipledDiffuseBsdf),
               "PrincipledDiffuseBsdf is too large!");
 
-ccl_device float3 calculate_principled_diffuse_brdf(
+ccl_device int bsdf_principled_diffuse_setup(PrincipledDiffuseBsdf *bsdf)
+{
+  bsdf->type = CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID;
+  return SD_BSDF | SD_BSDF_HAS_EVAL;
+}
+
+ccl_device float3 bsdf_principled_diffuse_compute_brdf(
     const PrincipledDiffuseBsdf *bsdf, float3 N, float3 V, float3 L, float *pdf)
 {
-  float NdotL = dot(N, L);
+  const float NdotL = dot(N, L);
 
   if (NdotL <= 0) {
     return make_float3(0.0f, 0.0f, 0.0f);
   }
 
-  float NdotV = dot(N, V);
+  const float NdotV = dot(N, V);
+
+  const float FV = schlick_fresnel(NdotV);
+  const float FL = schlick_fresnel(NdotL);
+
+  float f = 0.0f;
 
-  /* H = normalize(L + V);  // Bisector of an angle between L and V.
-   * LH2 = 2 * dot(L, H)^2 = 2cos(x)^2 = cos(2x) + 1 = dot(L, V) + 1,
-   * half-angle x between L and V is at most 90 deg
-   */
-  float LH2 = dot(L, V) + 1;
+  /* Lambertian component. */
+  if (bsdf->components & (PRINCIPLED_DIFFUSE_FULL | PRINCIPLED_DIFFUSE_LAMBERT)) {
+    f += (1.0f - 0.5f * FV) * (1.0f - 0.5f * FL);
+  }
+  else if (bsdf->components & PRINCIPLED_DIFFUSE_LAMBERT_EXIT) {
+    f += (1.0f - 0.5f * FL);
+  }
 
-  float FL = schlick_fresnel(NdotL), FV = schlick_fresnel(NdotV);
-  const float Fd90 = 0.5f + LH2 * bsdf->roughness;
-  float Fd = (1.0f - FL + Fd90 * FL) * (1.0f - FV + Fd90 * FV);
+  /* Retro-reflection component. */
+  if (bsdf->components & (PRINCIPLED_DIFFUSE_FULL | PRINCIPLED_DIFFUSE_RETRO_REFLECTION)) {
+    /* H = normalize(L + V);  // Bisector of an angle between L and V
+     * LH2 = 2 * dot(L, H)^2 = 2cos(x)^2 = cos(2x) + 1 = dot(L, V) + 1,
+     * half-angle x between L and V is at most 90 deg. */
+    const float LH2 = dot(L, V) + 1;
+    const float RR = bsdf->roughness * LH2;
+    f += RR * (FL + FV + FL * FV * (RR - 1.0f));
+  }
 
-  float value = M_1_PI_F * NdotL * Fd;
+  float value = M_1_PI_F * NdotL * f;
 
   return make_float3(value, value, value);
 }
 
-ccl_device int bsdf_principled_diffuse_setup(PrincipledDiffuseBsdf *bsdf)
+/* Compute Fresnel at entry point, to be compbined with PRINCIPLED_DIFFUSE_LAMBERT_EXIT
+ * at the exit point to get the complete BSDF. */
+ccl_device_inline float bsdf_principled_diffuse_compute_entry_fresnel(const float NdotV)
+{
+  const float FV = schlick_fresnel(NdotV);
+  return (1.0f - 0.5f * FV);
+}
+
+/* Ad-hoc weight adjusment to avoid retro-reflection taking away half the
+ * samples from BSSRDF. */
+ccl_device_inline float bsdf_principled_diffuse_retro_reflection_sample_weight(
+    PrincipledDiffuseBsdf *bsdf, const float3 I)
+{
+  return bsdf->roughness * schlick_fresnel(dot(bsdf->N, I));
+}
+
+ccl_device int bsdf_principled_diffuse_setup(PrincipledDiffuseBsdf *bsdf, int components)
 {
   bsdf->type = CLOSURE_BSDF_PRINCIPLED_DIFFUSE_ID;
+  bsdf->components = components;
   return SD_BSDF | SD_BSDF_HAS_EVAL;
 }
 
@@ -79,7 +127,7 @@ ccl_device float3 bsdf_principled_diffuse_eval_reflect(const ShaderClosure *sc,
 
   if (dot(N, omega_in) > 0.0f) {
     *pdf = fmaxf(dot(N, omega_in), 0.0f) * M_1_PI_F;
-    return calculate_principled_diffuse_brdf(bsdf, N, V, L, pdf);
+    return bsdf_principled_diffuse_compute_brdf(bsdf, N, V, L, pdf);
   }
   else {
     *pdf = 0.0f;
@@ -115,7 +163,7 @@ ccl_device int bsdf_principled_diffuse_sample(const ShaderClosure *sc,
   sample_cos_hemisphere(N, randu, randv, omega_in, pdf);
 
   if (dot(Ng, *omega_in) > 0) {
-    *eval = calculate_principled_diffuse_brdf(bsdf, N, I, *omega_in, pdf);
+    *eval = bsdf_principled_diffuse_compute_brdf(bsdf, N, I, *omega_in, pdf);
 
 #ifdef __RAY_DIFFERENTIALS__
     // TODO: find a better approximation for the diffuse bounce
diff --git a/intern/cycles/kernel/closure/bssrdf.h b/intern/cycles/kernel/closure/bssrdf.h
index db183887018..d2f8af7910c 100644
--- a/intern/cycles/kernel/closure/bssrdf.h
+++ b/intern/cycles/kernel/closure/bssrdf.h
@@ -277,10 +277,27 @@ ccl_device_inline Bssrdf *bssrdf_alloc(ShaderData *sd, float3 weight)
 ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type, const float ior)
 {
   int flag = 0;
+
+  /* Add retro-reflection component as separate diffuse BSDF. */
+  if (bssrdf->roughness != FLT_MAX) {
+    PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf *)bsdf_alloc(
+        sd, sizeof(PrincipledDiffuseBsdf), bssrdf->weight);
+
+    if (bsdf) {
+      bsdf->N = bssrdf->N;
+      bsdf->roughness = bssrdf->roughness;
+      flag |= bsdf_principled_diffuse_setup(bsdf, PRINCIPLED_DIFFUSE_RETRO_REFLECTION);
+
+      /* Ad-hoc weight adjusment to avoid retro-reflection taking away half the
+       * samples from BSSRDF. */
+      bsdf->sample_weight *= bsdf_principled_diffuse_retro_reflection_sample_weight(bsdf, sd->I);
+    }
+  }
+
+  /* Verify if the radii are large enough to sample without precision issues. */
   int bssrdf_channels = 3;
   float3 diffuse_weight = make_float3(0.0f, 0.0f, 0.0f);
 
-  /* Verify if the radii are large enough to sample without precision issues. */
   if (bssrdf->radius.x < BSSRDF_MIN_RADIUS) {
     diffuse_weight.x = bssrdf->weight.x;
     bssrdf->weight.x = 0.0f;
@@ -304,17 +321,13 @@ ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type, co
     /* Add diffuse BSDF if any radius too small. */
 #ifdef __PRINCIPLED__
     if (bssrdf->roughness != FLT_MAX) {
-      float roughness = bssrdf->roughness;
-      float3 N = bssrdf->N;
-
       PrincipledDiffuseBsdf *bsdf = (PrincipledDiffuseBsdf *)bsdf_alloc(
           sd, sizeof(PrincipledDiffuseBsdf), diffuse_weight);
 
       if (bsdf) {
-        bsdf->type = CLOSURE_BSDF_BSSRDF_PRINCIPLED_ID;
-        bsdf->N = N;
-        bsdf->roughness = roughness;
-        flag |= bsdf_principled_diffuse_setup(bsdf);
+        bsdf->N = bssrdf->N;
+        bsdf->roughness = bssrdf->roughness;
+        flag |= bsdf_principled_diffuse_setup(bsdf, PRINCIPLED_DIFFUSE_LAMBERT);
       }
     }
     else
@@ -323,7 +336,6 @@ ccl_device int bssrdf_setup(ShaderData *sd, Bssrdf *bssrdf, ClosureType type, co
       DiffuseBsdf *bsdf = (DiffuseBsdf *)bsdf_alloc(sd, sizeof(DiffuseBsdf), diffuse_weight);
 
       if (bsdf) {
-        bsdf->type = CLOSURE_BSDF_BSSRDF_ID;
         bsdf->N = bssrdf->N;
         flag |= bsdf_diffuse_setup(bsdf);
       }