Improve Fresnel weighting in GGX energy compensation (AcademySoftwareFoundation#2958)

This changelist improves the multiple-scattering energy compensation of the GGX microfacet model in hardware shading languages, tinting the recovered energy by the average Fresnel reflectance over the hemisphere rather than a single-direction Fresnel value.  The following specific changes are included:

- Update `mx_ggx_energy_compensation` to take a `FresnelData` argument and compute the cosine-weighted average of the Fresnel reflectance over the hemisphere internally, replacing the prior overloads that accepted a precomputed single-angle value.
- Update the conductor, dielectric, generalized Schlick, and Chiang hair closures to pass `FresnelData` directly, removing the redundant Fresnel evaluations that existed only to feed energy compensation.

As measured by our render test suite in GitHub CI, this change consistently improves the visual parity between MSL and OSL renders, with a few typical examples below (per-pixel RMS error, before -> after):

- `standard_surface_brick_procedural`: 0.00925 -> 0.00916
- `standard_surface_marble_procedural`: 0.00571 -> 0.00558
- `open_pbr_carpaint`: 0.00734 -> 0.00720
- `open_pbr_default`: 0.00627 -> 0.00619

Author: jstone-lucasfilm

libraries/pbrlib/genglsl/lib/mx_microfacet_specular.glsl

@@ -173,19 +173,6 @@ float mx_ggx_dir_albedo(float NdotV, float alpha, float F0, float F90)
     return mx_ggx_dir_albedo(NdotV, alpha, vec3(F0), vec3(F90)).x;
 }
 
-// https://blog.selfshadow.com/publications/turquin/ms_comp_final.pdf
-// Equations 14 and 16
-vec3 mx_ggx_energy_compensation(float NdotV, float alpha, vec3 Fss)
-{
-    float Ess = mx_ggx_dir_albedo(NdotV, alpha, 1.0, 1.0);
-    return 1.0 + Fss * (1.0 - Ess) / Ess;
-}
-
-float mx_ggx_energy_compensation(float NdotV, float alpha, float Fss)
-{
-    return mx_ggx_energy_compensation(NdotV, alpha, vec3(Fss)).x;
-}
-
 // Compute the average of an anisotropic alpha pair.
 float mx_average_alpha(vec2 alpha)
 {
@@ -511,6 +498,28 @@ vec3 mx_ggx_dir_albedo(float NdotV, float alpha, FresnelData fd)
     }
 }
 
+// Compute the cosine-weighted average of the Fresnel reflectance over the hemisphere.
+// https://blog.selfshadow.com/publications/s2017-shading-course/imageworks/s2017_pbs_imageworks_slides_v2.pdf
+vec3 mx_fresnel_average(FresnelData fd)
+{
+    vec3 F0 = mx_compute_fresnel(1.0, fd);
+    vec3 F90 = (fd.model == FRESNEL_MODEL_SCHLICK && !fd.airy) ? fd.F90 : vec3(1.0);
+
+    // The constant 1/21 is exact for a Schlick term with an exponent of 5, while for
+    // a generalized Schlick exponent n it would be 2 / ((n + 1) * (n + 2)).
+    return F0 + (F90 - F0) * (1.0 / 21.0);
+}
+
+// Multiple-scattering energy compensation for the GGX microfacet model.
+// https://blog.selfshadow.com/publications/turquin/ms_comp_final.pdf
+// Equations 14 and 16
+vec3 mx_ggx_energy_compensation(float NdotV, float alpha, FresnelData fd)
+{
+    vec3 Fss = mx_fresnel_average(fd);
+    float Ess = mx_ggx_dir_albedo(NdotV, alpha, 1.0, 1.0);
+    return 1.0 + Fss * (1.0 - Ess) / Ess;
+}
+
 // Compute the refraction of a ray through a solid sphere.
 vec3 mx_refraction_solid_sphere(vec3 R, vec3 N, float ior)
 {

libraries/pbrlib/genglsl/mx_chiang_hair_bsdf.glsl

@@ -272,15 +272,14 @@ void mx_chiang_hair_bsdf(ClosureData closureData, vec3 tint_R, vec3 tint_TT, vec
 
         float NdotV = clamp(dot(N, V), M_FLOAT_EPS, 1.0);
         FresnelData fd = mx_init_fresnel_dielectric(ior, 0.0, 1.0);
-        vec3 F = mx_compute_fresnel(NdotV, fd);
 
         vec2 roughness = (roughness_R + roughness_TT + roughness_TRT) / vec2(3.0);  // ?
         vec2 safeAlpha = clamp(roughness, M_FLOAT_EPS, 1.0);
         float avgAlpha = mx_average_alpha(safeAlpha);
 
         // Use GGX to match the behavior of mx_environment_radiance.
         float F0 = mx_ior_to_f0(ior);
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, F0, 1.0) * comp;
 
         vec3 Li = mx_environment_radiance(N, V, X, safeAlpha, 0, fd);

libraries/pbrlib/genglsl/mx_conductor_bsdf.glsl

@@ -37,15 +37,14 @@ void mx_conductor_bsdf(ClosureData closureData, float weight, vec3 ior_n, vec3 i
         float D = mx_ggx_NDF(Ht, safeAlpha);
         float G = mx_ggx_smith_G2(NdotL, NdotV, avgAlpha);
 
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
 
         // Note: NdotL is cancelled out
         bsdf.response = D * F * G * comp * closureData.occlusion * weight / (4.0 * NdotV);
     }
     else if (closureData.closureType == CLOSURE_TYPE_INDIRECT)
     {
-        vec3 F = mx_compute_fresnel(NdotV, fd);
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 Li = mx_environment_radiance(N, V, X, safeAlpha, distribution, fd);
         bsdf.response = Li * comp * weight;
     }

libraries/pbrlib/genglsl/mx_dielectric_bsdf.glsl

@@ -44,17 +44,15 @@ void mx_dielectric_bsdf(ClosureData closureData, float weight, vec3 tint, float
         float D = mx_ggx_NDF(Ht, safeAlpha);
         float G = mx_ggx_smith_G2(NdotL, NdotV, avgAlpha);
 
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, F0, 1.0) * comp;
         bsdf.throughput = 1.0 - dirAlbedo * weight;
 
         bsdf.response = D * F * G * comp * safeTint * closureData.occlusion * weight / (4.0 * NdotV);
     }
     else if (closureData.closureType == CLOSURE_TYPE_TRANSMISSION)
     {
-        vec3 F = mx_compute_fresnel(NdotV, fd);
-
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, F0, 1.0) * comp;
         bsdf.throughput = 1.0 - dirAlbedo * weight;
 
@@ -65,9 +63,7 @@ void mx_dielectric_bsdf(ClosureData closureData, float weight, vec3 tint, float
     }
     else if (closureData.closureType == CLOSURE_TYPE_INDIRECT)
     {
-        vec3 F = mx_compute_fresnel(NdotV, fd);
-
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, F0, 1.0) * comp;
         bsdf.throughput = 1.0 - dirAlbedo * weight;
 

libraries/pbrlib/genglsl/mx_generalized_schlick_bsdf.glsl

@@ -45,7 +45,7 @@ void mx_generalized_schlick_bsdf(ClosureData closureData, float weight, vec3 col
         float D = mx_ggx_NDF(Ht, safeAlpha);
         float G = mx_ggx_smith_G2(NdotL, NdotV, avgAlpha);
 
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, safeColor0, safeColor90) * comp;
         float avgDirAlbedo = dot(dirAlbedo, vec3(1.0 / 3.0));
         bsdf.throughput = vec3(1.0 - avgDirAlbedo * weight);
@@ -55,9 +55,7 @@ void mx_generalized_schlick_bsdf(ClosureData closureData, float weight, vec3 col
     }
     else if (closureData.closureType == CLOSURE_TYPE_TRANSMISSION)
     {
-        vec3 F = mx_compute_fresnel(NdotV, fd);
-
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, safeColor0, safeColor90) * comp;
         float avgDirAlbedo = dot(dirAlbedo, vec3(1.0 / 3.0));
         bsdf.throughput = vec3(1.0 - avgDirAlbedo * weight);
@@ -71,9 +69,7 @@ void mx_generalized_schlick_bsdf(ClosureData closureData, float weight, vec3 col
     }
     else if (closureData.closureType == CLOSURE_TYPE_INDIRECT)
     {
-        vec3 F = mx_compute_fresnel(NdotV, fd);
-
-        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, F);
+        vec3 comp = mx_ggx_energy_compensation(NdotV, avgAlpha, fd);
         vec3 dirAlbedo = mx_ggx_dir_albedo(NdotV, avgAlpha, safeColor0, safeColor90) * comp;
         float avgDirAlbedo = dot(dirAlbedo, vec3(1.0 / 3.0));
         bsdf.throughput = vec3(1.0 - avgDirAlbedo * weight);