Improvements to GGX importance sampling

- Implement the paper "Sampling Visible GGX Normals with Spherical Caps" by Jonathan Dupuy and Anis Benyoub, which improves the performance and spatial continuity of VNDF sampling.
- Switch to VNDF sampling for FIS environment lights, improving the convergence of this lighting path for lower sample counts.
- Additional optimizations for FIS environment lights, leveraging the improved term cancellation in VNDF sampling.

Test results:
- GLSL render performance is improved for all tested materials, e.g. an increase from 205 fps to 214 fps for standard_surface_default.mtlx with 16 environment samples on an NVIDIA RTX A6000.
- Convergence of FIS environment lights is improved for all tested materials, with the maximum visual error between 16 and 4096 environment samples reduced from 0.11 to 0.07 for a rough gold material.

libraries/pbrlib/genglsl/lib/mx_environment_fis.glsl

@@ -23,6 +23,7 @@ vec3 mx_environment_radiance(vec3 N, vec3 V, vec3 X, vec2 alpha, int distributio
     // Compute derived properties.
     float NdotV = clamp(V.z, M_FLOAT_EPS, 1.0);
     float avgAlpha = mx_average_alpha(alpha);
+    float G1V = mx_ggx_smith_G1(NdotV, avgAlpha);
 
     // Integrate outgoing radiance using filtered importance sampling.
     // http://cgg.mff.cuni.cz/~jaroslav/papers/2008-egsr-fis/2008-egsr-fis-final-embedded.pdf
@@ -33,18 +34,16 @@ vec3 mx_environment_radiance(vec3 N, vec3 V, vec3 X, vec2 alpha, int distributio
         vec2 Xi = mx_spherical_fibonacci(i, envRadianceSamples);
 
         // Compute the half vector and incoming light direction.
-        vec3 H = mx_ggx_importance_sample_NDF(Xi, alpha);
+        vec3 H = mx_ggx_importance_sample_VNDF(Xi, V, alpha);
         vec3 L = fd.refraction ? mx_refraction_solid_sphere(-V, H, fd.ior.x) : -reflect(V, H);
 
         // Compute dot products for this sample.
-        float NdotH = clamp(H.z, M_FLOAT_EPS, 1.0);
         float NdotL = clamp(L.z, M_FLOAT_EPS, 1.0);
         float VdotH = clamp(dot(V, H), M_FLOAT_EPS, 1.0);
-        float LdotH = VdotH;
 
         // Sample the environment light from the given direction.
         vec3 Lw = tangentToWorld * L;
-        float pdf = mx_ggx_PDF(H, LdotH, alpha);
+        float pdf = mx_ggx_NDF(H, alpha) * G1V / (4.0 * NdotV);
         float lod = mx_latlong_compute_lod(Lw, pdf, float($envRadianceMips - 1), envRadianceSamples);
         vec3 sampleColor = mx_latlong_map_lookup(Lw, $envMatrix, lod, $envRadiance);
 
@@ -61,13 +60,15 @@ vec3 mx_environment_radiance(vec3 N, vec3 V, vec3 X, vec2 alpha, int distributio
         // From https://cdn2.unrealengine.com/Resources/files/2013SiggraphPresentationsNotes-26915738.pdf
         //   incidentLight = sampleColor * NdotL
         //   microfacetSpecular = D * F * G / (4 * NdotL * NdotV)
-        //   pdf = D * NdotH / (4 * VdotH)
+        //   pdf = D * G1V / (4 * NdotV);
         //   radiance = incidentLight * microfacetSpecular / pdf
-        radiance += sampleColor * FG * VdotH / (NdotV * NdotH);
+        radiance += sampleColor * FG;
     }
 
-    // Normalize and return the final radiance.
-    radiance /= float(envRadianceSamples);
+    // Apply the global component of the geometric term and normalize.
+    radiance /= G1V * float(envRadianceSamples);
+
+    // Return the final radiance.
     return radiance;
 }
 

libraries/pbrlib/genglsl/lib/mx_microfacet_specular.glsl

@@ -61,48 +61,22 @@ float mx_ggx_NDF(vec3 H, vec2 alpha)
     return 1.0 / (M_PI * alpha.x * alpha.y * mx_square(denom));
 }
 
-// https://media.disneyanimation.com/uploads/production/publication_asset/48/asset/s2012_pbs_disney_brdf_notes_v3.pdf
-// Appendix B.1 Equation 3
-float mx_ggx_PDF(vec3 H, float LdotH, vec2 alpha)
-{
-    float NdotH = H.z;
-    return mx_ggx_NDF(H, alpha) * NdotH / (4.0 * LdotH);
-}
-
-// https://media.disneyanimation.com/uploads/production/publication_asset/48/asset/s2012_pbs_disney_brdf_notes_v3.pdf
-// Appendix B.2 Equation 15
-vec3 mx_ggx_importance_sample_NDF(vec2 Xi, vec2 alpha)
-{
-    float phi = 2.0 * M_PI * Xi.x;
-    float tanTheta = sqrt(Xi.y / (1.0 - Xi.y));
-    vec3 H = vec3(tanTheta * alpha.x * cos(phi),
-                  tanTheta * alpha.y * sin(phi),
-                  1.0);
-    return normalize(H);
-}
-
-// http://jcgt.org/published/0007/04/01/paper.pdf
-// Appendix A Listing 1
+// https://ggx-research.github.io/publication/2023/06/09/publication-ggx.html
 vec3 mx_ggx_importance_sample_VNDF(vec2 Xi, vec3 V, vec2 alpha)
 {
     // Transform the view direction to the hemisphere configuration.
     V = normalize(vec3(V.xy * alpha, V.z));
 
-    // Construct an orthonormal basis from the view direction.
-    float len = length(V.xy);
-    vec3 T1 = (len > 0.0) ? vec3(-V.y, V.x, 0.0) / len : vec3(1.0, 0.0, 0.0);
-    vec3 T2 = cross(V, T1);
-
-    // Parameterization of the projected area.
-    float r = sqrt(Xi.y);
+    // Sample a spherical cap in (-V.z, 1].
     float phi = 2.0 * M_PI * Xi.x;
-    float t1 = r * cos(phi);
-    float t2 = r * sin(phi);
-    float s = 0.5 * (1.0 + V.z);
-    t2 = (1.0 - s) * sqrt(1.0 - mx_square(t1)) + s * t2;
-
-    // Reprojection onto hemisphere.
-    vec3 H = t1 * T1 + t2 * T2 + sqrt(max(0.0, 1.0 - mx_square(t1) - mx_square(t2))) * V;
+    float z = (1.0 - Xi.y) * (1.0 + V.z) - V.z;
+    float sinTheta = sqrt(clamp(1.0 - z * z, 0.0, 1.0));
+    float x = sinTheta * cos(phi);
+    float y = sinTheta * sin(phi);
+    vec3 c = vec3(x, y, z);
+
+    // Compute the microfacet normal.
+    vec3 H = c + V;
 
     // Transform the microfacet normal back to the ellipsoid configuration.
     H = normalize(vec3(H.xy * alpha, max(H.z, 0.0)));