Proposed improvements to energy conservation

src/testrender/shading.cpp

@@ -885,10 +885,8 @@ struct MxMicrofacet final : public BSDF, MxMicrofacetParams {
         // if transmission is enabled, punt on
         if (EnableTransmissionLobe)
             return Color3(1.0f);
-        // FIXME: this heuristic is not particularly good, and looses energy
-        // compared to the reference solution
 
-        return MxMicrofacetParams::evalR(
+        return MxMicrofacetParams::dirAlbedoR(
             get_fresnel_angle(MxMicrofacetParams::N.dot(wo)));
     }
 

src/testrender/shading.h

@@ -135,6 +135,35 @@ struct MxDielectricParams : public MxMicrofacetBaseParams {
     {
         return transmission_tint * (1.0f - fresnel_dielectric(cos_theta, ior));
     }
+
+    OSL_HOSTDEVICE Color3 dirAlbedoR(float cos_theta) const
+    {
+        float iorRatio = (ior - 1.0f) / (ior + 1.0f);
+        Color3 f0(iorRatio * iorRatio);
+        Color3 f90(1.0f);
+
+        // Rational quadratic fit for GGX directional albedo
+        // https://github.com/AcademySoftwareFoundation/MaterialX/blob/main/libraries/pbrlib/genglsl/lib/mx_microfacet_specular.glsl
+        float x  = OIIO::clamp(cos_theta, 0.0f, 1.0f);
+        float y  = sqrtf(roughness_x * roughness_y);  // average alpha
+        float x2 = x * x;
+        float y2 = y * y;
+        Vec2 num = Vec2(0.1003f, 0.9345f) + Vec2(-0.6303f, -2.323f) * x
+                   + Vec2(9.748f, 2.229f) * y + Vec2(-2.038f, -3.748f) * x * y
+                   + Vec2(29.34f, 1.424f) * x2 + Vec2(-8.245f, -0.7684f) * y2
+                   + Vec2(-26.44f, 1.436f) * x2 * y
+                   + Vec2(19.99f, 0.2913f) * x * y2
+                   + Vec2(-5.448f, 0.6286f) * x2 * y2;
+        Vec2 den = Vec2(1.0f, 1.0f) + Vec2(-1.765f, 0.2281f) * x
+                   + Vec2(8.263f, 15.94f) * y + Vec2(11.53f, -55.83f) * x * y
+                   + Vec2(28.96f, 13.08f) * x2 + Vec2(-7.507f, 41.26f) * y2
+                   + Vec2(-36.11f, 54.9f) * x2 * y
+                   + Vec2(15.86f, 300.2f) * x * y2
+                   + Vec2(33.37f, -285.1f) * x2 * y2;
+        float a = OIIO::clamp(num.x / den.x, 0.0f, 1.0f);
+        float b = OIIO::clamp(num.y / den.y, 0.0f, 1.0f);
+        return reflection_tint * (f0 * a + f90 * b);
+    }
 };
 
 struct MxConductorParams : public MxMicrofacetBaseParams {
@@ -151,6 +180,13 @@ struct MxConductorParams : public MxMicrofacetBaseParams {
 
     OSL_HOSTDEVICE Color3 evalT(float cos_theta) const { return Color3(0.0f); }
 
+    OSL_HOSTDEVICE Color3 dirAlbedoR(float cos_theta) const
+    {
+        // TODO: Integrate the MaterialX fit for GGX directional albedo, which
+        //       may improve multiscatter compensation for conductors.
+        return evalR(cos_theta);
+    }
+
     // Avoid function was declared but never referenced
     // float get_ior() const
     // {
@@ -180,6 +216,31 @@ struct MxGeneralizedSchlickParams : public MxMicrofacetBaseParams {
                * (Color3(1.0f)
                   - fresnel_generalized_schlick(cos_theta, f0, f90, exponent));
     }
+
+    OSL_HOSTDEVICE Color3 dirAlbedoR(float cos_theta) const
+    {
+        // Rational quadratic fit for GGX directional albedo
+        // https://github.com/AcademySoftwareFoundation/MaterialX/blob/main/libraries/pbrlib/genglsl/lib/mx_microfacet_specular.glsl
+        float x  = OIIO::clamp(cos_theta, 0.0f, 1.0f);
+        float y  = sqrtf(roughness_x * roughness_y);  // average alpha
+        float x2 = x * x;
+        float y2 = y * y;
+        Vec2 num = Vec2(0.1003f, 0.9345f) + Vec2(-0.6303f, -2.323f) * x
+                   + Vec2(9.748f, 2.229f) * y + Vec2(-2.038f, -3.748f) * x * y
+                   + Vec2(29.34f, 1.424f) * x2 + Vec2(-8.245f, -0.7684f) * y2
+                   + Vec2(-26.44f, 1.436f) * x2 * y
+                   + Vec2(19.99f, 0.2913f) * x * y2
+                   + Vec2(-5.448f, 0.6286f) * x2 * y2;
+        Vec2 den = Vec2(1.0f, 1.0f) + Vec2(-1.765f, 0.2281f) * x
+                   + Vec2(8.263f, 15.94f) * y + Vec2(11.53f, -55.83f) * x * y
+                   + Vec2(28.96f, 13.08f) * x2 + Vec2(-7.507f, 41.26f) * y2
+                   + Vec2(-36.11f, 54.9f) * x2 * y
+                   + Vec2(15.86f, 300.2f) * x * y2
+                   + Vec2(33.37f, -285.1f) * x2 * y2;
+        float a = OIIO::clamp(num.x / den.x, 0.0f, 1.0f);
+        float b = OIIO::clamp(num.y / den.y, 0.0f, 1.0f);
+        return reflection_tint * (f0 * a + f90 * b);
+    }
 };
 
 struct MxTranslucentParams {