MultipleScattering BRDF
In 2014 and 2015, many groundbreaking papers from Heitz et al. ^{1} ^{2} ^{3} ^{4} started to shed a new light (haha) on alternative models to the classical microfacet model and introduced a more general volumetric framework called the "microflakes model" that tried to unify surfacic and volumetric materials into a unique solution (and indeed, in the real world, there is no such distinction between a "surface" and a "volume").
With this new microflakes model and adequate phase functions, Heitz et al. began to reflect (haha) on the reintroduction of the energy lost with classical BRDF models because of shadowing and masking.
Indeed, although the shadowing term correctly prevents a microsurface element to be directly lit by a light source, the masking on the other hand is incorrectly removing lighting that, in the real world, would otherwise be indirectly perceived by reflection.
Suffice it to say that classical BRDF models up until 2014 were singlescattering models, while in reality, light bounces of the microstructures of a surface/volume and is reinjected multiple times, showing a new gain in energy: realworld BRDFs are multiplescattering models.
(Source: "MultipleScattering Microfacet BSDFs with the Smith Model" ^{4})
Simulation
As soon as I had read the papers from Heitz et al. I immediately got to work to try and be one of the first to get a working MSBRDF model for runtime game use.
My idea was simple:
 Simulate the probabilistic reflection of many rays of light hitting a rough surface (first the diffuse case, then the metal case, then the general dielectric case with refraction)
 Analyze the shape of the lobe that reflects off from it
 Fit a simple analytical lobe model with few free parameters
 Find an analytical fit of these free parameters
 Use this analytical fit to express the lobes at runtime for many light bounces
 Add this quantity to the existing BRDF to obtain additional, multiplyscattered bounces
I was planning on having enough time to write 2 methods: one that uses many rays, the other one implementing the statistical model from Heitz et al. ^{4} and also to be able to simulate diffuse, metallic and dielectric surfaces.
In practice, I only had time to write the bruteforce raycasting method for diffuse, metallic and dielectrics, and do the fitting for the diffuse lobes only.
You can find the latest stage of the project from late 2015, early 2016 as I left it in my God Complex Repository.
Using Brute Force RayCasting
So I went on writing a small test application that would be using Compute Shaders to cast many rays — actually 500 million rays — on a tiny patch of rough microsurface (the size of the patch is not relevant, only the distribution of microfacet slopes is important).
My little simulator, bombarding 500 million rays on a tiny surface patch and accounting for 4 orders of scattering.
Basically, the algorithm when no refraction is involved goes like this:
 Create a random ray coming from the userspecified direction (i.e. an incidence angle \theta)

For each scattering event S
3.1. Shoot the ray across a carefully generated heightfield whose height distribution obeys a Beckmann distribution ^{9}
3.2. If the ray exits the surface, store it in the S order scattering histogram bin and exit
3.3. Else, reflect the ray across the perfectly specular microsurface and continue the trace
Here is a view of the resulting histogram bins for the 4^{th} order of scattering over a metallic surface of roughness \alpha = 0.8:
Each pixel represents a bin for a directional vector with spherical coordinates (\theta,\phi), the image covers the entire upper hemisphere (or lower hemisphere when we are dealing with refraction).
Lobe Model
I went on and implemented several microfacet lobe models, a model being the NDF and masking/shadowing terms. The Fresnel term is the same one for all the models.

Beckmann Model, that follows the classical Beckmann distribution

GGX Model, that follows the GGX distribution described by Walter et al. ^{8}

Modified Phong (isotropic and anisotropic), that follows the classical Phong distribution
Note
It is important to understand that these models are purely analytical models of normal distributions and shadowing/masking that were devised to work in the limited framework of the singlescattering CookTorrance microfacet model.
Customization
The goal of my application was to fit new lobes from simulated, empirical data thus I wasn't shy of adding new alien parameters to these models in order to make the resulting lobes more "bendable" and capable of fitting a larger set of shapes that could come out from multiple bounces of light, as can be seen in this image:
Multiple orders of scattering on a diffuse surface with roughness \alpha = 0.8. Lobes are scaled to roughly the same size each time, otherwise their volume collapses to 0 very rapidly with each new bounce.
We see that the diffuse lobes can be pretty "squashed". We get even worse kinds of shapes when dealing with dielectric materials.
So in addition to the default parameters of these "classical" models:
 The lobe's deviation angle \theta, giving its deviation from the macroscopic normal axis Z
 The roughness \alpha, in [0,1]
I added:
 A "global scale" factor in [0,\infty] that is applied to the entire lobe
 An "anisotropic flattening" factor in [1,1] that allows to squash the lobes along the tangential X & Y axes
 A "masking importance" factor in [0,1] that allows to ignore the influence of the shadowing/masking term entirely. It turns out this parameter is not very significant after all.
Lobe Fitting
The raycasting was quite fast (about ½ a second per simulation, for 4 scattering orders), now I needed to fit the lobe models with their 5 free parameters.
I revamped the BFGS implementation I wrote a few years back (it can be found here).
This BFGS minimizer is really easy to use: feed it a distance function and a list of free parameters then let it run for a few iterations.
So in order to make it work with my simulation, I only needed to compute the square difference of each simulated directional bin with my analytical lobe model and let BFGS find the parameter values for which this difference is minimal.
This is a classical minimization problem that has now gained tremendous momentum due to the blooming field of Machine Learning.
Anyway, minimization is the most timeconsuming part of the application though, as it sometimes takes a hundred iterations to fit a lobe, and that can take quite a while!
That is why I wrote an automation tool to let it work during the night...
Automation
I wrote an automation form that allows to configure a "simulation session":
You need to specify how many sampling directions you need and the parameters of the surface (i.e. its roughness, whether it's a metal, a dielectric or a diffuse material, its albedo or F_0 fresnel term) and let it run.
It's interactive in the sense you can see what it's doing in realtime, you can pause, stop, restart, redo one specific direction you are not satisfied with, etc.
The resulting lobe parameters for each simulation and each scattering order are dumped into a file for manual exploitation later.
Analytical Fit of Diffuse Lobes
Next, I had a lot of fun with Mathematica trying to find a suitable analytical expression that corresponds to my experimental lobes.
Here are the conclusions, straight from 2 years ago (january 2016), I'm not sure my results are accurate or useable anymore (I should have written this page when I was developping it, now it's hard to remember where I left it at):
Analytical Lobe Expression
The intensity of the lobe in a specific direction is given by:
Where:
 \sigma is the global scale factor
 m is the importance for the masking/shadowing term, which will be later set to 0 as we will see
 G(\mu,\alpha) is the masking/shadowing term for the Phong model which is actually that of Beckmann
 \mu_i = \boldsymbol{\omega_i} \cdot \boldsymbol{Z} is the cosine of the angle between the incoming direction and the macroscopic surface normal
 \mu_o = \boldsymbol{\omega_o} \cdot \boldsymbol{Z} is the cosine of the angle between the outgoing direction and the macroscopic surface normal
 \alpha is the surface roughness
And the BlinnPhong normal distribution:
with \eta(\alpha) = 2^{10(1\alpha)}1 defining the exponent based on the surface's roughness \alpha (notice the 1 in the end that allows use to have a 0 exponent to make constant lobes)
Warning
Don't expect the regular BlinnPhong model for microfacet models here: I wrote my own to fit my needs in this whole "fitting business"!
After fitting each parameter one after another, I noticed that:
 Incident light angle \theta has no effect on fitted lobe, assuming we ignore the backscattering that is visible at highly grazing angles and that would be better fitted using maybe a GGX lobe that features a nice backscatter property.
 Final masking importance m is 0 after all
 There is only a dependency on albedo \rho for the scale factor (that was expected) and it is proportional to \rho^2 for the 2^{nd} order, and to \rho^3 for the 3^{rd} order, which was also expected.
NOTE: We can safely assume there should be a \rho^N dependency for the Nth scattering order...
Parameters Fitting
Finally, we obtain the following generic analytical model of a rough diffuse surface for all scattering orders S > 1:
The exponent \eta\left(\alpha\right) is given as a function of surface roughness by:
The generic scale factor \sigma used for all scattering orders is given by:
Where: $$ \begin{align} a( \alpha ) &= 0.02881326115  0.92153748116 \alpha + 6.63272611438 \alpha^2  4.595702230 \alpha^3 \\ b( \alpha ) &= 0.09663259042 + 7.21414360220 \alpha  19.7868451171 \alpha^2 + 11.04205888 \alpha^3 \\ c( \alpha ) &= 0.10935692546  10.7904051575 \alpha + 28.5080366763 \alpha^2  15.66525827 \alpha^3 \\ d( \alpha ) &= 0.04376425480 + 5.24919600918 \alpha  13.5827073397 \alpha^2 + 7.348408854 \alpha^3 \\ \end{align} $$
The flattening factor \sigma_n along the main lobe direction Z is given by:
Where: $$ \begin{align} a(\alpha) &= 0.9136430  1.655480 \alpha + 1.39617 \alpha^2  0.320331 \alpha^3 \\ b(\alpha) &= 0.0447239 + 0.624740 \alpha \\ c(\alpha) &= 0.1188440  0.973213 \alpha + 0.36902 \alpha^2 \\ d(\alpha) &= 0.1325770 + 0.169750 \alpha \\ \end{align} $$
So the worldspace intensity of the fitted lobe is finally obtained by multiplying the lobespace intensity with the scale factor:
Scale factor for order 2
And the main takeaway here is the global scale factor for scattering order 2:
Scale factor for order 3
Identically, using the same generic parameters as order 2 and fitting the scale factor for order 3, we get:
General rule
Maybe there is a simple general rule to obtain the factor for any scattering order S, it would seem the general rule on \rho^S is quite clear but anyway, any scattering order above 3 is completely negligible so I basically stopped there...
Example Code
All this seems really complex but we eventually get the new code which ends up being "quite simple":
float3 ComputeDiffuseModel( float3 _wsIncomingDirection, float3 _wsOutgoingDirection, float _roughness, float3 _albedo ) { // Reorder components _wsIncomingDirection = float3( _wsIncomingDirection.x, _wsIncomingDirection.z, _wsIncomingDirection.y ); _wsOutgoingDirection = float3( _wsOutgoingDirection.x, _wsOutgoingDirection.z, _wsOutgoingDirection.y ); float cosTheta = saturate( _wsOutgoingDirection.z ); // Compute lobe scale, exponent and flattening factor based on incoming direction and roughness float mu = saturate( _wsIncomingDirection.z ); float mu2 = mu*mu; float mu3 = mu*mu2; float r = _roughness; float r2 = r*r; float r3 = r*r2; float4 abcd = float4( 0.028813261153483097  0.9215374811620882 * r + 6.632726114385572 * r2  4.5957022306534 * r3, 0.09663259042197028 + 7.214143602200921 * r  19.786845117100626 * r2 + 11.042058883797509 * r3, 0.10935692546815767  10.790405157520944 * r + 28.50803667636733 * r2  15.665258273262731 * r3, 0.04376425480146207 + 5.2491960091879 * r  13.582707339717146 * r2 + 7.348408854602616 * r3 ); float sigma2 = abcd.x + abcd.y * mu + abcd.z * mu2 + abcd.w * mu3; // 2nd order scattering float sigma3 = 0.363902052363025 * sigma2; // 3rd order scattering // Compute lobe exponent float eta = 2.588380909161985 * r  1.3549594389004276 * r2; // Compute unscaled lobe intensity float intensity = (eta+2) * pow( cosTheta, eta ) / PI; // Compute flattening factor abcd = float4( 0.8850557867448499  1.2109761138443194 * r + 0.22569832413951335 * r2 + 0.4498256199595464 * r3, 0.0856807009397115 + 0.5659031384072539 * r, 0.07707463071513312  1.384614678037336 * r + 0.8565888280926491 * r2, 0.010423083821992304 + 0.8525591060832015 * r  0.6844738691665317 * r2 ); float sigma_n = abcd.x + abcd.y * mu + abcd.z * mu2 + abcd.w * mu3; float L = rsqrt( 1.0 + cosTheta*cosTheta * (1.0 / (sigma_n * sigma_n)  1.0) ); // Add albedodependency return L * intensity * _albedo * _albedo * (sigma2 + _albedo * sigma3); }
The new code is used like this:
float3 diffuseTerm = (albedo / PI) * LdotN * shadow * lightColor; // Add multiplescattering term float shadowMS = ContrastShadow( shadow, LdotN ); // If not needed, just return "shadow" diffuseTerm += (ComputeDiffuseModel( _light, _view, roughness, albedo ) / PI) * shadowMS * lightColor;
The ContrastShadow()
should either return "shadow" in the basic case, or you can use the one I describe in the article about Color Shadows to give it a little coloring!
Result
You can see below the effect of multiplescattering on shadows and transition areas when the roughness increases:
This is a live demo of what's happening when we increase the roughness:
Conclusion
Regarding this lobe fitting business, you may understandably question the complexity of the computation of the multiplescattering term considering the low visual impact it's bringing to the table, and I would completely agree with you!
Unfortunately, I never had the time to finish this project due to the time constraint of working on the production of Dishonored 2 but I would have loved to continue experimenting, especially reusing Heitz's results instead of casting millions of rays, and fitting better lobe models or even find a much simpler way to add back the energy lost by singlescattering models.
Of course, people didn't stop investigate like I did, especially in large companies like Disney, Dreamworks, Weta or ImageWorks. And what had to happen did happen...
Energy Compensation
In their 2017 talk ^{5} at the now famous Siggraph's Physically Based Shading in Theory and Practice courses, Kulla and Conti reintroduced a computation ^{6} devised long ago by Kelemen and SzirmayKalos in 2001.
The Original Paper
In section 2.2 of their paper ^{6}, 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)} $$
Applications to existing BRDF models
We start by precomputing the "complement albedo table" for all possible viewing angle \theta_o and all roughness values \alpha for the specular BRDF:
GGX Specular Model
We use the (now) very common GGX normal distribution and Smith GGX shadowing/masking term:
With F_0 = 1 \Rightarrow F( \boldsymbol{\omega_o} \cdot \boldsymbol{\omega_h}, F_0 ) = 1 at the moment (i.e. perfectly reflective case).
\boldsymbol{\omega_h} is the normalized halfvector.
Rewritten in terms of \mu we have:
With \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) } } and \phi is the azimutal angle between \boldsymbol{\omega_o} and \boldsymbol{\omega_i}.
You can see the resulting table below:
Warning
Obviously, don't use this awful JPG image directly! Use 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!)
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)
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( \boldsymbol{\omega_h}, F_0 ) = 1, the tables we just calculated can only give us the perfectly reflective 100% white metal BRDF case:
Todo
TODO
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^{10}, then published again^{11} 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:
\sigma \in [0,\frac{\pi}{2}] is the standard angle deviation for the microfacets' slope distribution and 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 ComputeOrenNayar( in float3 _normal, in float3 _light, in float3 _view, in float _roughness ) { float3 n = _normal; float3 l = _light; float3 v = _view; float LdotN = dot( l, n ); float VdotN = dot( v, n ); // I realize that this doesn't give cosine phi, we need to divide by sqrt( 1VdotN*VdotN ) * sqrt( 1LdotN*LdotN ) // but I couldn't distinguish any difference from the actual formula so I just left that as it is... float gamma = dot( v  n * VdotN, l  n * 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 saturate( LdotN ) * (A + B * max( 0.0, gamma ) * C); }
Irradiance Table
Once again, 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 16x16 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.
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.
With varying Fresnel
References

Heitz, E. 2014 "Understanding the MaskingShadowing Function in Microfacetsbased BRDFs " ↩

Heitz, E. Dupuy, J. Crassin, C. Dachsbacher, C. 2014 "The SGGX Microflake Distribution" ↩

Heitz, E. d'Eon, E. 2014 "Importance Sampling MicrofacetBased BSDFs using the Distribution of Visible Normals" ↩

Heitz, E. Hanika, J. d'Eon, E. Dachsbacher, C. 2015 "MultipleScattering Microfacet BSDFs with the Smith Model" ↩↩↩

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" ↩↩

d'Eon, E. 2016 "A Hitchhiker’s Guide to Multiple Scattering" ↩

Walter, B. Marschner, S. R. Li, H. Torrance, K. E. 2007 "Microfacet Models for Refraction through Rough Surfaces" ↩

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

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" ↩