Energy Compensation Pt. 2
As we saw in the previous part, the energy compensation method can be very efficient at restituting the missing energy from multiple scattering events.
For real time purposes though, 3 texture taps per pixel and per light can be quite expensive, especially if the concerned light source is not covering a large solid angle, in which case the effect of multiple scattering is quite lost and isn't really worth the cost.
Todo
Image showing a simple point/directional with and without MS, because I'm not even sure of what I'm writing is actually true!
This is why we will focus on light sources that cover an extended area, namely:
 The far field environment, encoded as Spherical Harmonics
 Area lights, encoded as linearly transformed cosines ^{1}
In this section, we will use the fact that the multiple scattering part of the full BRDF only manifests itself in the case of very rough surfaces which means that almost the entire hemisphere of directions above the surface normal is involved in the computation of multiplyscattered radiance and we can thus easily make the important assumption that the overall response of the MSBRDF to a large light source must be low frequency.
Preliminary Remark
Before tackling the actual lighting problems, just a quick remark on the irradiance tables discussed in part 1: they might seem like a costly additional texture to add to the rendering pipeline but it's important to realize that these textures are probably already present in your PBR renderer!
Indeed, Brian Karis introduced a now wellknown splitsum approximation to the lighting equation for imagebased lighting in his 2013 Siggraph presentation^{2}:
The first term (teal) of the right expression represents the preconvolved environment map, while the second term (orange) represents the preintegrated BRDF.
Karis assumes the BRDF is based on the microfacet model:
He further assumes the Fresnel term is Schlick's simplified model:
That allows us to rewrite the orange BRDF integral as:
We see that the integral in green is exactly the irradiance table described in part 1 and if you implemented Karis's method then there's a strong chance you're already sampling this table, although only in the direction of the light.
The multiplescattering energy compensation method really just requires an additional sample of that table in the view direction... No big deal.
Far Field Environment
Since we made the assumption that the MSBRDF response to the far field environment is necessarily a low frequency phenomenon, we will only focus our efforts on computing the response of the MSBRDF to a lowfrequency encoding of the environment: this is why we will ignore cube map representations traditionnally used for (single scattering) specular reflections and directly use the SH representation of these cube maps.
Assuming we have a representation of the far field environment given as a set of spherical harmonic coefficients L_{lm}, Ramamoorthi and Hanrahan ^{3} showed that a simple 2^{nd} order representation (i.e. 9 coefficients) is enough to properly recover the irradiance (i.e. the integral of the environment's radiance over an entire hemisphere).
From these coefficients, we can obtain the bandlimited directional value of radiance:
Where N is the SH order.
We compute the response of the MSBRDF to the environment by writing:
We first notice that this can be easily split into N^2 integral expressions for each index of l, m:
Now, if we use the energy compensation formulation in the integral, complete with the roughness value \alpha, we get:
We see that the result is quite simple and is nicely split into 3 distinct parts:
 \frac{1  E(\boldsymbol{\omega_o}, \alpha)}{\pi  E_{avg}( \alpha )} has already been covered in part 1 and gives the viewdependent part of the equation
 L_{lm} are the environmentdependent SH coefficients (think diffuse environment probes here)
 E_{lm}( \boldsymbol{n}, \alpha ) = \int_{\Omega_+} Y_{lm}(\boldsymbol{\omega_i}) \left( 1  E(\boldsymbol{\omega_i}, \alpha)) \right) \left( \boldsymbol{\omega_i} \cdot \boldsymbol{n} \right) d\omega_i is certainly the most interesting part of all as it represents the encoding in SH of the MSBRDF response. We see that it's both view and environmentagnostic and can thus be precomputed once per type of BRDF (i.e. one for GGX, one for OrenNayar, etc.) but there is still a dependence on the surface normal orientation \boldsymbol{n} and the roughness coefficient \alpha that needs to be worked out.
Simplification
First, we begin by noticing that due to the isotropic nature of the MSBRDF, the MSBRDF is radially symmetric thus only the Zonal Harmonics coefficients should be non zero. This allows us to rewrite:
Furthermore, we could decide to either:
 Perform the integral over all possible values of \boldsymbol{n} (much too expensive!),
 Perform the integral in a single direction \boldsymbol{n} = (0,0,1) and rotate the environment light's L_{lm} coefficients to align them on the surface normal (better but still expensive to rotate the SH),
 Or simply perform the integral in a single direction \boldsymbol{n} = (0,0,1) and rotate the ZH coefficients to align them on the surface normal.
The 3^{rd} option is obviously much cheaper since ZH coefficients can easily be rotated into any direction to obtain a full set of SH coefficients, as given by Sloan ^{5}:
You can find the code for the ZH rotation below:
Implementation of order 2 ZH coefficients rotation in any direction (HLSL)
// Computes the Ylm coefficients in the requested direction // void Ylm( float3 _direction, out float _SH[9] ) { const float c0 = 0.28209479177387814347403972578039; // 1/2 sqrt(1/pi) const float c1 = 0.48860251190291992158638462283835; // 1/2 sqrt(3/pi) const float c2 = 1.09254843059207907054338570580270; // 1/2 sqrt(15/pi) const float c3 = 0.31539156525252000603089369029571; // 1/4 sqrt(5/pi) float x = _direction.x; float y = _direction.y; float z = _direction.z; _SH[0] = c0; _SH[1] = c1*y; _SH[2] = c1*z; _SH[3] = c1*x; _SH[4] = c2*x*y; _SH[5] = c2*y*z; _SH[6] = c3*(3.0*z*z  1.0); _SH[7] = c2*x*z; _SH[8] = 0.5*c2*(x*x  y*y); } // Rotate ZH cosine lobe into specific direction void RotateZH( float3 _A, float3 _wsDirection, out float _SH[9] ) { _A *= float3( 3.5449077018110320545963349666823, 2.0466534158929769769591032497785, 1.5853309190424044053380115060481 ); // Multiply by sqrt( 4 * PI / (2*l+1) ) as by eq. 26 in "On the relationship between radiance and irradiance" by Ramamoorthi Ylm( _wsDirection, _SH ); _SH[0] *= _A.x; _SH[1] *= _A.y; _SH[2] *= _A.y; _SH[3] *= _A.y; _SH[4] *= _A.z; _SH[5] *= _A.z; _SH[6] *= _A.z; _SH[7] *= _A.z; _SH[8] *= _A.z; }
So overall, we only need to precompute the ZH coefficients for various roughness values for all our BRDFs:
Below we see the appearance of each of the 3 ZH coefficients E_0( \alpha ), E_1( \alpha ) and E_2( \alpha ) in red, green and blue respectively for the GGX and OrenNayar MSBRDFs:
Fitting
We find a very close fit for the GGX coefficients:
As well as for the OrenNayar coefficients:
Validation
Here is the comparison for the GGX BRDF with the IOR for gold:
Comparison of "ground truth" against simplified spherical harmonics environment, for various values of roughness from 0.25 to 1. Only the multiplescattering component is shown.
Top Row: Ground truth sampling of the GGX MSBRDF.
Middle Row: A single sample of the GGX MSBRDF is used as well as a 2^{nd} order SH representation of the environment (i.e. 9 coefficients).
Bottom Row: Difference between the two, amplified 4x.
And the comparison for the OrenNayar BRDF with a blue reflectance:
Comparison of "ground truth" against simplified spherical harmonics environment, for various values of roughness from 0.25 to 1. Only the multiplescattering component is shown.
Top Row: Ground truth sampling of the OrenNayar MSBRDF.
Middle Row: A single sample of the OrenNayar MSBRDF is used as well as a 2^{nd} order SH representation of the environment (i.e. 9 coefficients).
Bottom Row: Difference between the two, amplified 4x.
The code used for these tests is as follows:
(HLSL)
float3 EstimateMSIrradiance_SH_GGX( float _roughness, float3 _wsNormal, float _mu_o, float3 _environmentSH[9], Texture2D<float> _Eo, Texture2D<float> _Eavg ) { // Estimate the MSBRDF response in the normal direction const float3x3 fit = float3x3( 0.01792303243636725, 1.0561278339405598, 0.4865495717038784, 0.06127443169094851, 1.3380225947779523, 0.6195823982255909, 0.10732852337149004, 0.8686198207608287, 0.3980009298364805 ); float3 roughness = sqrt( _roughness ) * float3( 1, _roughness, _roughness*_roughness ); float3 ZH = saturate( mul( fit, roughness ) ); float SH[9]; RotateZH( ZH, _wsNormal, SH ); // Estimate MSBRDF irradiance float3 result = 0.0; [unroll] for ( uint i=0; i < 9; i++ ) result += SH[i] * _environmentSH[i]; result = max( 0.0, result ); // Finish by estimating the viewdependent part float E_o = _Eo.SampleLevel( LinearClamp, float2( _mu_o, _roughness ), 0.0 ); float E_avg = _Eavg.SampleLevel( LinearClamp, float2( _roughness, 0.5 ), 0.0 ); result *= (1.0  E_o) / max( 0.001, PI  E_avg ); return result; } float3 EstimateMSIrradiance_SH_OrenNayar( float _roughness, float3 _wsNormal, float _mu_o, float3 _environmentSH[9], Texture2D<float> _Eo, Texture2D<float> _Eavg ) { // Estimate the MSBRDF response in the normal direction const float3x4 fit = float3x4( 0.0919559140506979, 1.467037714315657, 1.673544888379740, 0.607800523815945, 0.1136684128860008, 1.901273744271233, 2.322322430339633, 0.909815621695672, 0.0412482175221291, 1.093354950053632, 1.417191923789875, 0.581084435989362 ); float4 roughness = sqrt( _roughness ) * float4( 1, _roughness, _roughness*_roughness, _roughness*_roughness*_roughness ); float3 ZH = saturate( mul( fit, roughness ) ); float SH[9]; RotateZH( ZH, _wsNormal, SH ); // Estimate MSBRDF irradiance float3 result = 0.0; [unroll] for ( uint i=0; i < 9; i++ ) result += SH[i] * _environmentSH[i]; result = max( 0.0, result ); // Finish by estimating the viewdependent part float E_o = _Eo.SampleLevel( LinearClamp, float2( _mu_o, _roughness ), 0.0 ); float E_avg = _Eavg.SampleLevel( LinearClamp, float2( _roughness, 0.5 ), 0.0 ); result *= (1.0  E_o) / max( 0.001, PI  E_avg ); return result; } float3 EstimateMSIrradiance_SH( float3 _wsNormal, float3 _wsReflected, float _mu_o, float _roughnessSpecular, float3 _IOR, float _roughnessDiffuse, float3 _albedo, float3 _environmentSH[9] ) { // Estimate specular irradiance float3 F0 = Fresnel_F0FromIOR( _IOR ); float3 MSFactor_spec = F0 * (0.04 + F0 * (0.66 + F0 * 0.3)); // From http://patapom.com/blog/BRDF/MSBRDFEnergyCompensation/#varyingthefresnelreflectancef_0f_0 float3 Favg = FresnelAverage( _IOR ); float3 E_spec = MSFactor_spec * EstimateMSIrradiance_SH_GGX( _roughnessSpecular, _wsReflected, _mu_o, _environmentSH, _tex_GGX_Eo, _tex_GGX_Eavg ); // Estimate diffuse irradiance const float tau = 0.28430405702379613; const float A1 = (1.0  tau) / pow2( tau ); float3 rho = tau * _albedo; float3 MSFactor_diff = A1 * pow2( rho ) / (1.0  rho); // From http://patapom.com/blog/BRDF/MSBRDFEnergyCompensation/#varyingdiffusereflectancerhorho float3 E_diff = MSFactor_diff * EstimateMSIrradiance_SH_OrenNayar( _roughnessDiffuse, _wsNormal, _mu_o, _environmentSH, _tex_OrenNayar_Eo, _tex_OrenNayar_Eavg ); // Attenuate diffuse contribution float a = _roughnessSpecular; float E_o = _tex_GGX_Eo.SampleLevel( LinearClamp, float2( _mu_o, a ), 0.0 ); // Already sampled by MSBRDF earlier, optimize! float3 kappa = 1  (Favg * E_o + MSFactor_spec * (1.0  E_o)); return E_spec + kappa * E_diff; }
You then simply need to add a single call to "EstimateMSIrradiance_SH( ... )" at the end of your lighting pipeline, like this:
Retrieving the MSBRDF response from a SHencoded environment (HLSL)
// Here: // wsView is the worldspace camera vector pointing toward the camera // wsNormal is the worldspace normal vector // alphaS is the specular roughness // IOR is the index of refraction of the surface // alphaD is the diffuse roughness // rho is the diffuse reflectance in [0,1] // envSH is the SH representation of the surrounding environment (my environment values are usually filtered with a Hanning filter with window size 2.8) // #if USE_REFLECTED_VIEW // I just give that as an idea but I actually get worse results than using the plain normal float3 wsReflectedView = reflect( wsView, wsNormal ); float3 wsReflected = normalize( lerp( wsReflectedView, wsNormal, alphaS ) ); // Go more toward prefectly reflected direction when roughness drops to 0 #else float3 wsReflected = wsNormal; #endif sumIrradiance += EstimateMSIrradiance_SH( wsNormal, wsReflected, saturate( dot( wsView, wsNormal ) ), alphaS, IOR, alphaD, rho, envSH );
The entirety of the project is available at my God Complex GitHub repository.
Area Lights
In the paper "RealTime PolygonalLight Shading with Linearly Transformed Cosines" ^{1}, Heitz et al. introduced an important tool to represent easily integrable and configurable linear distribution transformations.
The main idea is to transport difficult integration situations involving an arbitrary distribution D(\boldsymbol{\omega}) into a canonical domain where the new original distribution D_o(\boldsymbol{\omega_o}) is easily characterized.
The paper specially focuses on the clamped cosine distribution D_o(\boldsymbol{\omega_o}) = \frac{1}{\pi} max( 0, z ), ~~ \boldsymbol{\omega_o}=(x,y,z).
The source area light polygon P is the domain where the arbitrary distribution D(\boldsymbol{\omega}) should initialy be integrated.
Instead, the transformed polygon P_o is integrated with the simpler canonical original distribution D_o(\boldsymbol{\omega_o})
The linearly transformed distribution has the nice property that the integration over polygons doesn't change:
And thus, we can write the classical lighting equation as:
Here, \boldsymbol{\omega_v} and \boldsymbol{\omega_l} respectively replace the vectors \boldsymbol{\omega_o} and \boldsymbol{\omega_i} used everywhere else in this dossier due to the choice of the o subscript used by Heitz et al. to denote the canonical basis.
Next, if we make the assumption that the radiance L(\boldsymbol{\omega}) is constant over the entire area light then we can write:
Where E(\boldsymbol{P_o}) is the irradiance over the area of polygon \boldsymbol{P_o}. Such irradiance has a closed form solution as given by Baum et al. ^{4} and the time cost grows linearly with the amount of edges of the polygon.
Basically then, all that is required is to transform the vertices of the area light polygon \boldsymbol{P} into the canonical domain using a precomputed transform M so that \boldsymbol{P_o} = M^{1} \boldsymbol{P}, then use the analytical expression for the irradiance.
The essential work of this technique resides in finding the proper value of the transform matrix M depending on our BRDF. Heitz et al. provide the source code they used to fit the GGX BRDF for various values of elevation \cos(\theta_o) and surface roughness \alpha that they stored into 2 textures (because the matrix M is described by 5 values).
ViewDependence
We take model on our work for the SHencoded environment to once again separate the viewdependent part out of the integral:
So the part we need to fit with a clampedcosine distribution is only this expression:
What this means for us is that, similarly to regular BRDFs, we will have to precompute transform matrices for our MSBRDF but also that:
 There is no need to store a 2D table since we factoredout the viewdependent part of the BRDF and only a dependence on roughness \alpha remains
 Due to the isotropic and smooth nature of the MSBRDF, there shouldn't be a need to handle anisotropy or complicated coefficients so the matrix M should be simpler
Fitting
Todo
Fit MSBRDF GGX / Oren MSBRDF
Use the fact that our MSBRDFs are isotropic. Should require less coefficients for M!
Can we also factor out the viewdependent part of regular BRDFs? (i.e. that would be the masking term \frac{G( \boldsymbol{\omega_v} )}{\boldsymbol{\omega_v} \cdot \boldsymbol{n}}). Does it yield a smoother, more accurate result?
Validation
Todo
Show the "ground truth" (i.e. many environment rays) against a single precomputed estimate...
Simulation de ClearCoated Diffuse Surface
Todo
TODO Pour voir si mon approximation tient la route...
References

Heitz, E. Dupuy, J. Hill, S. Neubelt, D. 2016 "RealTime PolygonalLight Shading with Linearly Transformed Cosines" ↩↩

Karis, B. 2013 "Real Shading in Unreal Engine 4" ↩

Ramamoorthi, R. Hanrahan, P. 2001 "On the relationship between radiance and irradiance: determining the illumination from images of a convex Lambertian object" ↩

Baum, D. R. Rushmeier, H. E. Winget, J. M. 1989 "Improving radiosity solutions through the use of analytically determined formfactors" ↩

Sloan, PP. 2008 "Stupid Spherical Harmonics Tricks" ↩