Energy Compensation Pt. 1
In their 2017 talk ^{1} at the now famous Siggraph's Physically Based Shading in Theory and Practice courses, Kulla and Conti reintroduced a computation ^{2} devised long ago by Kelemen and SzirmayKalos in 2001.
The Original Paper
In section 2.2 of their paper ^{2}, Kelemen et al. wrote about coupling the matte (i.e. diffuse) and specular parts of the BRDF.
They write the complete BRDF as: $$ f_r(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) = f_{r,spec}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) + f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) $$
And claim that, although f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) is difficult to estimate due to the many scattering events that occur when light is not specularly reflected but rather diffused through the material, they can safely wager about the fact that f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) is:
 Energypreserving
 Symmetrical
 Somehow, the complement (that's the operative word here) of the specular f_{r,spec}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) part
I remember being very impressed by the images produced by this paper by the time it was published (yes! I'm old!).
I believe even now there is a very strong "ground truth" flavor emanating from these images.
Prototyping the diffuse BRDF
They write the prototype for such an (isotropic) BRDF as:
Where:
 k(\lambda) is a wavelengthdependent factor (namely, the surface's reflectance in [0,1] for different wavelengths)
 s is a normalization factor yet to be determined
 r(\mu) is some unknown "appropriate scalar function", yet to be determined too
 \mu_i and \mu_o are the \boldsymbol{\omega_i}\cdot\boldsymbol{n} = \cos(\theta_i) and \boldsymbol{\omega_o}\cdot\boldsymbol{n} = \cos(\theta_o) respectively, \boldsymbol{n} being the surface's normal
Solving for unknowns
Kelemen et al. continue by writing the albedo (i.e. total reflectance for a particular viewing direction) for such a diffuse material:
Note
Nowadays, we call this total reflectance integral the "white furnace test" as it simply integrates the BRDF against a unit radiance over the entire hemisphere.
The albedo can thus simply be viewed as a measure of irradiance against a totally white ambient background and we will now write:
(E being the symbol usually used for the irradiance)
Since E_{diff}(\mu_o) + E_{spec}(\mu_o) \le 1, we can conclude that necessarily: $$ E_{diff}(\mu_o) \le 1E_{spec}(\mu_o) $$
Moreover, in the perfectly reflecting case where the total albedo E_{diff}(\mu_o) + E_{spec}(\mu_o) = 1 and k(\lambda)=1 then strictly: $$ E_{diff}(\mu_o) = 1E_{spec}(\mu_o) \tag{3}\label{(3)} $$
Equation \eqref{(2)} shows that the diffuse albedo is proportional to r(\mu_o) and, symmetrically, r(\mu_o) is thus proportional to E_{diff}(\mu_o) = 1E_{spec}(\mu_o).
The important takeway remark here is that:
For the perfectly reflecting case, we can thus rewrite equation \eqref{(1)} as: $$ f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) = s.(1E_{spec}(\mu_o)).(1E_{spec}(\mu_i)) $$
Substituting \eqref{(3)} and this new BRDF into \eqref{(2)} we get:
Where the specular albedo averaged over all possible view directions on the hemisphere is represented by:
Warning
Notice that in the Kelemen and Kulla notations for f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}), they factorized the \pi out of the denominator so they write it as \pi \left( 1  E_{spec}^{avg} \right) and they have E_{spec}^{avg} = 2 \int_0^1 E_{spec}(\mu_i) \mu_i d\mu_i instead but I find that highly disturbing so I didn't follow their example (I like to imagine the E_{spec}^{avg} integral converging to a maximum of \pi instead of obfuscating that fact for the sake of a "nicer way of writing the result").
Finally we write: $$ f_{r,diff}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) = \frac{(1E_{spec}(\mu_o)).(1E_{spec}(\mu_i))}{\pi  E_{spec}^{avg}} \tag{5}\label{(5)} $$
Proof of concept
Armed with this new expression for the diffuse BRDF, what happens if we integrate against a unit radiance over the entire hemisphere (i.e. the white furnace test again)?
We see that E_{diff}(\mu_o) ends up being the exact complement of E_{spec}(\mu_o)!
So Kelemen et al. already had the key in 2001 but they apparently failed to notice the importance of their result (or did they? ).
The Revised Usage
After all, isn't that result what we're looking for when looking to compute the multiplescattering term? Missing energy due to a singlescattering BRDF term that is often too simple?
This is indeed exactly what Kulla et al. very cleverly noticed in their new way of viewing of the problem!
Let E_{diff}(\mu_o) get rewritten as E_{ms}(\mu_o) instead and we get: $$ E_{ms}(\mu_o) = 1  E(\mu_o) $$
With E_{ms}(\mu_o) the irradiance from the multiplyscattered BRDF and E(\mu_o) the irradiance from our classical singlescattered BRDF.
And there you have it:
And we have our new expression for the (ideally reflecting) multiplescattering BRDF: $$ f_{r,ms}(\boldsymbol{\omega_o},\boldsymbol{\omega_i}) = \frac{(1E(\mu_o)).(1E(\mu_i))}{\pi  E_{avg}} \tag{6}\label{(6)} $$
White furnace test without and with energy compensation for various values of surface roughness.
(the background environment intensity is 1 while the scene intensity is 0.9 to make sure there is no overshoot in energy reflection).
Top Row: The classical singlescattering BRDF clearly fails at restituting the entire energy received by the material
Middle Row: The multiplescattering BRDF term from the energy compensation technique successfully bounces back the entirety of the energy received by the material.
Bottom Row: Difference.
Applications to existing BRDF models
GGX Specular Model
We use the (now) very common GGX normal distribution and Smith GGX shadowing/masking term:
With F_0 the reflectance at normal incidence, \alpha \in [0,1] the surface roughness and \boldsymbol{\omega_h} the normalized halfvector.
Rewritten in terms of \mu we have:
With:
 \mu_d = \boldsymbol{\omega_i} \cdot \boldsymbol{\omega_h} the cosine of the angle between the light and the half vector
 \mu_h(\mu_o, \mu_i, \phi) = \boldsymbol{\omega_h} \cdot \boldsymbol{n} = \frac{ \mu_o + \mu_i } { \sqrt{ 2 \left( 1 + \mu_o \mu_i + (1\mu_o^2)(1\mu_i^2) \cos(\phi) \right) } }
 \phi the azimutal angle between \boldsymbol{\omega_o} and \boldsymbol{\omega_i}.
Irradiance Table
With F_0 = 1 \Rightarrow F( \mu_d, F_0 ) = 1 at the moment (i.e. perfectly reflective case), we precompute the "complement albedo table" for all possible viewing angle \mu_o=\cos(\theta_o) and all roughness values \alpha for the specular BRDF:
You can see the resulting table below:
Warning
Obviously, don't use this awful image directly but this 128x128 table instead! (I provide a 128x128 version although, as noted by Kulla et al., the function is very smooth and a 32x32 texture is more than enough).
The 1^{st} float is \mu = \cos(\theta) of the incident or outgoing ray direction, the 2^{nd} float is the roughness \alpha and the 3^{rd} float is E\left( \mu, \alpha \right) (caution, not 1E!)
If you're using the correlated Smith Shadowing/Masking function instead of the uncorrelated one described by the BRDF equation above, please use this table instead.
Average Irradiance
Armed with this table, we can obtain the average irradiance table that only depends on roughness by computing \eqref{(4)}:
Info
You can download this table representing the E_{avg} for different values of roughness.
The 1^{st} float is the roughness \alpha and the 2^{nd} float is E_{avg}\left( \alpha \right)
Alternatively, you can use the following analytical fit: E_{avg}(\alpha) = \pi  0.446898 \cdot \alpha  5.72019 \cdot \alpha^2 + 6.61848 \cdot \alpha^3  2.41727 \cdot \alpha^4
If you're using the correlated Smith Shadowing/Masking function instead of the uncorrelated one described by the BRDF equation above, please use this table instead.
Energy Conservation Check
We quickly check the energy is conserved by ensuring that for all roughness values of \alpha \in [0,1] we have:
Gray curve is the GGX specular BRDF, blue curve is the "energy compensation BRDF", red curve is their sum that always yield \pi, thus ensuring the conservation of energy.
We also quickly notice that the multiple scattering BRDF term becomes preponderant over the single scattering term when \alpha > 0.8, so for very rough materials where shadowing and masking are playing a major
role in energy loss.
The case of perfectly reflective rough metal
Remembering that we fixed the Fresnel term to be F( \mu_d, F_0 ) = 1, the tables we just calculated can only give us the perfectly reflective 100% white metal BRDF case:
Top Row: The classical singlescattering GGX BRDF shows a loss of energy at high roughness
Middle Row: The multiplescattering BRDF term with the energy compensation technique successfully restores the lost energy.
Bottom Row: Difference.
OrenNayar Diffuse Model
We now concentrate on another BRDF, the OrenNayar diffuse model.
First introduced in 1992 by Michael Oren and Shree Nayar in the context of machine vision^{3}, then published again^{4} in 1994 for CGI this time, this model is similar to the CookTorrance microfacet model used for specular reflection as it uses vcavities to represent the roughness of a surface but this time for the diffuse case (i.e. the microfacets are purely lambertian, not pure mirrors like with the CookTorrance model), the Lambertian reflectance being only a special configuration of this more general model.
Oren and Nayar also provide an additional interreflection term, as a complement to the single scattering BRDF but it is finally not included in their formulation.
It's this additional term that I propose to provide here using the multiplescattering integral instead.
In section 4.4 of their paper, we find the expression for the BRDF of a generalized rough diffuse surface that is written as:
With \rho the surface reflectance and \sigma \in [0,\frac{\pi}{2}] the standard angle deviation for the microfacets' slope distribution that represents the roughness of the surface. Note that \sigma = 0 falls back to the standard Lambertian reflection.
You can find below a simple HLSL implementation for the OrenNayar diffuse model:
OrenNayar Implementation (HLSL)
// Simple OrenNayar implementation // _normal, unit surface normal // _light, unit vector pointing toward the light // _view, unit vector pointing toward the view // _roughness, OrenNayar roughness parameter in [0,PI/2] // float BRDF_OrenNayar( in float3 _normal, in float3 _view, in float3 _light, in float _roughness ) { float3 n = _normal; float3 l = _light; float3 v = _view; float LdotN = dot( l, n ); float VdotN = dot( v, n ); float gamma = dot( v  n * VdotN, l  n * LdotN ) / (sqrt( saturate( 1.0  VdotN*VdotN ) ) * sqrt( saturate( 1.0  LdotN*LdotN ) )); float rough_sq = _roughness * _roughness; float A = 1.0  0.5 * (rough_sq / (rough_sq + 0.33)); // You can replace 0.33 by 0.57 to simulate the missing interreflection term, as specified in footnote of page 22 of the 1992 paper float B = 0.45 * (rough_sq / (rough_sq + 0.09)); // Original formulation // float angle_vn = acos( VdotN ); // float angle_ln = acos( LdotN ); // float alpha = max( angle_vn, angle_ln ); // float beta = min( angle_vn, angle_ln ); // float C = sin(alpha) * tan(beta); // Optimized formulation (without tangents, arccos or sines) float2 cos_alpha_beta = VdotN < LdotN ? float2( VdotN, LdotN ) : float2( LdotN, VdotN ); // Here we reverse the min/max since cos() is a monotonically decreasing function float2 sin_alpha_beta = sqrt( saturate( 1.0  cos_alpha_beta*cos_alpha_beta ) ); // Saturate to avoid NaN if ever cos_alpha > 1 (it happens with floatingpoint precision) float C = sin_alpha_beta.x * sin_alpha_beta.y / (1e6 + cos_alpha_beta.y); return A + B * max( 0.0, gamma ) * C; }
Irradiance Table
Once again, keeping \rho = 1 at the moment (i.e. perfectly reflective case), we compute the irradiance table E(\mu_o,\alpha) using the OrenNayar BRDF model, with \sigma = \frac{\pi}{2} * \alpha and we obtain the following table:
Info
The function is so smooth that this time I only give a 32x32 table with the same formatting as the GGX tables.
Average Irradiance
And again, we can obtain the average irradiance table that only depends on roughness:
Info
You can download this table representing the E_{avg} for different values of roughness.
Alternatively, you can use the following analytical fit: E_{avg}(\alpha) = \pi  7.49231 \cdot \alpha^2 + 11.4409 \cdot \alpha^3  5.05903 \cdot \alpha^4
Energy Conservation Check
We quickly check the energy is conserved by ensuring that for all roughness values of \alpha \in [0,1] we have:
Gray curve is the OrenNayar diffuse BRDF, blue curve is the "energy compensation BRDF", red curve is their sum that always yield \pi, thus ensuring the conservation of energy.
The case of perfectly reflective rough diffuse
Remembering that we fixed the diffuse reflectance term to be \rho = 1, the tables we just calculated can only give us the perfectly reflective 100% white diffuse BRDF case:
Top Row: The classical singlescattering OrenNayar BRDF shows a loss of energy at high roughness
Middle Row: The multiplescattering BRDF term with the energy compensation technique successfully restores the lost energy.
Bottom Row: Difference.
Dealing with Varying Reflectances
Up until this point we have made 2 important assumptions regarding the BRDFs we have seen:
 The Fresnel reflectance F( \mu_d, F_0 ) = 1 for the GGX BRDF
 The diffuse reflectance \rho = 1 for the OrenNayar BRDF.
What this means is that the total irradiance that is lost when the light hits the microsurface multiple times, which is \pi  E_{avg}, is entirely redistributed into the multiplescattering term.
Obviously, with real life materials the light gets more or less transmitted through the surface and continues its path:
 For metals, most of the photons get absorbed, transformed into electrons and they never exit the surface as photons again (or maybe super weakened as infrared? Not sure about that).
 For dielectrics:
 If the material has a crystaline structure then photons just don't interact with it unless they have sufficient energy to change the orbits of electrons (that's quantum mechanic domain here)
 For regular materials though (or crystals with impurities), it's a mixture of scattering and absorption events not unlike what happens in a participating medium but here the medium density
is extremely high and extinction happens very quickly.
Nevertheless after many scattering events, the photons adopt a rather isotropic distribution (hence the mainly diffuse aspect of reflected light) and some of it may exit the surface again, maybe at the same location it entered the surface (the BRDF hypothesis) but more plausibly at a nearby location (the BSSRDF hypothesis), and possibly even through the surface (i.e. translucency)
In order to apply the multiplescattering part of the BRDF, we need to know how much is actually redistributed in reflection depending on the reflectance associated to the BRDF, be it the Fresnel reflectance for the specular BRDFs, or the diffuse reflectance for the diffuse BRDFs.
Varying the Fresnel Reflectance F_0
Obviously, when the F_0 term is not 1 anymore then the tiny microfacets composing the microscopic surface stop being perfect mirrors: some of the incoming energy gets reflected off the surface (in green), and the remaining energy gets refracted below the surface (in blue).
The macroscopic result of combining many microfacets gives a certain amount of the incoming energy that gets reflected (in green), and a certain amount gets transmitted (in blue).
BRDF models provide a way of estimating the reflected and refracted amounts of energy based on a statistical analysis of the properties of idealized microsurfaces.
Inset: The ideal microfacet mirror configuration, with a perfect reflection direction (in green) and perfect refraction direction (in blue).
We are thus looking for a factor f(F_0) \in [0,1] that will be applied to the multiple scattering BRDF term, the complement of this factor 1  f(F_0) times the MSBRDF term should be redistributed
to the part below the specular interface (e.g. the diffuse BRDF).
In order to compute such a factor I used my heavyduty microfacet raytracer to simulate the total irradiance outgoing the microsurface for various values of incident angle, surface roughness and F_0.
Note
I verified that the dielectric and metallic Fresnel expressions almost match for all possible F_0 so the provided expression works regardless of the nature of the material,
as can be seen below where the blue curve is the "ArtistFriendly Metal Fresnel" formulation by Ole Gulbrandsen^{5} and the red curve is the classical dielectric Fresnel formulation
given by Cook & Torrance^{6}.
We see the curves almost completely match as long as the grazing angle reflectance F_{90}=1. The dielectric formulation is obviously prefered for low values of F_0 < 0.2:
Then, I performed the white furnace integration for scattering orders from 2 to 6. So basically what we have is the experimental data that allows us to compute this expression for each scattering order:
Note
The white furnace from experimental data (in red, rescaled by a factor 2.1 to match) differs quite a lot from the one we got from the analytical GGX integral (in blue). Despite the factor 2.1 in amplitude, the shape of the experimental white furnace test is similar to the one from the analytical GGX integral but not quite the same, I believe that is expected because I used a Beckmann distribution for the simulation heightfield^{7} and not a Smith distribution.
Anyway, we will assume the anaytical formulation behaves in the same manner as the experimental data even though the resulting shapes and amplitudes don't quite match... We are essentially interested in the relative relationships between various scattering orders rather than the absolute magnitude of the results.
In any case, the total energy is conserved and the sum of experimental white furnace tests for each order > 2 (in red) is the exact complement of the singlescattered order 1 (in blue), their sum yielding the constant radiated energy \pi (in green).
For each scattering order, we are looking for a polynomial of the form a_n(F_0) = A_n \cdot F_0^{n'}, with n' a variation on the scattering order and A_n an amplitude factor depending on the order.
It's not as easy as with a diffuse BRDF because F_0 is not a linear value used by the Fresnel equation but is instead converted into an index of refraction (IOR) given by:
And indeed, by using the expression:
We obtain quite a good fit as can be seen below where each curve represents the white furnace integral as a function of specular reflectance F_0, the red curve is the polynomial fit for the scattering order:
By cheating a little more, we know the A_n terms should have the form of a decreasing geometric series A_n = A_1 \cdot \left(\sqrt{\tau}\right)^{n^2} where A_1 is the initial factor and \tau < 1 the factor to apply to the specular reflectance so that it matches our experimental data.
After fitting, we indeed obtain the value for \tau:
We then obtain the following complicated power series that we require to sum to 1:
It seems such series are special kinds of functions called Jacobi theta functions and especially the Jacobi elliptic θ_3(z,q) function with the specific value θ_3(0,\sqrt{\tau}).
Fortunately, with \sqrt{\tau} < 1, the series converges and we can write:
And thus:
We see that we have a somewhat okay fit to the white furnace tests for each order using our newlyfitted function a_n(\rho) = A_1 \left(\sqrt{\tau F_0}\right)^{n^2}:
Our remaining issue is to compute the \theta_3\left(0,\sqrt{\tau F_0}\right) function for various values of F_0. Fortunately, the theta function is very smooth and can be easily fitted using regular loworder polynomials.
And the final factor to apply to the GGX multiple scattering BRDF is thus simply:
Which gives this pretty uninteresting function that nevertheless makes for a nice saturation in color due to its nonlinear nature (blue is the theta function, red is fitting):
From Left to Right: Varying values of F_0 from 0.04 to F_{0_{gold}}=\left( 1, 0.765557, 0.336057 \right) with constant roughness of 0.5.
Previous to last image shows no sturation on F_0 (i.e. F_0 is used directly as a multiplier for the MS term).
Last image on the right shows amplified difference of the saturated colors obtained with the F_{ms}(F_0) discussed above.
Varying Diffuse Reflectance \rho
Identically to Fresnel reflectance, when the diffuse reflectance term \rho is not 1 anymore then the tiny microfacets composing the microscopic surface stop being perfect lambert reflectors: some of the incoming energy gets reflected off the surface (in green), and the remaining energy gets absorbed by the surface (in blue).
Identical to the specular case, the diffuse macrosurface reflects and transmits a certain amount of the energy it receives.
Inset: The ideal microfacet lambert reflector configuration, with a perfectly lambertian probability of reflecting the incoming energy (in green) and a simplified absorption model where the energy simply "disappears" (in blue).
This time, my heavyduty microfacet raytracer used various values of incident angle, surface roughness and albedo, and I performed the white furnace integration for scattering orders from 2 to 6. So basically what we have is the experimental data that allows us to compute this expression for each scattering order:
Note
Once again, the white furnace from experimental data (in red) differs from the one we got from the analytical OrenNayar integral (in blue) by a factor 1.17 and not quite the same shape, but I guess that's to be expected from an approximate model and an empirical simulation?
Anyway, we will assume once again the anaytical formulation behaves in the same manner as the experimental data even though the resulting shapes and amplitudes don't quite match...
By cheating a little and knowing full well we should expect a polynomial of the form a_n(\rho) = A_n \cdot \rho^n, with n the scattering order, we obtain an excellent fit indeed, as can be seen below where each curve represents the white furnace integral as a function of surface albedo \rho, the red curve is the polynomial fit for the scattering order:
By cheating even more, we know the A_n terms should have the form of a decreasing geometric series A_n = A_1 \cdot \tau^n where A_1 is the initial factor and \tau < 1 the factor to apply to the diffuse reflectance so that it matches our experimental data.
After fitting, we first obtain the value for tau:
Now we are looking for the global factor A_1 so that the sum of all the amplitudes over all the scattering orders to sum to 1 when \rho=1:
Since \tau < 1 we are dealing with a wellbehaved geometric series that can be simplified into:
We see that we have a nice fit to the white furnace tests for each order using our newlyfitted function a_n(\rho) = A_1 \left(\tau \rho\right)^n:
And the final factor to apply to the OrenNayar multiple scattering BRDF is thus:
Which gives this pretty uninteresting function that nevertheless makes for a nice saturation in diffuse color due to its nonlinear nature:
From Left to Right: Varying values of \rho from 0.2 to 1 with constant roughness of 1.
Last image on the right shows no sturation on \rho (i.e. \rho is used directly as a multiplier for the MS term instead of using the F_{ms}(\rho) term discussed above).
Complete Approximate Model
We can finally put together what we saw into an approximate model that should be good enough for real time game engine (maybe not so much for productiongrade rendering ).
For a metal, we can easily obtain the specularly reflected radiance as:
Where dE_i(\boldsymbol{\omega_i}) = L_i(\boldsymbol{\omega_i}) (\boldsymbol{\omega_i} \cdot \boldsymbol{n}) d\omega_i is the small irradiance element contributed by the projected solid angle (\boldsymbol{\omega_i} \cdot \boldsymbol{n}) d\omega_i.
As for dielectrics, the expression is a little bit more involved:
The main issue here is to determine the factor \kappa( \boldsymbol{\omega_o}, F_0 ) that gives the ratio of the incoming irradiance that is not specularly reflected.
If we remember the assumption originally made by Kelemen in equation \eqref{(3)}, we had in the perfectly reflecting case:
We know the diffuse irradiance E_{diff} is bounded by the complement to the specular irradiance E_{spec}, i.e. the diffusely reflected energy is what's left of the energy that hasn't been specularly reflected.
We can thus choose \kappa( \boldsymbol{\omega_o}, F_0 ) = E_{diff}( \boldsymbol{\omega_o} ) as our diffuse factor, all we need to find is the expression for E_{spec}( \boldsymbol{\omega_o} ) that we can now fully write as:
With E(\boldsymbol{\omega_o}) the precomputed GGX irradiance table we saw earlier.
If we make the (wrong) assumption that:
Then we finally have:
We can then use the expression for F_{avg}(\theta_o, F_0) given by Kulla & Conty (slide 18) for dielectrics that makes an additional simplification by ignoring the dependance on \theta_o:
Where \eta = \frac{1+\sqrt{F_0}}{1\sqrt{F_0}} is the relative Index Of Refraction (IOR) of the surface.
We can see that:
 \kappa( \boldsymbol{\omega_o}, F_0 ) \to 0 when F_0 \to 1, which is the case of metals
 \kappa( \boldsymbol{\omega_o}, F_0 ) \to 1 when F_0 \to 0, which is the case of dielectrics
 \kappa( \boldsymbol{\omega_o}, F_0 ) = 1, for perfectly diffuse surfaces
The full dielectric model with varying values for the specular roughness \alpha and reflectance F_0
Here is the code for the full MSBRDF model that you can also find in my God Complex repository on GitHub.
Full MSBRDF Model Implementation (HLSL)
static const bool _enableMultipleScattering = true; // ====== Multiplescattering BRDF computed from energy compensation ====== // _roughness, roughness for the BRDF model (in [0,1]) // _tsNormal, unit surface normal vector (in tangent space) // _tsView, unit view vector, pointing toward the camera (in tangent space) // _tsLight, unit light vector, pointing toward the light (in tangent space) // _Eo, the precomputed 2D table of values for various elevations and roughnesses // _Eavg, the precomputed 1D tables of total reflectance for various roughnesses // float MSBRDF( float _roughness, float3 _tsNormal, float3 _tsView, float3 _tsLight, Texture2D<float> _Eo, Texture2D<float> _Eavg ) { float mu_o = saturate( dot( _tsView, _tsNormal ) ); float mu_i = saturate( dot( _tsLight, _tsNormal ) ); float a = _roughness; float E_o = 1.0  _Eo.SampleLevel( LinearClamp, float2( mu_o, a ), 0.0 ); // 1  E_o float E_i = 1.0  _Eo.SampleLevel( LinearClamp, float2( mu_i, a ), 0.0 ); // 1  E_i float E_avg = _Eavg.SampleLevel( LinearClamp, float2( a, 0.5 ), 0.0 ); // E_avg return E_o * E_i / max( 0.001, PI  E_avg ); } // ====== Computes the full dielectric BRDF model ====== // _tsNormal, unit surface normal vector (in tangent space) // _tsView, unit view vector, pointing toward the camera (in tangent space) // _tsLight, unit light vector, pointing toward the light (in tangent space) // _roughnessSpecular, roughness for the GGX specular model (in [0,1]) // _IOR, surface Index Of Refraction for the GGX specular model (in [1,+oo]) // _roughnessDiffuse, roughness for the OrenNayar diffuse model (in [0,1]) // _albedo, diffuse albedo for the OrenNayar diffuse model (in [0,1]) // float3 ComputeBRDF_Full( float3 _tsNormal, float3 _tsView, float3 _tsLight, float _roughnessSpecular, float3 _IOR, float _roughnessDiffuse, float3 _albedo ) { // Compute specular BRDF float3 F0 = Fresnel_F0FromIOR( _IOR ); float3 MSFactor_spec = (_flags & 2) ? F0 * (0.04 + F0 * (0.66 + F0 * 0.3)) : F0; // From http://patapom.com/blog/BRDF/MSBRDFEnergyCompensation/#varyingthefresnelreflectancef_0f_0 float3 Favg = FresnelAverage( _IOR ); float3 BRDF_spec = BRDF_GGX( _tsNormal, _tsView, _tsLight, _roughnessSpecular, _IOR ); if ( _enableMultipleScattering ) { BRDF_spec += MSFactor_spec * MSBRDF( _roughnessSpecular, _tsNormal, _tsView, _tsLight, _tex_GGX_Eo, _tex_GGX_Eavg ); // _tex_GGX_Eo, _tex_GGX_Eavg are precomputed textures with the tables given earlier } // Compute diffuse contribution float3 BRDF_diff = _albedo * BRDF_OrenNayar( _tsNormal, _tsView, _tsLight, _roughnessDiffuse ); if ( _enableMultipleScattering ) { const float tau = 0.28430405702379613; const float A1 = (1.0  tau) / pow2( tau ); float3 rho = tau * _albedo; float3 MSFactor_diff = (_flags & 2) ? A1 * pow2( rho ) / (1.0  rho) : rho; // From http://patapom.com/blog/BRDF/MSBRDFEnergyCompensation/#varyingdiffusereflectancerhorho BRDF_diff += MSFactor_diff * MSBRDF( _roughnessDiffuse, _tsNormal, _tsView, _tsLight, _tex_OrenNayar_Eo, _tex_OrenNayar_Eavg ); // _tex_OrenNayar_Eo, _tex_OrenNayar_Eavg are precomputed textures with the tables given earlier } // Attenuate diffuse contribution float mu_o = saturate( dot( _tsView, _tsNormal ) ); 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 BRDF_spec + kappa * BRDF_diff; }
Discussion
We saw interesting improvements on classical BRDF formulations that can be pretty cheap to evaluate, even for game engines.
Unfortunately, multiplescattering BRDFs only reveal their true potential when dealing with advanced light sources like area lights or imagebased lighting (i.e. lights that cover a sufficiently large solid angle),
otherwise the small contribution of the MS term is quite "lost in the wind" and the cost of sampling 3 (even super low resolution) precomputed textures per pixel and per light can rapidly become prohibitive.
That is why we will discuss how the multiplescattering BRDF can be truly harnessed for a real added benefit in the 2^{nd} part of this series about energy compensation! (coming soon) (hopefully )
References

Kulla, C. Conty, A. 2017 "Revisiting Physically Based Shading at Imageworks" ↩

Kelemen, C. SzirmayKalos, L. 2001 "A Microfact Based Coupled SpecularMatte BRDF Model with Importance Sampling" ↩↩

Oren, M. Nayar, S. 1992 "Generalization of the Lambertian Model and Implications for Machine Vision" ↩

Oren, M. Nayar, S. 1994 "Generalization of Lambert's Reflectance Model" ↩

Gulbrandsen, O. 2014 "Artist Friendly Metallic Fresnel" ↩

Cook, R. L. Torrance, K. E. 1981 "A Reflectance Model for Computer Graphics" ↩

Heitz, E. 2014 "Generating Procedural Beckmann Surfaces" ↩