Practical multiple scattering compensation for microfacet models
Emmanuel Turquin∗
Industrial Light & Magic

Figure 1: Increasingly rough GGX conductor widgets. Top: single-scattering only. Bottom: with our multiple scattering compensation.
Notice the boosted intensity and saturation as the surface becomes rougher, which offers a visually more satisfying and uniform gradation.

Abstract

A particularly popular choice is the use of the GGX distribution as
D (and its associated term G; see [Heitz 2014] on how they should
relate), along with a conductor Fresnel term, to reproduce roughmetals (as in Figure 1, top row). Alternatively, adding a BTDF
combined with a dieletric F can replicate refractive materials suchas frosted glass.

With the adoption by the game, animation and VFX industries ofmore physically grounded light transport methods and material representations, an effort has been put on energy conservation, to ensure both numerical stability of the computations and realistic looking images. However, this effort has been much more focused on
avoiding undue light emission than ensuring no energy is lost in
the process. Microfacet models have become a standard building
block of surface materials for representing specular components
of varying roughness; and yet, while they possess many desirable
properties in addition to producing convincing results, their very
design neglects an important source of scattering, that can causea significant loss of energy. Specifically, they only model singlescattering across the microfacets, and ignore the subsequent interactions, which get more important as roughness increases. Froma user standpoint, this results in an unexpected darkening of roughspecular lobes, which typically has to be accounted for in ad hocways. In this document, we present and compare different approaches aiming to address this flaw and ensure energy preservation, including one that has been developed at ILM.

Figure 2: White furnace test on a microfacet BRDF with F = 1,
and roughness going from 0 (left) to 1 (right). From top to bottom:
Beckmann, GGX, GTRγ=1 and STDγ=1.55 distributions. If multiple
scattering was properly accounted for, the albedo would be 1.

Keywords: microfacets, multiple scattering, energy compensation

1

Unfortunately, this model also suffers from an important limitation:
it only simulates a single interaction of the incoming ray of lighton the microfacets, after which the scattered ray is either visible or
not from ωi /ωo , as determined by G. The multiple scattering thatwould occur when the ray is not visible (i.e., hits another microfacet, once or more) is not accounted for (see Section 3 of [Heitz
2014] for more details). As one would expect, there is no missing
energy when the roughness is zero, as this corresponds to a perfectly flat surface with no possible occlusion. But the rougher thesurface (or the wider the spread of microfacet orientations is), themore this lack of multiple scattering energy becomes noticeable.
Looking back at Figure 1, top row, the right-most objects appear
overly dark and dull, in comparison to the shinier ones. As a matterof fact, in the case of GGX and when ignoring Fresnel absorption,
the loss of energy reaches about 60% for a roughness value α = 1.
It gets even worst when using distributions with longer tails, suchas GTR [Burley 2012] or STD [Ribardière et al. 2017], where it canreach in excess of 90%, as illustrated by Figure 2. In a production
environment, this calls for an eye-balled, manual compensation of
the albedo, applied typically at the look development stage, and potentially breaking energy conservation in other ways. Thankfully,
easier and sounder methods have been devised.

Introduction

Since its introduction by Cook and Torrance [1982], the microfacet
model has become an ubiquitous way of repesenting rough specularreflections and refractions [Walter et al. 2007], both in real-time andoffline rendering. If we consider mirror-like reflecting facets, a microfacet BRDF ρ can be expressed using the well known formula:. F (ωo , h) G (ωo , ωi , h) D (h)ρ (ωo , ωi ) =4|ωi · n||ωo · n|(1)

with n being the geometric normal, h the microfacet normal, ωiand ωo the incoming and outgoing directions, F the Fresnel term,
G the shadowing-masking function and D the microfacet distribution function. This formulation, amongst other advantages, offersa great modularity as all three terms F , G and D can be chosen
between candidates with different characteristics.
∗ e-mail: eturquin@ilm.com

1

2

Previous work

“If it turns out to be simple, then as a first approximation
we could model it with an analytical function (e.g. as asingle lobe).” (Q2)

The various ways of reinjecting the missing multiple scattering contribution range from heuristic and artistically controlled, to fully
physically based and automated. They also vastly differ in terms ofcomplexity, practicality and speed.

About Fresnel, he suggests to precompute an average term Fss for
one bounce and rescale Ems by it. If deemed not accurate enough,“perhaps the average value Fms after multiple bounces could beprecomputed as well”. Concluding with this promising observation:

Burley [2015] introduces a simple, physically justified, user-driven
Sheen component using the Schlick Fresnel profile (1 − cos θo )5 ,that adds grazing forward reflection to better match observed materials, and complements the diffuse retro-reflection already proposed in [Burley 2012]. The motivation for this term is to approximately compensate for the assumed missing “multiple-scattering
effects between and through microsurface features”.

“Since multiple scattering tends to smooth out functions,
one can reasonably expect it to be efficiently represented
and stored with simple analytical functions or small precomputed look-up textures.” (Q3)
In a subsequent effort, Heitz et al. [2016] follow up on the idea
expressed in Q1, and propose a full-on stochastic model to simulate multiple scattering, using the microflake theory for volumes
and random walks on the microsurface. The model is shown to be
extremely accurate, not only for conductor reflections but also dielectric scattering, and would be ideal if not for its involved naturethat makes it not so straightforward to plug into an existing production rendering system (notably with the use of additional random
numbers), and more importantly computationally very expensive.

Figure 3: Effect of Sheen, ranging from 0 to 1, as in [Burley 2015].
Unfortunately, this qualitative compensation (illustrated in Figure 3) fails to take into account important factors, such as the effectof roughness variation, and thus needs to be readjusted wheneverother attributes of the BSDF are edited. In Section 5.1 of thosesame course notes, an emphasis is put on how much of a practical issue it still is in production, while a more exact and automatic
method is mentioned as future work.

An interesting and somewhat surprising take from this work (visible
in Figure 15 of their paper, and partially reproduced here in Figure 4
for convenience), is that secondary lobes do not look diffuse but
rather like scaled-down versions of the primary one.

In the Multiple Scattering sub-section of Section 7 of [Heitz 2014],
a suggestion is made to “combine the knowledge of energy conservation and empirical observations”, and to consider a new BRDF
model that would be expressed as:.
ρ (ωo , ωi ) = ρss (ωo , ωi ) + ρms (ωo , ωi )(2)
where ρss is the usual, single-scattering only term of Equation 1,
and ρms a new multiple scattering term to be defined, that accounts
for all secondary lobes. A first important constraint on ρms is naturally imposed by energy preservation. Let us introduce the albedo
as:Z.
ρ (ωo , ωi ) |ωi · n|dωi(3)
E (ωo ) =Ωi(similarly, Ess and Ems ). Then, if we momentarily ignore Fresnel(i.e. we consider that F = 1), the following equality must hold:E (ωo ) = Ess (ωo ) + Ems (ωo ) = 1

Figure 4: Resulting lobes from a simulated microsurface using the
Beckmann distribution, courtesy of Heitz et al. [2016].

(4)

expressing the fact that without any absorption, all incoming energy
is reflected back, after one or more bounce(s) on the microfacets.
Put differently:
Ems (ωo ) = 1 − Ess (ωo )In an effort concurrent to ours, Kulla and Conty [2017] exploit observations made in Q2 and Q3 and propose a significantly less accurate, but much simpler and faster solution. They adapt the work of
Kelemen and Szirmay-Kalos [2001] to the present problem, and use
their diffuse-looking matte component as multiple scattering lobe,
which by construction respects Equation 5:(5)

Given that Ess is completely defined, we have different options to(pre-)compute it, or Ems . Respecting Equation 5 prevents energy
loss, and would naturally lead to a perfectly constant Figure 2, as
desired. Now, what remains to be explored is:(1 − Ess (ωo ))(1 − Ess (ωi )).
ρms (ωo , ωi ) = Fmsπ(1 − Eavg )• The shape of ρms(6)with the Fresnel term defined as:• How to incorporate Fresnel absorption.

.
Fms =

While Heitz [2014] doesn’t offer a direct solution to either problem,
he notes that:Fss Eavg
1 − Fss (1 − Eavg )(7)

R1
and with Fss = 2 0 F (µ) µdµ computed assuming diffuse reR1
flection and Eavg = 2 0 Ess (µ) µdµ. Note that they choose to
remap F and Ess with µ = |ω · n| instead of ω, which is perfectly
fine as both terms are azimuthally invariant. The detailed derivation of the Fms term, obtained using a geometric series, is given“The shape of ρms could be investigated, for instance,by computing Monte-Carlo simulations on rough surface samples.” (Q1)adding right after:2

Beckmannin [Jakob et al. 2014], Section 5.6. The overall method relies on
small 1D (α) and 2D (α, µ) look-up tables for the smoothly varying Eavg and Ess (with dimensions of 32 and 32×32, respectively,in their implementation), which have to be precomputed for eachgiven microfacet distribution. By design the lobe is also reciprocal(i.e. ρms (ωo , ωi ) = ρms (ωi , ωo )), which is crucial if used in abi-directional setting.

GTRγ=1STDγ=1.55

Figure 5: Ess for different microfacet models,quantized at 32×32.
√
Left to right: cos θo , and top to bottom: α ∈ [0, 1].

As an additional contribution, given that F also depends on its IOR(two wavelength-dependent values for conductors), they offer simple analytical fits of Fss for various Fresnel formulations (Schlick
approximation, conductor, artist-friendly conductor). They also
support rough dielectric reflection and refraction, this time using
a couple of tables, one for each type of interface (exterior towards
interior, and the other way around). In this case, Fresnel has to
be incorporated into the energy LUTs, as it now represents a ratio
between reflected and refracted rays, thus adding a dimension totables of otherwise matching resolution.

3.1.2

Fresnel term

Before going any further, we have to acknowledge that the decision
to model multiple scattering with a single lobe, and to consider Fms
as separable and directionally invariant, is a very coarse approximation. Within this design, we can only expect a partial match at best.
However, a simple qualitative match will likely be good enoughfor our production needs, and for a user it makes the behaviour ofthe compensation more easily predictable, with a saturation directly
proportional to the amount of energy being added.

All in all, this method is very fast and robust, and can be added toa renderer with relatively minor modifications to its BSDFs. However, it elects to use a diffuse-looking ρms lobe, hence contradicting previous results from Heitz et al. [2016], which showed that thesecondary lobes are rather similar in shape to the primary one.

3

GGX

Our method

Whereas Kulla and Conty’s main concern and design constraint
was to preserve reciprocity, our method stems from this observed
similarity of shapes between lobes, while also aiming at extreme
minimalism (unless more complexity would prove beneficial to theuser). This brought us to simply try and use a scaled ρss as ρms :.
ρ (ωo , ωi ) = ρss (ωo , ωi ) + Fms kms (ωo ) ρss (ωo , ωi )(8)

Figure 6: Comparison between no mult. scattering compensation,compensation without and with Fms , on a rough GGX conductor.
Note the difference in color saturation when Fresnel is applied.

with kms a factor to be defined, accounting for the missing energy.

A great option for Fms is to simply reuse the term defined in Eq. 7.
This function goes from being exactly equal to Fss when α = 0(or Eavg = 1), corresponding indeed to a case where there is only2
single-scattering, to being close to Fss
for α = 1 (or Eavg ≃ 0.4).

In this section, we use a GGX BRDF, unless mentioned otherwise.

3.1

Conductors

3.1.1

Energy term

2
Looking at Figure 7 left1 , we see that Fsswould give slightly moresaturation than even this formulation, which happens to be visuallya bit closer to ground truth results (more on that in Section 4).

If we push aside the effect of Fresnel on the multiple scattering
term, we can actually see ρ as a normalized version of ρss :ρss (ωo , ωi )
ρ (ωo , ωi ) = (1 + kms (ωo )) ρss (ωo , ωi ) =Ess (ωo )We use it with a minor alteration though, as ρss already incorporates a directional single-scattering Fresnel coefficient F (seeEquation 1), effectively replacing Fss at the numerator. We canjustify this as a way to give more directionality to our overall effect, counter-balancing the diffuse assumption made in the original
derivation.
Eavg
.
Fms =
(12)
1 − Fss (1 − Eavg )(9)Naturally giving us:
kms (ωo ) =1 − Ess (ωo )Ess (ωo )And by construction, such a BRDF verifies Equation 4:Z
Ess (ωo )ρss (ωo , ωi )|ωi · n|dωi ==1E (ωo ) =
E
(ω
)
E
ss
o
ss (ωo )Ωi(10)

The function’s behaviour corroborates measurements made with the
Heitz [2016] model, showing that the average random walk depth(which can be seen as an exponent for Fss ) varies from 1 to almost
2, when α goes from 0 to 1 (see Figure 7 right).

(11)

Due to its smooth variation, we found like Kulla and Conty that
Ess can be precomputed (again, with F = 1) and stored in a small
LUT, parameterized by outgoing angle (or cosine) and roughness.
In our implementation, we too use 32 × 32 tables with cos θo and√α (see Figure 5). For distributions with an additional γ tail-lengthparameter, such as GTR or STD, we extend the tables to 32 × 32 ×
32 with a γ-mapping tailored to each distribution.

Exploiting these results, another (even simpler) suitable formulation for Fms would be:. α√αFms = Fss(13)

We can simplify the Fresnel term further still, by observing that, asa multiplier of kms in Eq. 8, it makes more of an impact when the
kms value is high, which corresponds to large roughness values.

But while the energy term is probably the most critical componentof the overall compensation, the Fresnel absorption happening at
each bounce on a microfacet is far from negligible, as shown in
Figure 6. We still need to find a satisfying Fms term.

1 A dynamic version of these curves is available online at [Turquin 2017].

3

S
Given that Fresnel is included in the Essterm, our LUT becomes
3D, with the IOR as an additional dimension. On the other hand,
this means that we do not need to worry about a separate multiple
scattering Fresnel term anymore. As Kulla and Conty [2017], we
compute tables separately for the two types of interface, leading
to two 32 × 32 × 32 tables (note that the tables would even be
4D if computed for a tail-length parametrized model), and use thebounded F0 ∈ [0, 1], instead of η for the extra dimension.

Figure 7: Left: Fms as defined in Eq. 7, for Eavg = 0.4 (i.e.
2α ≃ 1), as a function of Fss . We can see that it gets close to Fss.
Right: In blue: Measured average random walk depth
for
GGX,
as
√
a function of α. In red: fitting the data with 1 + α α.
Since the influence of Fms is most visible when α = 1, we canstrip it down to the value it takes for such a roughness:.
Fms = Fss
(14)
We are not completely done yet: a last simplification can be applied. Fss , as defined in [Kulla and Conty 2017], is expressed asa polynomial fit with wavelength-dependent parameters η and κ,or alternatively F0 and Fedge (F0 being the color at normal angle;
see [Gulbrandsen 2014] for details on this parameterization). Its
evaluation has a cost, especially if targetting real-time applications.

Figure 8: Comparison between rough (α = 0.5) GGX glass object
without (left) and with (right) multiple scattering compensation.
As can be seen in Figure 8, the effect of this compensation is quite
noticeable – maybe even more so than with conductors.

In our tests, for strongly tinted, physically plausible conductors
such as gold or copper, Fss and F0 values are so close that they
can be interchangeably used with no visual difference in the results. Note that this observation only holds when using a physicalFresnel; if we pick a more artist-driven one where grazing values
can be arbitrarily colored, it is recommended to stick with Eq. 14.
This leads us to the ultimate, bare-bones formula:.
Fms = F04

The method presented in Section 3 has been implemented in
Isotropix Clarisse, Pixar RenderMan’s RIS and the Mitsuba pathtracer, but could easily be added to any other renderer (including
real-time rasterization engines).
It is close to Kulla and Conty’s, essentially with a different choiceof shape for ρms : the same lobe as ρss in our case, as opposed toa diffuse one in theirs. The former is closer to the ground truth (asshown in [Heitz et al. 2016]), but not reciprocal, which could bea deal breaker when the resulting BSDF is used with bidirectional
light transport. However, we believe that this lack of formal reciprocity is unlikely to produce any significant artefacts.

(15)

and if we rewrite our BRDF defined in Eq. 8 more explicitly:`
1 − Ess (ωo ) ´ρss (ωo , ωi )
ρ (ωo , ωi ) = 1 + F0
Ess (ωo )(16)This is the formulation that has been used to generate the results of
Figure 1, bottom row. Interestingly, no matter what Fresnel term
is chosen between 12, 13, 14 and 15, the results are visually very
close, and they all achieve the saturated effect we were originallylooking for. It is likely that anything close to an F 2 shape towards
high roughnesses would do the trick for our production needs.

3.2

Because our method depends only on ωo , rather than (ωo , ωi ), it iseven simpler to integrate into an existing renderer, as the compensation can be applied as a gain to the closure, instead of having tomodify the closure itself. Being simpler, it is marginally faster too.
Another beneficial side effect is that, as our result can be seen asa scaling of ρss , it perfectly maintains its shape, meaning we canreuse the exact same sampling and evaluation routines, with thesame PDF. In contrast, Kulla and Conty’s lobe is a mixture of ρss
and an azimuthally invariant lobe, for which optimal analytical importance sampling cannot be achieved.

Dielectrics

For dielectrics, the approach is similar with few key differences. We
cannot neglect Fresnel when determining albedo normalization theway we did with conductors in Eq. 9, as it is indicative of the ratioRT
between reflected Ess
and transmitted Ess
energies, and not justreflected vs absorbed anymore. What can be normalized, however,S
R
T
is the sum Ess
= Ess
+ Ess
, each time an interface is hit.
If we define our single-scattering dielectric BSDF as:. R
TρSss (ωo , ωi ) = ρss (ωo , ωi ) + ρss (ωo , ωi )Both methods, while theoretically less elegant or accurate, and far
from representing the same step forward scientifically, are markedly
easier to implement than Heitz et al. [2016], and much faster at runtime. In tests performed in Mitsuba using their original implementation, a 7× slowdown was observed on a simple conductor test,
and as high as 15× on a dielectric transmission test, while the visual difference appeared minimal (see results in Figure 9).

(17)

and consider that the same reflection/transmission ratio applies tomultiple scattering, we can obtain our energy preserving BSDF:
R

The most notable difference is the amount of saturation on conductors, with the reference looking slightly more saturated at highroughnesses, even with us using Fms of Eq. 15, rather than the morejustified one of Eq. 12, that saturates the result a bit less (Fig. 7 left).

T. ρ (ωo , ωi ) + ρss (ωo , ωi )ρS (ωo , ωi ) = ss
S (ω )EssoComparison with other approaches

(18)

4

Method
Heitz et al.
Kulla and Conty
Ours

Shape of secondary lobes
Similar to primary
Azimuthally invariant
Same as primary

Reciprocity
Yes
Yes
No

Optimal sampling
No
No
Yes

Plausibility
A+
B
B+

Simplicity
C
A
A+

Speed
C
A
A+

Flexibility
B
A+A+

Table 1: Qualitative summary of the pros and cons of each method.

Acknowledgements

In any case, we cannot expect a perfect match, given that Fms
should not be separated from ρms as we did, and would definitely
vary directionally. We believe, however, that the results are closeenough to warrant keeping the approach as simple as it is.

Big thanks are in order to the ILM Rendering R&D team, Stephen
Hill and Naty Hoffman from Lucasfilm, Brent Burley and Matt
Chiang from Disney Animation, and Chris Kulla from Sony
Imageworks, for insightful discussions and feedback that greatly
helped shape this document and figure out its trickiest aspects.

References
B URLEY, B., 2012. Physically based shading at Disney. ACM SIGGRAPH course notes, http://blog.selfshadow.com/
publications/s2012-shading-course/.
B URLEY, B., 2015.
Extending the Disney BRDF to a
BSDF with integrated subsurface scattering.
ACM SIGGRAPH course notes, http://blog.selfshadow.com/
publications/s2015-shading-course/.
C OOK , R. L., AND T ORRANCE , K. E. 1982. A reflectance model
for computer graphics. ACM Trans. Graph. 1, 1 (Jan.), 7–24.

Figure 9: Comparisons for rough GGX reflection and transmission.
Left to right: no compensation, ours and [Heitz et al. 2016].
Top to bottom: medium, high and max roughness.

G ULBRANDSEN , O. 2014. Artist Friendly Metallic Fresnel. Journal of Computer Graphics Techniques (JCGT) 3, 4 (December),
64–72.

Finally, both [Kulla and Conty 2017] and our method can accommodate any sort of BSDF, as they only rely on the computation of
the single-scattering albedo Ess , whereas [Heitz et al. 2016] imposes considerably stronger constraints, making it even more difficult to use for models other than Beckmann and GGX.

H EITZ , E., H ANIKA , J., D ’E ON , E., AND DACHSBACHER , C.
2016. Multiple-scattering microfacet BSDFs with the Smith
model. ACM Trans. Graph. 35, 4 (July), 58:1–58:14.

For a quick summary of all pros and cons, please refer to Table 1.

H EITZ , E. 2014. Understanding the masking-shadowing function in microfacet-based BRDFs. Journal of Computer Graphics
Techniques (JCGT) 3, 2 (June), 48–107.

5

Conclusion and future work

JAKOB , W., D’E ON , E., JAKOB , O., AND M ARSCHNER , S.,
2014. A comprehensive framework for rendering layered materials. Expanded technical report, http://www.cs.cornell.
edu/projects/layered-sg14/.

We have presented different approaches to compensate for the missing multiple scattering energy in rough reflections and transmissions. For a production use, the most suitable candidates appear
to be either [Kulla and Conty 2017] or the one introduced here,
depending on how critical reciprocity is for the renderer it is integrated into. We are convinced that either one will be of great help
to artists, and remove the need for more error-prone, ad hoc fixes.

K ELEMEN , C., AND S ZIRMAY- KALOS , L., 2001. A Microfacet
based coupled specular-matte BRDF model with importancesampling. Short presentation, Eurographics.
K ULLA , C., AND C ONTY, A., 2017. Revisiting Physically Based
Shading at Imageworks. ACM SIGGRAPH course notes,
http://blog.selfshadow.com/publications/s2017-shading-course/.

Comparisons in this report have been largely of a qualitative nature. A more thorough and quantitative analysis of the differences
between the three main methods is left for future work.

R IBARDI ÈRE , M., B RINGIER , B., M ENEVEAUX , D., AND S I MONOT, L. 2017. STD: Student’s t-Distribution of slopes for
microfacet based BSDFs. Computer Graphics Forum.

Furthermore, it would be interesting to find good analytical fits forall of the smooth look-up tables, which might be beneficial to a realtime implementation, by virtue of avoiding costly texture accesses.

T URQUIN , E., 2017. Multiple scattering Fresnel curves. https://www.desmos.com/calculator/0mig1qjjvf.

We could also consider not separating Fresnel for conductors anymore, and instead extend the LUT to 3D (or 4D if we really wantto take the imaginary part of the IOR into account), computing itwith a reference algorithm such as [Heitz et al. 2016] (at least for
the Beckmann and GGX distributions).

WALTER , B., M ARSCHNER , S. R., L I , H., AND T ORRANCE ,K. E. 2007. Microfacet models for refraction through roughsurfaces. In Proceedings of the 18th Eurographics Conference
on Rendering Techniques, EGSR’07, Eurographics, 195–206.

Speaking of which, as a much more ambitious future endeavour,
it would be fantastic to figure out a way to turn their method into
something more production friendly – the critical element being
speed here, with ease of use coming second.

5