Brute Force Processing
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. ^{1} 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 ^{3}
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. ^{2}

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 Colored 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...
References

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

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