Implement new Oren-Nayar model from OpenPBR (#1817)

* Implement new EON mode for Oren-Nayar closures (via an opt-in keyword arg for backwards compatibility)
* Add test for energy conserving Oren-Nayar (on half the sphere so we keep checking the old mode as well)

Signed-off-by: Chris Kulla <[email protected]>

Author: fpsunflower

src/doc/languagespec.tex

@@ -4671,10 +4671,17 @@ \subsection{Surface BSDF closures}
   % Optional string parameter to name this component. For use in AOVs / LPEs.
   % \apiend
 
+  \apiitem{energy_compensation}
+  \vspace{12pt}
+  Optional int parameter to select if energy compensation should be applied.
+  \apiend
+
 The Oren-Nayar reflection model is described in  M.\ Oren and S.\ K.\
 Nayar, ``Generalization of Lambert's Reflectance Model,'' Proceedings of
 SIGGRAPH 1994, pp.239-246 (July, 1994).
 
+The energy compensated model is described in the white paper: ``An energy-preserving Qualitative Oren-Nayar model'' by Jamie Portsmouth.
+
 \apiend
 
 

src/testrender/shading.cpp

@@ -101,6 +101,7 @@ struct MxOrenNayarDiffuseParams {
     float roughness;
     // optional
     ustring label;
+    int energy_compensation;
 };
 
 struct MxBurleyDiffuseParams {
@@ -316,6 +317,8 @@ register_closures(OSL::ShadingSystem* shadingsys)
             CLOSURE_COLOR_PARAM(MxOrenNayarDiffuseParams, albedo),
             CLOSURE_FLOAT_PARAM(MxOrenNayarDiffuseParams, roughness),
             CLOSURE_STRING_KEYPARAM(MxOrenNayarDiffuseParams, label, "label"),
+            CLOSURE_INT_KEYPARAM(MxOrenNayarDiffuseParams, energy_compensation,
+                                 "energy_compensation"),
             CLOSURE_FINISH_PARAM(MxOrenNayarDiffuseParams) } },
         { "burley_diffuse_bsdf",
           MX_BURLEY_DIFFUSE_ID,
@@ -459,23 +462,24 @@ template<int trans> struct Diffuse final : public BSDF, DiffuseParams {
 struct OrenNayar final : public BSDF, OrenNayarParams {
     OrenNayar(const OrenNayarParams& params) : BSDF(), OrenNayarParams(params)
     {
-        // precompute some constants
-        float s2 = sigma * sigma;
-        A        = 1 - 0.50f * s2 / (s2 + 0.33f);
-        B        = 0.45f * s2 / (s2 + 0.09f);
     }
     Sample eval(const Vec3& wo, const OSL::Vec3& wi) const override
     {
         float NL = N.dot(wi);
         float NV = N.dot(wo);
         if (NL > 0 && NV > 0) {
+            float LV = wo.dot(wi);
+            float s  = LV - NL * NV;
             // Simplified math from: "A tiny improvement of Oren-Nayar reflectance model"
             // by Yasuhiro Fujii
             // http://mimosa-pudica.net/improved-oren-nayar.html
-            // NOTE: This is using the math to match the original ON model, not the tweak
-            // proposed in the text which is a slightly different BRDF
-            float LV    = wo.dot(wi);
-            float s     = LV - NL * NV;
+            // NOTE: This is using the math to match the original qualitative ON model
+            // (QON in the paper above) and not the tweak proposed in the text which
+            // is a slightly different BRDF (FON in the paper above). This is done for
+            // backwards compatibility purposes only.
+            float s2    = sigma * sigma;
+            float A     = 1 - 0.50f * s2 / (s2 + 0.33f);
+            float B     = 0.45f * s2 / (s2 + 0.09f);
             float stinv = s > 0 ? s / std::max(NL, NV) : 0.0f;
             return { wi, Color3(A + B * stinv), NL * float(M_1_PI), 1.0f };
         }
@@ -489,9 +493,77 @@ struct OrenNayar final : public BSDF, OrenNayarParams {
         Sampling::sample_cosine_hemisphere(N, rx, ry, out_dir, pdf);
         return eval(wo, out_dir);
     }
+};
+
+struct EnergyCompensatedOrenNayar : public BSDF, MxOrenNayarDiffuseParams {
+    EnergyCompensatedOrenNayar(const MxOrenNayarDiffuseParams& params)
+        : BSDF(), MxOrenNayarDiffuseParams(params)
+    {
+    }
+    Sample eval(const Vec3& wo, const OSL::Vec3& wi) const override
+    {
+        float NL = N.dot(wi);
+        float NV = N.dot(wo);
+        if (NL > 0 && NV > 0) {
+            float LV = wo.dot(wi);
+            float s  = LV - NL * NV;
+            // Code below from Jamie Portsmouth's tech report on Energy conversion Oren-Nayar
+            // See slack thread for whitepaper:
+            // https://academysoftwarefdn.slack.com/files/U03SWQFPD08/F06S50CUKV1/oren_nayar.pdf
+
+            // TODO: rho should be the albedo which is a parameter of the closure in the Mx parameters
+            // This only matters for the color-saturation aspect of the BRDF which is rather subtle anyway
+            // and not always desireable for artists. Hardcoding to 1 leaves the coloring entirely up to the
+            // closure weight.
+
+            const Color3 rho  = albedo;
+            const float sigma = roughness;
+
+            float AF      = 1.0f / (1.0f + constant1_FON * sigma);
+            float stinv   = s > 0 ? s / std::max(NL, NV) : s;
+            float f_ss    = AF * (1.0 + sigma * stinv);  // single-scatt. BRDF
+            float EFo     = E_FON_analytic(NV);  // EFo at rho=1 (analytic)
+            float EFi     = E_FON_analytic(NL);  // EFi at rho=1 (analytic)
+            float avgEF   = AF * (1.0f + constant2_FON * sigma);  // avg. albedo
+            Color3 rho_ms = (rho * rho) * avgEF
+                            / (Color3(1.0f)
+                               - rho * std::max(0.0f, 1.0f - avgEF));
+            float f_ms = std::max(1e-7f, 1.0f - EFo)
+                         * std::max(1e-7f, 1.0f - EFi)
+                         / std::max(1e-7f, 1.0f - avgEF);  // multi-scatter lobe
+            return { wi, Color3(rho * f_ss + rho_ms * f_ms), NL * float(M_1_PI),
+                     1.0f };
+        }
+        return {};
+    }
+
+    Sample sample(const Vec3& wo, float rx, float ry,
+                  float /*rz*/) const override
+    {
+        Vec3 out_dir;
+        float pdf;
+        Sampling::sample_cosine_hemisphere(N, rx, ry, out_dir, pdf);
+        return eval(wo, out_dir);
+    }
 
 private:
-    float A, B;
+    static constexpr float constant1_FON = float(0.5 - 2.0 / (3.0 * M_PI));
+    static constexpr float constant2_FON = float(2.0 / 3.0
+                                                 - 28.0 / (15.0 * M_PI));
+
+    float E_FON_analytic(float mu) const
+    {
+        const float sigma = roughness;
+        float AF          = 1.0f
+                   / (1.0f
+                      + constant1_FON * sigma);  // Fujii model A coefficient
+        float BF = sigma * AF;                   // Fujii model B coefficient
+        float Si = sqrtf(std::max(0.0f, 1.0f - mu * mu));
+        float G  = Si * (OIIO::fast_acos(mu) - Si * mu)
+                  + 2.0 * ((Si / mu) * (1.0 - Si * Si * Si) - Si) / 3.0f;
+        float E = AF + (BF * float(M_1_PI)) * G;
+        return E;
+    }
 };
 
 struct Phong final : public BSDF, PhongParams {
@@ -1615,14 +1687,20 @@ process_bsdf_closure(const OSL::ShaderGlobals& sg, ShadingResult& result,
                 ok = result.bsdf.add_bsdf<Transparent>(cw);
                 break;
             case MX_OREN_NAYAR_DIFFUSE_ID: {
-                // translate MaterialX parameters into existing closure
                 const MxOrenNayarDiffuseParams* srcparams
                     = comp->as<MxOrenNayarDiffuseParams>();
-                OrenNayarParams params = {};
-                params.N               = srcparams->N;
-                params.sigma           = srcparams->roughness;
-                ok = result.bsdf.add_bsdf<OrenNayar>(cw * srcparams->albedo,
-                                                     params);
+                if (srcparams->energy_compensation) {
+                    // energy compensation handled by its own BSDF
+                    ok = result.bsdf.add_bsdf<EnergyCompensatedOrenNayar>(
+                        cw, *srcparams);
+                } else {
+                    // translate MaterialX parameters into existing closure
+                    OrenNayarParams params = {};
+                    params.N               = srcparams->N;
+                    params.sigma           = srcparams->roughness;
+                    ok = result.bsdf.add_bsdf<OrenNayar>(cw * srcparams->albedo,
+                                                         params);
+                }
                 break;
             }
             case MX_BURLEY_DIFFUSE_ID: {

testsuite/render-mx-furnace-oren-nayar/matte.osl

@@ -18,5 +18,5 @@ matte
             float UImin = 0, float UImax = 1 ]]
   )
 {
-    Ci = Kd * oren_nayar_diffuse_bsdf (N, Cs, roughness);
+    Ci = Kd * oren_nayar_diffuse_bsdf (N, Cs, roughness, "energy_compensation", N.y < 0);
 }