Cycles: Switch to energy-preserving multiscattering Oren-Nayar BSDF

This multiscattering term comes from the OpenPBR specification and nicely
preserves energy while correctly modeling increased saturation at high
roughness.

Preparation for adding a diffuse roughness option to the Principled BSDF.

To me, the difference in output and computation seems small enough to
not need an enum for the old behavior.

Note that this also switches sampling to cosine-weighted, in my tests this
gives lower noise. I also checked doing MIS between cosine and uniform,
using the A term as a weight for how often to use cosine (since that term
is Lambertian diffuse), but always using cosine was better.
A nice consequence of that is that you don't get a huge noise jump when
going from 0.0 to 0.01 roughness.

Pull Request: https://projects.blender.org/blender/blender/pulls/123345

intern/cycles/kernel/closure/bsdf_oren_nayar.h

@@ -12,10 +12,24 @@ typedef struct OrenNayarBsdf {
   float roughness;
   float a;
   float b;
+  Spectrum multiscatter_term;
 } OrenNayarBsdf;
 
 static_assert(sizeof(ShaderClosure) >= sizeof(OrenNayarBsdf), "OrenNayarBsdf is too large!");
 
+/* NOTE: This implements the improved Oren-Nayar model by Yasuhiro Fujii
+ * (https://mimosa-pudica.net/improved-oren-nayar.html), plus an
+ * energy-preserving multiscattering term based on the OpenPBR specification
+ * (https://academysoftwarefoundation.github.io/OpenPBR). */
+
+ccl_device_inline float bsdf_oren_nayar_G(const float cosTheta)
+{
+  const float sinTheta = sin_from_cos(cosTheta);
+  const float theta = safe_acosf(cosTheta);
+  return sinTheta * (theta - 2.0f / 3.0f - sinTheta * cosTheta) +
+         2.0f / 3.0f * (sinTheta / cosTheta) * (1.0f - sqr(sinTheta) * sinTheta);
+}
+
 ccl_device Spectrum bsdf_oren_nayar_get_intensity(ccl_private const ShaderClosure *sc,
                                                   float3 n,
                                                   float3 v,
@@ -28,22 +42,32 @@ ccl_device Spectrum bsdf_oren_nayar_get_intensity(ccl_private const ShaderClosur
 
   if (t > 0.0f)
     t /= max(nl, nv) + FLT_MIN;
-  float is = nl * (bsdf->a + bsdf->b * t);
-  return make_spectrum(is);
+
+  const float single_scatter = bsdf->a + bsdf->b * t;
+
+  const float El = bsdf->a * M_PI_F + bsdf->b * bsdf_oren_nayar_G(nl);
+  const Spectrum multi_scatter = bsdf->multiscatter_term * (1.0f - El);
+
+  return nl * (make_spectrum(single_scatter) + multi_scatter);
 }
 
-ccl_device int bsdf_oren_nayar_setup(ccl_private OrenNayarBsdf *bsdf)
+ccl_device int bsdf_oren_nayar_setup(ccl_private const ShaderData *sd,
+                                     ccl_private OrenNayarBsdf *bsdf,
+                                     const Spectrum albedo)
 {
-  float sigma = bsdf->roughness;
-
   bsdf->type = CLOSURE_BSDF_OREN_NAYAR_ID;
 
-  sigma = saturatef(sigma);
+  const float sigma = saturatef(bsdf->roughness);
+  bsdf->a = 1.0f / (M_PI_F + sigma * (M_PI_2_F - 2.0f / 3.0f));
+  bsdf->b = sigma * bsdf->a;
 
-  float div = 1.0f / (M_PI_F + ((3.0f * M_PI_F - 4.0f) / 6.0f) * sigma);
-
-  bsdf->a = 1.0f * div;
-  bsdf->b = sigma * div;
+  /* Compute energy compensation term (except for (1.0f - El) factor since it depends on wo). */
+  const float Eavg = bsdf->a * M_PI_F + ((M_2PI_F - 5.6f) / 3.0f) * bsdf->b;
+  const Spectrum Ems = M_1_PI_F * sqr(albedo) * (Eavg / (1.0f - Eavg)) /
+                       (one_spectrum() - albedo * (1.0f - Eavg));
+  const float nv = max(dot(bsdf->N, sd->wi), 0.0f);
+  const float Ev = bsdf->a * M_PI_F + bsdf->b * bsdf_oren_nayar_G(nv);
+  bsdf->multiscatter_term = Ems * (1.0f - Ev);
 
   return SD_BSDF | SD_BSDF_HAS_EVAL;
 }
@@ -54,8 +78,9 @@ ccl_device Spectrum bsdf_oren_nayar_eval(ccl_private const ShaderClosure *sc,
                                          ccl_private float *pdf)
 {
   ccl_private const OrenNayarBsdf *bsdf = (ccl_private const OrenNayarBsdf *)sc;
-  if (dot(bsdf->N, wo) > 0.0f) {
-    *pdf = 0.5f * M_1_PI_F;
+  const float cosNO = dot(bsdf->N, wo);
+  if (cosNO > 0.0f) {
+    *pdf = cosNO * M_1_PI_F;
     return bsdf_oren_nayar_get_intensity(sc, bsdf->N, wi, wo);
   }
   else {
@@ -73,7 +98,8 @@ ccl_device int bsdf_oren_nayar_sample(ccl_private const ShaderClosure *sc,
                                       ccl_private float *pdf)
 {
   ccl_private const OrenNayarBsdf *bsdf = (ccl_private const OrenNayarBsdf *)sc;
-  sample_uniform_hemisphere(bsdf->N, rand, wo, pdf);
+
+  sample_cos_hemisphere(bsdf->N, rand, wo, pdf);
 
   if (dot(Ng, *wo) > 0.0f) {
     *eval = bsdf_oren_nayar_get_intensity(sc, bsdf->N, wi, *wo);

intern/cycles/kernel/osl/closures_setup.h

@@ -95,6 +95,7 @@ ccl_device void osl_closure_diffuse_setup(KernelGlobals kg,
   sd->flag |= bsdf_diffuse_setup(bsdf);
 }
 
+/* Deprecated form, will be removed in OSL 2.0. */
 ccl_device void osl_closure_oren_nayar_setup(KernelGlobals kg,
                                              ccl_private ShaderData *sd,
                                              uint32_t path_flag,
@@ -115,7 +116,31 @@ ccl_device void osl_closure_oren_nayar_setup(KernelGlobals kg,
   bsdf->N = safe_normalize_fallback(closure->N, sd->N);
   bsdf->roughness = closure->roughness;
 
-  sd->flag |= bsdf_oren_nayar_setup(bsdf);
+  sd->flag |= bsdf_oren_nayar_setup(sd, bsdf, rgb_to_spectrum(weight));
+}
+
+ccl_device void osl_closure_oren_nayar_diffuse_bsdf_setup(
+    KernelGlobals kg,
+    ccl_private ShaderData *sd,
+    uint32_t path_flag,
+    float3 weight,
+    ccl_private const OrenNayarDiffuseBSDFClosure *closure,
+    float3 *layer_albedo)
+{
+  if (osl_closure_skip(kg, sd, path_flag, LABEL_DIFFUSE)) {
+    return;
+  }
+
+  ccl_private OrenNayarBsdf *bsdf = (ccl_private OrenNayarBsdf *)bsdf_alloc(
+      sd, sizeof(OrenNayarBsdf), rgb_to_spectrum(weight));
+  if (!bsdf) {
+    return;
+  }
+
+  bsdf->N = safe_normalize_fallback(closure->N, sd->N);
+  bsdf->roughness = closure->roughness;
+
+  sd->flag |= bsdf_oren_nayar_setup(sd, bsdf, rgb_to_spectrum(closure->albedo));
 }
 
 ccl_device void osl_closure_translucent_setup(KernelGlobals kg,

intern/cycles/kernel/osl/closures_template.h

@@ -19,11 +19,18 @@ OSL_CLOSURE_STRUCT_BEGIN(Diffuse, diffuse)
   OSL_CLOSURE_STRUCT_MEMBER(Diffuse, VECTOR, packed_float3, N, NULL)
 OSL_CLOSURE_STRUCT_END(Diffuse, diffuse)
 
+/* Deprecated form, will be removed in OSL 2.0. */
 OSL_CLOSURE_STRUCT_BEGIN(OrenNayar, oren_nayar)
   OSL_CLOSURE_STRUCT_MEMBER(OrenNayar, VECTOR, packed_float3, N, NULL)
   OSL_CLOSURE_STRUCT_MEMBER(OrenNayar, FLOAT, float, roughness, NULL)
 OSL_CLOSURE_STRUCT_END(OrenNayar, oren_nayar)
 
+OSL_CLOSURE_STRUCT_BEGIN(OrenNayarDiffuseBSDF, oren_nayar_diffuse_bsdf)
+  OSL_CLOSURE_STRUCT_MEMBER(OrenNayarDiffuseBSDF, VECTOR, packed_float3, N, NULL)
+  OSL_CLOSURE_STRUCT_MEMBER(OrenNayarDiffuseBSDF, VECTOR, packed_float3, albedo, NULL)
+  OSL_CLOSURE_STRUCT_MEMBER(OrenNayarDiffuseBSDF, FLOAT, float, roughness, NULL)
+OSL_CLOSURE_STRUCT_END(OrenNayarDiffuseBSDF, oren_nayar)
+
 OSL_CLOSURE_STRUCT_BEGIN(Translucent, translucent)
   OSL_CLOSURE_STRUCT_MEMBER(Translucent, VECTOR, packed_float3, N, NULL)
 OSL_CLOSURE_STRUCT_END(Translucent, translucent)

intern/cycles/kernel/osl/shaders/node_diffuse_bsdf.osl

@@ -12,5 +12,5 @@ shader node_diffuse_bsdf(color Color = 0.8,
   if (Roughness == 0.0)
     BSDF = Color * diffuse(Normal);
   else
-    BSDF = Color * oren_nayar(Normal, Roughness);
+    BSDF = Color * oren_nayar_diffuse_bsdf(Normal, Color, Roughness);
 }

intern/cycles/kernel/svm/closure.h

@@ -445,7 +445,8 @@ ccl_device
         }
         else {
           bsdf->roughness = roughness;
-          sd->flag |= bsdf_oren_nayar_setup(bsdf);
+          const Spectrum color = saturate(rgb_to_spectrum(stack_load_float3(stack, data_node.y)));
+          sd->flag |= bsdf_oren_nayar_setup(sd, bsdf, color);
         }
       }
       break;

intern/cycles/scene/shader_nodes.cpp

@@ -2613,7 +2613,7 @@ DiffuseBsdfNode::DiffuseBsdfNode() : BsdfNode(get_node_type())
 
 void DiffuseBsdfNode::compile(SVMCompiler &compiler)
 {
-  BsdfNode::compile(compiler, input("Roughness"), NULL);
+  BsdfNode::compile(compiler, input("Roughness"), nullptr, input("Color"));
 }
 
 void DiffuseBsdfNode::compile(OSLCompiler &compiler)