Implement GGX multiscatter energy compensation using Cycles approach

Author: mrdoob

src/shader/bsdf/bsdf_functions.glsl.js

@@ -72,11 +72,12 @@ export const bsdf_functions = /* glsl */`
 		// Single-scatter term (standard Cook-Torrance microfacet BRDF)
 		vec3 singleScatter = wi.z * F * G * D / ( 4.0 * abs( wi.z * wo.z ) );
 
-		// Multi-scatter energy compensation (Kulla-Conty 2017)
-		// This accounts for energy lost due to multiple bounces within the microfacet structure
-		// The multiscatter term is already divided by PI and accounts for cosine weighting
+		// Multiscatter energy compensation
+		// Adds back energy lost to multiple bounces within the microfacet structure
+		// Returns a diffuse-like lobe scaled by the missing energy and average Fresnel
 		vec3 multiScatter = ggxMultiScatterCompensation( wo, wi, roughness, f0Color ) * wi.z;
 
+		// Total specular reflection = single-scatter + multiscatter
 		color = singleScatter + multiScatter;
 		return ggxPdf / ( 4.0 * dot( wo, wh ) );
 
@@ -207,10 +208,12 @@ export const bsdf_functions = /* glsl */`
 		// Single-scatter clearcoat term
 		float fClearcoatSingle = F * D * G / ( 4.0 * abs( wi.z * wo.z ) );
 
-		// Multi-scatter compensation for clearcoat layer
+		// Multiscatter energy compensation for clearcoat layer
+		// Clearcoat is a dielectric layer (IOR 1.5), so we use its F0 for compensation
 		vec3 f0ColorClearcoat = vec3( f0 );
 		vec3 clearcoatMultiScatter = ggxMultiScatterCompensation( wo, wi, roughness, f0ColorClearcoat );
 
+		// Total clearcoat reflection = single-scatter + multiscatter
 		float fClearcoat = fClearcoatSingle + clearcoatMultiScatter.r;
 		color = color * ( 1.0 - surf.clearcoat * F ) + fClearcoat * surf.clearcoat * wi.z;
 

src/shader/bsdf/multiscatter_functions.glsl.js

@@ -1,26 +1,53 @@
 export const multiscatter_functions = /* glsl */`
 
-// Analytical multiscatter energy compensation for GGX BRDF
-// Compensates for energy loss due to multiple bounces within the microfacet structure
-// Based on observations that rough surfaces at grazing angles lose the most energy
+// GGX Multiscatter Energy Compensation
+// Implementation based on Blender Cycles' approach
+//
+// References:
+// - "Multiple-Scattering Microfacet BSDFs with the Smith Model" (Heitz et al. 2016)
+// - Blender Cycles: intern/cycles/kernel/closure/bsdf_microfacet_multi.h
+//
+// Single-scatter GGX loses energy due to rays bouncing multiple times within
+// the microfacet structure before escaping. This compensation adds back the
+// missing energy as a diffuse-like multiscatter lobe.
+//
+// The approach uses a fitted albedo approximation to estimate how much energy
+// single-scatter captures, then adds the remainder back. This is simpler and
+// faster than full random-walk multiscatter simulation while providing good
+// energy conservation for path tracers.
+
+// Directional albedo approximation for single-scatter GGX
+// Returns the fraction of energy captured by single-scatter as a function of roughness
+// Fitted curve from Blender Cycles based on precomputed ground truth data
+float ggxAlbedo( float roughness ) {
+	float r2 = roughness * roughness;
+	return 0.806495 * exp( -1.98712 * r2 ) + 0.199531;
+}
+
+// GGX multiscatter energy compensation term
+// wo: outgoing direction (view direction)
+// wi: incident direction (light direction)
+// roughness: surface roughness [0, 1]
+// F0: Fresnel reflectance at normal incidence
+// Returns: Additional BRDF contribution to compensate for multiscatter energy loss
 vec3 ggxMultiScatterCompensation( vec3 wo, vec3 wi, float roughness, vec3 F0 ) {
-	float NdotV = abs( wo.z );
-	float NdotL = abs( wi.z );
+	// Estimate the fraction of energy captured by single-scatter GGX
+	float singleScatterAlbedo = ggxAlbedo( roughness );
 
-	// Energy compensation increases with roughness
-	// At roughness=0, no compensation needed (perfect mirror)
-	// At roughness=1, significant compensation needed (very rough)
-	float a = roughness * roughness;
-	float energyFactor = a * sqrt( a );  // Scales as roughness^1.5
+	// The missing energy that needs compensation
+	float missingEnergy = 1.0 - singleScatterAlbedo;
 
-	// Angular dependence - more energy lost at grazing angles
-	float angularLoss = ( 1.0 - NdotV * 0.9 ) * ( 1.0 - NdotL * 0.9 );
+	// Average Fresnel reflectance over all directions (spherical albedo)
+	// Approximation: F_avg ≈ F0 + (1 - F0) / 21
+	vec3 Favg = F0 + ( 1.0 - F0 ) / 21.0;
 
-	// Combined energy compensation
-	vec3 compensation = F0 * energyFactor * angularLoss;
+	// Multiscatter contribution: diffuse-like lobe scaled by average Fresnel
+	// This represents energy that bounced multiple times before escaping
+	vec3 Fms = Favg * missingEnergy;
 
-	// Conservative global scale to avoid over-brightening
-	return compensation * 0.25;
+	// Return as a Lambertian BRDF (energy / π)
+	// The π accounts for the hemispherical integral in the rendering equation
+	return Fms / PI;
 }