Backport of OpenPBR Surface (GLSL) to 1.38.10

libraries/pbrlib/genglsl/lib/mx39_microfacet.glsl

@@ -0,0 +1,30 @@
+#include "libraries/pbrlib/genglsl/lib/mx_microfacet.glsl"
+
+float mx39_pow6(float x)
+{
+    float x2 = mx_square(x);
+    return mx_square(x2) * x2;
+}
+
+// Generate a cosine-weighted sample on the unit hemisphere.
+vec3 mx39_cosine_sample_hemisphere(vec2 Xi)
+{
+    float phi = 2.0 * M_PI * Xi.x;
+    float cosTheta = sqrt(Xi.y);
+    float sinTheta = sqrt(1.0 - Xi.y);
+    return vec3(cos(phi) * sinTheta,
+                sin(phi) * sinTheta,
+                cosTheta);
+}
+
+// Construct an orthonormal basis from a unit vector.
+// https://graphics.pixar.com/library/OrthonormalB/paper.pdf
+mat3 mx39_orthonormal_basis(vec3 N)
+{
+    float sign = (N.z < 0.0) ? -1.0 : 1.0;
+    float a = -1.0 / (sign + N.z);
+    float b = N.x * N.y * a;
+    vec3 X = vec3(1.0 + sign * N.x * N.x * a, sign * b, -sign * N.x);
+    vec3 Y = vec3(b, sign + N.y * N.y * a, -N.y);
+    return mat3(X, Y, N);
+}

libraries/pbrlib/genglsl/lib/mx39_microfacet_diffuse.glsl

@@ -0,0 +1,88 @@
+#include "mx39_microfacet.glsl"
+#include "libraries/pbrlib/genglsl/lib/mx_microfacet_diffuse.glsl"
+
+const float FUJII_CONSTANT_1 = 0.5 - 2.0 / (3.0 * M_PI);
+const float FUJII_CONSTANT_2 = 2.0 / 3.0 - 28.0 / (15.0 * M_PI);
+
+// Qualitative Oren-Nayar diffuse with simplified math:
+// https://www1.cs.columbia.edu/CAVE/publications/pdfs/Oren_SIGGRAPH94.pdf
+float mx39_oren_nayar_diffuse(float NdotV, float NdotL, float LdotV, float roughness)
+{
+    float s = LdotV - NdotL * NdotV;
+    float stinv = (s > 0.0) ? s / max(NdotL, NdotV) : 0.0;
+
+    float sigma2 = mx_square(roughness);
+    float A = 1.0 - 0.5 * (sigma2 / (sigma2 + 0.33));
+    float B = 0.45 * sigma2 / (sigma2 + 0.09);
+
+    return A + B * stinv;
+}
+
+// Rational quadratic fit to Monte Carlo data for Oren-Nayar directional albedo.
+float mx39_oren_nayar_diffuse_dir_albedo_analytic(float NdotV, float roughness)
+{
+    vec2 r = vec2(1.0, 1.0) +
+             vec2(-0.4297, -0.6076) * roughness +
+             vec2(-0.7632, -0.4993) * NdotV * roughness +
+             vec2(1.4385, 2.0315) * mx_square(roughness);
+    return r.x / r.y;
+}
+
+float mx39_oren_nayar_diffuse_dir_albedo(float NdotV, float roughness)
+{
+    float dirAlbedo = mx39_oren_nayar_diffuse_dir_albedo_analytic(NdotV, roughness);
+    return clamp(dirAlbedo, 0.0, 1.0);
+}
+
+// Improved Oren-Nayar diffuse from Fujii:
+// https://mimosa-pudica.net/improved-oren-nayar.html
+float mx39_oren_nayar_fujii_diffuse_dir_albedo(float cosTheta, float roughness)
+{
+    float A = 1.0 / (1.0 + FUJII_CONSTANT_1 * roughness);
+    float B = roughness * A;
+    float Si = sqrt(max(0.0, 1.0 - mx_square(cosTheta)));
+    float G = Si * (acos(clamp(cosTheta, -1.0, 1.0)) - Si * cosTheta) +
+              2.0 * ((Si / cosTheta) * (1.0 - Si * Si * Si) - Si) / 3.0;
+    return A + (B * G * M_PI_INV);
+}
+
+float mx39_oren_nayar_fujii_diffuse_avg_albedo(float roughness)
+{
+    float A = 1.0 / (1.0 + FUJII_CONSTANT_1 * roughness);
+    return A * (1.0 + FUJII_CONSTANT_2 * roughness);
+}   
+
+// Energy-compensated Oren-Nayar diffuse from OpenPBR Surface:
+// https://academysoftwarefoundation.github.io/OpenPBR/
+vec3 mx39_oren_nayar_compensated_diffuse(float NdotV, float NdotL, float LdotV, float roughness, vec3 color)
+{
+    float s = LdotV - NdotL * NdotV;
+    float stinv = (s > 0.0) ? s / max(NdotL, NdotV) : s;
+
+    // Compute the single-scatter lobe.
+    float A = 1.0 / (1.0 + FUJII_CONSTANT_1 * roughness);
+    vec3 lobeSingleScatter = color * A * (1.0 + roughness * stinv);
+
+    // Compute the multi-scatter lobe.
+    float dirAlbedoV = mx39_oren_nayar_fujii_diffuse_dir_albedo(NdotV, roughness);
+    float dirAlbedoL = mx39_oren_nayar_fujii_diffuse_dir_albedo(NdotL, roughness);
+    float avgAlbedo = mx39_oren_nayar_fujii_diffuse_avg_albedo(roughness);
+    vec3 colorMultiScatter = mx_square(color) * avgAlbedo /
+                             (vec3(1.0) - color * max(0.0, 1.0 - avgAlbedo));
+    vec3 lobeMultiScatter = colorMultiScatter *
+                            max(M_FLOAT_EPS, 1.0 - dirAlbedoV) *
+                            max(M_FLOAT_EPS, 1.0 - dirAlbedoL) /
+                            max(M_FLOAT_EPS, 1.0 - avgAlbedo);
+
+    // Return the sum.
+    return lobeSingleScatter + lobeMultiScatter;
+}
+
+vec3 mx39_oren_nayar_compensated_diffuse_dir_albedo(float cosTheta, float roughness, vec3 color)
+{
+    float dirAlbedo = mx39_oren_nayar_fujii_diffuse_dir_albedo(cosTheta, roughness);
+    float avgAlbedo = mx39_oren_nayar_fujii_diffuse_avg_albedo(roughness);
+    vec3 colorMultiScatter = mx_square(color) * avgAlbedo /
+                             (vec3(1.0) - color * max(0.0, 1.0 - avgAlbedo));
+    return mix(colorMultiScatter, color, dirAlbedo);
+}

libraries/pbrlib/genglsl/lib/mx39_microfacet_sheen.glsl

@@ -0,0 +1,99 @@
+#include "mx39_microfacet.glsl"
+#include "libraries/pbrlib/genglsl/lib/mx_microfacet_sheen.glsl"
+
+// The following functions are adapted from https://github.com/tizian/ltc-sheen.
+// "Practical Multiple-Scattering Sheen Using Linearly Transformed Cosines", Zeltner et al.
+
+// Gaussian fit to directional albedo table.
+float mx39_zeltner_sheen_dir_albedo(float x, float y)
+{
+    float s = y*(0.0206607 + 1.58491*y)/(0.0379424 + y*(1.32227 + y));
+    float m = y*(-0.193854 + y*(-1.14885 + y*(1.7932 - 0.95943*y*y)))/(0.046391 + y);
+    float o = y*(0.000654023 + (-0.0207818 + 0.119681*y)*y)/(1.26264 + y*(-1.92021 + y));
+    return exp(-0.5*mx_square((x - m)/s))/(s*sqrt(2.0*M_PI)) + o;
+}
+
+// Rational fits to LTC matrix coefficients.
+float mx39_zeltner_sheen_ltc_aInv(float x, float y)
+{
+    return (2.58126*x + 0.813703*y)*y/(1.0 + 0.310327*x*x + 2.60994*x*y);
+}
+
+float mx39_zeltner_sheen_ltc_bInv(float x, float y)
+{
+    return sqrt(1.0 - x)*(y - 1.0)*y*y*y/(0.0000254053 + 1.71228*x - 1.71506*x*y + 1.34174*y*y);
+}
+
+// V and N are assumed to be unit vectors.
+mat3 mx39_orthonormal_basis_ltc(vec3 V, vec3 N, float NdotV)
+{
+    // Generate a tangent vector in the plane of V and N.
+    // This required to correctly orient the LTC lobe.
+    vec3 X = V - N*NdotV;
+    float lenSqr = dot(X, X);
+    if (lenSqr > 0.0)
+    {
+        X *= inversesqrt(lenSqr);
+        vec3 Y = cross(N, X);
+        return mat3(X, Y, N);
+    }
+
+    // If lenSqr == 0, then V == N, so any orthonormal basis will do.
+    return mx39_orthonormal_basis(N);
+}
+
+// Multiplication by directional albedo is handled by the calling function.
+float mx39_zeltner_sheen_brdf(vec3 L, vec3 V, vec3 N, float NdotV, float roughness)
+{
+    mat3 toLTC = transpose(mx39_orthonormal_basis_ltc(V, N, NdotV));
+    vec3 w = toLTC * L;
+
+    float aInv = mx39_zeltner_sheen_ltc_aInv(NdotV, roughness);
+    float bInv = mx39_zeltner_sheen_ltc_bInv(NdotV, roughness);
+
+    // Transform w to original configuration (clamped cosine).
+    //                 |aInv    0 bInv|
+    // wo = M^-1 . w = |   0 aInv    0| . w
+    //                 |   0    0    1|
+    vec3 wo = vec3(aInv*w.x + bInv*w.z, aInv * w.y, w.z);
+    float lenSqr = dot(wo, wo);
+
+    // D(w) = Do(M^-1.w / ||M^-1.w||) . |M^-1| / ||M^-1.w||^3
+    //      = Do(M^-1.w) . |M^-1| / ||M^-1.w||^4
+    //      = Do(wo) . |M^-1| / dot(wo, wo)^2
+    //      = Do(wo) . aInv^2 / dot(wo, wo)^2
+    //      = Do(wo) . (aInv / dot(wo, wo))^2
+    return max(wo.z, 0.0) * M_PI_INV * mx_square(aInv / lenSqr);
+}
+
+vec3 mx39_zeltner_sheen_importance_sample(vec2 Xi, vec3 V, vec3 N, float roughness, out float pdf)
+{
+    float NdotV = clamp(dot(N, V), 0.0, 1.0);
+    roughness = clamp(roughness, 0.01, 1.0); // Clamp to range of original impl.
+
+    vec3 wo = mx39_cosine_sample_hemisphere(Xi);
+
+    float aInv = mx39_zeltner_sheen_ltc_aInv(NdotV, roughness);
+    float bInv = mx39_zeltner_sheen_ltc_bInv(NdotV, roughness);
+
+    // Transform wo from original configuration (clamped cosine).
+    //              |1/aInv      0 -bInv/aInv|
+    // w = M . wo = |     0 1/aInv          0| . wo
+    //              |     0      0          1|    
+    vec3 w = vec3(wo.x/aInv - wo.z*bInv/aInv, wo.y / aInv, wo.z);
+
+    float lenSqr = dot(w, w);
+    w *= inversesqrt(lenSqr);
+
+    // D(w) = Do(wo) . ||M.wo||^3 / |M|
+    //      = Do(wo / ||M.wo||) . ||M.wo||^4 / |M| 
+    //      = Do(w) . ||M.wo||^4 / |M| (possible because M doesn't change z component)
+    //      = Do(w) . dot(w, w)^2 * aInv^2
+    //      = Do(w) . (aInv * dot(w, w))^2
+    pdf = max(w.z, 0.0) * M_PI_INV * mx_square(aInv * lenSqr);
+
+    mat3 fromLTC = mx39_orthonormal_basis_ltc(V, N, NdotV);
+    w = fromLTC * w;
+
+    return w;
+}