diff --git a/intern/cycles/kernel/closure/bsdf_oren_nayar.h b/intern/cycles/kernel/closure/bsdf_oren_nayar.h
index 62c4955a397..c3c012043e7 100644
--- a/intern/cycles/kernel/closure/bsdf_oren_nayar.h
+++ b/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);
diff --git a/intern/cycles/kernel/osl/closures_setup.h b/intern/cycles/kernel/osl/closures_setup.h
index 970e4d8916c..c019753f55a 100644
--- a/intern/cycles/kernel/osl/closures_setup.h
+++ b/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,
diff --git a/intern/cycles/kernel/osl/closures_template.h b/intern/cycles/kernel/osl/closures_template.h
index a607700dfa3..11aac73ac27 100644
--- a/intern/cycles/kernel/osl/closures_template.h
+++ b/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)
diff --git a/intern/cycles/kernel/osl/shaders/node_diffuse_bsdf.osl b/intern/cycles/kernel/osl/shaders/node_diffuse_bsdf.osl
index 52db7f42d19..f6a073edb3e 100644
--- a/intern/cycles/kernel/osl/shaders/node_diffuse_bsdf.osl
+++ b/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);
 }
diff --git a/intern/cycles/kernel/svm/closure.h b/intern/cycles/kernel/svm/closure.h
index 24da8869a4d..cf8fc755c81 100644
--- a/intern/cycles/kernel/svm/closure.h
+++ b/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;
diff --git a/intern/cycles/scene/shader_nodes.cpp b/intern/cycles/scene/shader_nodes.cpp
index 6e14b27e1b7..c05f1de24aa 100644
--- a/intern/cycles/scene/shader_nodes.cpp
+++ b/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)