Implement MX conductor with multiple scattering

src/libbsdl/bsdl.cmake

@@ -21,7 +21,7 @@ function(ADD_BSDL_LIBRARY NAME)
 
     if (DEFINED bsdl_SPECTRAL_COLOR_SPACES)
         add_executable(jakobhanika_luts ${CMAKE_CURRENT_SOURCE_DIR}/${bsdl_SUBDIR}/src/jakobhanika_luts.cpp)
-        target_link_libraries(genluts PRIVATE Threads::Threads)
+        target_link_libraries(jakobhanika_luts PRIVATE Threads::Threads)
         foreach(CS ${bsdl_SPECTRAL_COLOR_SPACES})
             set(JACOBHANIKA_${CS} ${CMAKE_CURRENT_BINARY_DIR}/jakobhanika_${CS}.cpp)
             list(APPEND BSDL_LUTS_CPP ${JACOBHANIKA_${CS}})

src/libbsdl/include/BSDL/MTX/bsdf_conductor_decl.h

@@ -0,0 +1,76 @@
+// Copyright Contributors to the Open Shading Language project.
+// SPDX-License-Identifier: BSD-3-Clause
+// https://github.com/AcademySoftwareFoundation/OpenShadingLanguage
+
+
+#pragma once
+
+#include <BSDL/bsdf_decl.h>
+#include <BSDL/microfacet_tools_decl.h>
+
+BSDL_ENTER_NAMESPACE
+
+namespace mtx {
+
+struct ConductorFresnel {
+    BSDL_INLINE_METHOD ConductorFresnel() {}
+
+    BSDL_INLINE_METHOD
+    ConductorFresnel(Power IOR, Power extinction, float lambda_0);
+    BSDL_INLINE_METHOD Power eval(float cos_theta) const;
+    BSDL_INLINE_METHOD Power F0() const;
+
+    BSDL_INLINE_METHOD Power avg() const;
+
+private:
+    Power IOR, extinction;
+    float lambda_0;
+};
+
+template<typename BSDF_ROOT> struct ConductorLobe : public Lobe<BSDF_ROOT> {
+    using Base = Lobe<BSDF_ROOT>;
+    struct Data {
+        // microfacet params
+        Imath::V3f N, U;
+        float roughness_x;
+        float roughness_y;
+        // fresnel params
+        Imath::C3f IOR;
+        Imath::C3f extinction;
+        const char* distribution;
+        using lobe_type = ConductorLobe<BSDF_ROOT>;
+    };
+    template<typename D> static typename LobeRegistry<D>::Entry entry()
+    {
+        static_assert(
+            std::is_base_of<Data, D>::value);  // Make no other assumptions
+        using R = LobeRegistry<D>;
+        return { name(),
+                 { R::param(&D::N), R::param(&D::U), R::param(&D::roughness_x),
+                   R::param(&D::roughness_y), R::param(&D::IOR),
+                   R::param(&D::extinction), R::param(&D::distribution),
+                   R::close() } };
+    }
+
+    template<typename T>
+    BSDL_INLINE_METHOD ConductorLobe(T*, const BsdfGlobals& globals,
+                                     const Data& data);
+
+    static const char* name() { return "conductor_bsdf"; }
+
+    BSDL_INLINE_METHOD Power albedo_impl() const { return fresnel.avg(); }
+
+    BSDL_INLINE_METHOD Sample eval_impl(const Imath::V3f& wo,
+                                        const Imath::V3f& wi) const;
+    BSDL_INLINE_METHOD Sample sample_impl(const Imath::V3f& wo,
+                                          const Imath::V3f& rnd) const;
+
+private:
+    GGXDist dist;
+    ConductorFresnel fresnel;
+    float E_ms;
+};
+
+}  // namespace mtx
+
+BSDL_LEAVE_NAMESPACE

src/libbsdl/include/BSDL/MTX/bsdf_conductor_impl.h

@@ -0,0 +1,138 @@
+// Copyright Contributors to the Open Shading Language project.
+// SPDX-License-Identifier: BSD-3-Clause
+// https://github.com/AcademySoftwareFoundation/OpenShadingLanguage
+
+
+#pragma once
+
+#include <BSDL/MTX/bsdf_conductor_decl.h>
+
+BSDL_ENTER_NAMESPACE
+
+namespace mtx {
+
+BSDL_INLINE_METHOD
+ConductorFresnel::ConductorFresnel(Power IOR, Power extinction, float lambda_0)
+    : IOR(IOR), extinction(extinction), lambda_0(lambda_0)
+{
+}
+
+BSDL_INLINE_METHOD Power
+ConductorFresnel::eval(float cos_theta) const
+{
+    cos_theta = CLAMP(cos_theta, 0.0f, 1.0f);
+    const Power one(1, lambda_0);
+    const Power cosTheta2(cos_theta * cos_theta, lambda_0);
+    const Power sinTheta2 = one - cosTheta2;
+    const Power& n        = IOR;
+    const Power& k        = extinction;
+    const Power n2        = n * n;
+    const Power k2        = k * k;
+    const Power t0        = n2 - k2 - sinTheta2;
+    const Power a2plusb2  = sqrt(t0 * t0 + 4 * n2 * k2);
+    const Power t1        = a2plusb2 + cosTheta2;
+    const Power a         = sqrt(0.5f * (a2plusb2 + t0));
+    const Power t2        = (2.0f * cos_theta) * a;
+    const Power rs        = (t1 - t2) / (t1 + t2).clamped(FLOAT_MIN, BIG);
+
+    const Power t3 = cosTheta2 * a2plusb2 + sinTheta2 * sinTheta2;
+    const Power t4 = t2 * sinTheta2;
+    const Power rp = rs * (t3 - t4) / (t3 + t4).clamped(FLOAT_MIN, BIG);
+
+    return 0.5f * (rp + rs).clamped(0, 2);
+}
+
+BSDL_INLINE_METHOD Power
+ConductorFresnel::F0() const
+{
+    const Power one(1, lambda_0);
+    const Power& n       = IOR;
+    const Power& k       = extinction;
+    const Power n2       = n * n;
+    const Power k2       = k * k;
+    const Power t0       = n2 - k2;
+    const Power a2plusb2 = sqrt(t0 * t0 + 4 * n2 * k2);
+    const Power t1       = a2plusb2 + one;
+    const Power a        = sqrt(0.5f * (a2plusb2 + t0));
+    const Power t2       = 2.0f * a;
+    const Power rs       = (t1 - t2) / (t1 + t2).clamped(FLOAT_MIN, BIG);
+
+    return rs.clamped(0, 1);
+}
+
+BSDL_INLINE_METHOD Power
+ConductorFresnel::avg() const
+{
+    return Power(
+        [&](int i) {
+            // Very simple fit for cosine weighted average fresnel. Not very
+            // accurate but enough for albedo based sampling decisions.
+            constexpr float a = -0.32775145f, b = 0.18346033f, c = 0.61146583f,
+                            d = -0.07785134f;
+            const float x = IOR[i], y = extinction[i];
+            const float p = a + b * x + c * y + d * x * y;
+            return p / (1 + p);
+        },
+        lambda_0);
+}
+
+template<typename BSDF_ROOT>
+template<typename T>
+BSDL_INLINE_METHOD
+ConductorLobe<BSDF_ROOT>::ConductorLobe(T* lobe, const BsdfGlobals& globals,
+                                        const Data& data)
+    : Base(lobe, globals.visible_normal(data.N), data.U, 0.0f, globals.lambda_0,
+           false)
+{
+    Base::sample_filter = globals.get_sample_filter(Base::frame.Z, true);
+    // MaterialX expects the raw x/y roughness as input, but for albedo tables it
+    // is better to use the roughness/anisotropy parametrization so we can
+    // ignore roughness
+    const float rx    = CLAMP(data.roughness_x, EPSILON, 2.0f);
+    const float ry    = CLAMP(data.roughness_y, EPSILON, 2.0f);
+    const float ax    = std::max(rx, ry);
+    const float ay    = std::min(rx, ry);
+    const float b     = ay / ax;
+    const float aniso = (1 - b) / (1 + b);
+    // Also assume we square the roughness for linearity
+    const float roughness = globals.regularize_roughness(
+        sqrtf(ax / (1 + aniso)));
+    const float cosNO = Base::frame.Z.dot(globals.wo);
+    // Flip aniso if roughness_x < roughness_y
+    dist    = GGXDist(roughness, aniso, (rx < ry));
+    fresnel = ConductorFresnel(globals.wave(data.IOR),
+                               globals.wave(data.extinction), globals.lambda_0);
+    TabulatedEnergyCurve<spi::MiniMicrofacetGGX> curve(roughness, 0.0f);
+    E_ms = curve.Emiss_eval(cosNO);
+    Base::set_roughness(roughness);
+}
+
+template<typename BSDF_ROOT>
+BSDL_INLINE_METHOD Sample
+ConductorLobe<BSDF_ROOT>::eval_impl(const Imath::V3f& wo,
+                                    const Imath::V3f& wi) const
+{
+    Sample s    = eval_turquin_microms_reflection(dist, fresnel, E_ms, wo, wi);
+    s.roughness = Base::roughness();
+    return s;
+}
+
+template<typename BSDF_ROOT>
+BSDL_INLINE_METHOD Sample
+ConductorLobe<BSDF_ROOT>::sample_impl(const Imath::V3f& wo,
+                                      const Imath::V3f& rnd) const
+{
+    const float cosNO = wo.z;
+    if (cosNO <= 0)
+        return {};
+
+    // sample microfacet (half vector)
+    // generate outgoing direction
+    Imath::V3f wi = reflect(wo, dist.sample(wo, rnd.x, rnd.y));
+    // evaluate brdf on outgoing direction
+    return eval_impl(wo, wi);
+}
+
+}  // namespace mtx
+
+BSDL_LEAVE_NAMESPACE

src/libbsdl/include/BSDL/microfacet_tools_decl.h

@@ -14,11 +14,14 @@ BSDL_ENTER_NAMESPACE
 struct GGXDist {
     BSDL_INLINE_METHOD GGXDist() : ax(0), ay(0) {}
 
-    BSDL_INLINE_METHOD GGXDist(float rough, float aniso)
+    BSDL_INLINE_METHOD GGXDist(float rough, float aniso,
+                               bool flip_aniso = false)
         : ax(SQR(rough)), ay(ax)
     {
         assert(rough >= 0 && rough <= 1);
         assert(aniso >= 0 && aniso <= 1);
+        if (flip_aniso)
+            aniso = -aniso;
         constexpr float ALPHA_MIN = 1e-5f;
         ax                        = std::max(ax * (1 + aniso), ALPHA_MIN);
         ay                        = std::max(ay * (1 - aniso), ALPHA_MIN);
@@ -127,4 +130,11 @@ template<typename Fresnel> struct MicrofacetMS {
     float Eo_avg;
 };
 
+// Turquin style microfacet with multiple scattering
+template<typename Dist, typename Fresnel>
+BSDL_INLINE Sample
+eval_turquin_microms_reflection(const Dist& dist, const Fresnel& fresnel,
+                                float E_ms, const Imath::V3f& wo,
+                                const Imath::V3f& wi);
+
 BSDL_LEAVE_NAMESPACE

src/libbsdl/include/BSDL/microfacet_tools_impl.h

@@ -298,4 +298,35 @@ MicrofacetMS<Fresnel>::computeFmiss() const
     return Fmiss;
 }
 
+// Turquin style microfacet with multiple scattering
+template<typename Dist, typename Fresnel>
+BSDL_INLINE Sample
+eval_turquin_microms_reflection(const Dist& dist, const Fresnel& fresnel,
+                                float E_ms, const Imath::V3f& wo,
+                                const Imath::V3f& wi)
+{
+    const float cosNO = wo.z;
+    const float cosNI = wi.z;
+    if (cosNI <= 0 || cosNO <= 0)
+        return {};
+
+    // get half vector
+    Imath::V3f m    = (wo + wi).normalized();
+    float cosMO     = m.dot(wo);
+    const float D   = dist.D(m);
+    const float G1  = dist.G1(wo);
+    const float out = dist.G2_G1(wi, wo);
+    float s_pdf     = (G1 * D) / (4.0f * cosNO);
+    // fresnel term between outgoing direction and microfacet
+    const Power F = fresnel.eval(cosMO);
+    // From "Practical multiple scattering compensation for microfacet models" - Emmanuel Turquin
+    // equation 16. Should we use F0 for msf scaling? Doesn't make a big difference.
+    const Power F_ms = F;
+    const float msf  = E_ms / (1 - E_ms);
+    const Power one(1, 1);
+    const Power O = out * F * (one + F_ms * msf);
+
+    return { wi, O, s_pdf, -1 /* roughness set by caller */ };
+}
+
 BSDL_LEAVE_NAMESPACE

src/libbsdl/include/BSDL/spectrum_decl.h

@@ -307,6 +307,16 @@ struct Power {
     }
 };
 
+BSDL_INLINE Power
+sqrt(Power x)
+{
+    x = x.clamped(0.0f, std::numeric_limits<float>::max());
+    BSDL_UNROLL()
+    for (int i = 0; i != Power::N; ++i)
+        x.data[i] = sqrtf(x.data[i]);
+    return x;
+}
+
 BSDL_INLINE Power
 operator-(Power o)
 {