A Multi-Faceted Exploration (Part 4)

$$ \newcommand{\fss }{\color{gray}{f_\mathrm{ss}}} \newcommand{\fssp}{\color{gray}{f_\mathrm{ss}'}} \newcommand{\fssoi }{\color{gray}{f_\mathrm{ss}(\mu_o, \mu_i)}} \newcommand{\fms }{\color{purple}{f_\mathrm{ms}}} \newcommand{\E }{\color{olive}{E}} \newcommand{\Emo}{\color{olive}{E(\mu_o)}} \newcommand{\Emi}{\color{olive}{E(\mu_i)}} \newcommand{\Em }{\color{olive}{E(\mu)}} \newcommand{\Eavg}{\color{green}{E_\mathrm{avg}}} \newcommand{\Favg}{\color{teal}{F_\mathrm{avg}}} \newcommand{\Fms }{\color{brown}{F_\mathrm{ms}}} \newcommand{\Ess }{\color{olive}{E_\mathrm{Fss}}} \newcommand{\Esso}{\color{olive}{E_\mathrm{Fss}(\mu_o)}} \newcommand{\Ems }{\color{teal}{E_\mathrm{Fms}}} \newcommand{\Emso}{\color{teal}{E_\mathrm{Fms}(\mu_o)}} \newcommand{\Fr }{\color{purple}{F_0}} \newcommand{\Fs }{\color{brown}{F_\mathrm{schlick}}} \newcommand{\Fsr }{\color{brown}{F_\mathrm{schlick}(\omega_h, \omega_r)}} \newcommand{\ki }{\color{}{k}} \newcommand{\wa }{\color{brown}{w_0}} \newcommand{\wb }{\color{brown}{w_1}} \newcommand{\wc }{\color{brown}{w_6}} \newcommand{\wi }{\color{brown}{w_i}} \newcommand{\wn }{\color{brown}{w_n}} \newcommand{\wnm }{\color{brown}{w_{N - 1}}} \newcommand{\ea }{\color{teal}{e_0}} \newcommand{\eb }{\color{teal}{e_1}} \newcommand{\ec }{\color{teal}{e_2}} \newcommand{\ei }{\color{teal}{e_i}} \newcommand{\eim }{\color{teal}{e_{i - 1}}} \newcommand{\en }{\color{teal}{e_n}} \newcommand{\enm }{\color{teal}{e_{N - 1}}} \newcommand{\ss }{\color{green}{s}} \newcommand{\si }{\color{green}{(1 - s)}} \newcommand{\ssa }{\color{green}{s_{1}}} \newcommand{\sia }{\color{green}{(1 - s_{1})}} \newcommand{\ssb }{\color{olive}{s_2}} \newcommand{\sib }{\color{olive}{(1 - s_{2})}} $$

Bouncing Back

Last time around, I showed that we could improve upon Imageworks’ multiple-scattering approximation by precalculating a new term, $\Ems$, from the Heitz model itself. This is, of course, a well-trodden path in rendering: if you can’t use something directly then try to tabulate or fit it instead!

To make this more practical, I also outlined a multiple-scattering variant of the split integral trick already in common use in real-time rendering, which allows the specular colour, $\Fr$, to be factored out.

I’ll now go into the details of this precomputation process, and also show how, with further crunching, we can reduce the storage and run-time cost even more.

Walking the Walk

A practical approach for precomputing the directional albedo of a single-scattering BRDF is via Monte Carlo integration with BRDF importance sampling.¹ We can calculate $\Ems$ from the Heitz model in a similar way via its random-walk sampling process.

Conceptually, each random walk starts with a ray coming from the camera, which bounces one or more times on the microsurface before finally escaping it:

Figure 1: An illustration of the random-walk sampling process. (Courtesy of Eric Heitz.)

At each bounce, some energy is lost due to absorption and the rest is reflected.² In the case of our copper conductor test subject, the final energy (or energy throughput) is simply the product of Fresnel reflection at each bounce.

If we perform this process many times and average the final energy of every multi-bounce walk – i.e., treating single-scattering walks as having zero energy – then we obtain $\Ems$, as shown in the last post:

*Figure 2: $\Ems$ for copper GGX material.*

Doing the Splits

As I previously discussed, we don’t want to precompute $\Ems$ directly, since it bakes in $\Fr$. Instead, we’ll precalculate factors, $\wa \ldots \wnm$, which can be used to reconstruct $\Ems$ at run time for any $\Fr$:

$$ \begin{equation} \Emso = \sum_{i=0}^{N-1} \wi \Fr^i. \label{eq:ms_sum} \end{equation} $$

To achieve this, we’re going to make two very simple modifications the random walk process:

The walk will now keep track of a vector of energy throughputs, $\ea \ldots \enm$, rather than a single value. At the beginning of the walk, $\ea$, is initialised to 1 and $\eb \ldots \enm$ to 0.
At each scattering event, the energy throughput was previously multiplied by a Fresnel factor, $F(\omega_m, \omega_r)$, where $\omega_h$ is a sampled normal and $\omega_r$ is the outgoing ray direction. If we use the Schlick Fresnel approximation for $F$, this is

$$ \begin{align} \Fsr &= \Fr + (1 - \Fr)(1 - \omega_h \cdot \omega_r)^5 \\ &= \Fr \, \si + \ss, \quad \text{where } \ss = (1 - \omega_h \cdot \omega_r)^5.\label{eq:fres} \end{align} $$

$\quad$ We will now use $\ei$ to track the fraction of energy scaled by $\Fr^i$.

From Equation $\ref{eq:fres}$, we can see that that after the first bounce, a factor $\sia$ of the energy is multipled by $\Fr$ and the rest, $\ssa$, is left untinted, so we set $\ea = \ssa$ and $\eb = \sia$. On the next bounce, with a new Fresnel factor $\ssb$, we will have $\ea = \ssa \ssb$, $\eb = \sia \ssb + \ssa \sib$ and $\ec = \sia \sib$. More energy has moved from order 0 to order 1, and some from order 1 to order 2.

Working things though, the update to our energy throughput vector at each bounce can be written as

$$ \begin{align} \color{teal}{e'_0} &= \ss \, \ea \\ %\color{teal}{e'_1} &= \ss \, \eb + \si \, \ea \\ %\ldots \\ \color{teal}{e'_i} &= \ss \, \ei + \si \, \eim. \end{align} $$

Or in code:

1for (int i = N - 1; i > 0; i--)
2    e[i] = mix(e[i - 1], e[i], s);
3e[0] *= s;

The final $\wa \ldots \wnm$ factors that we’re after are simply the averages of $\ea \ldots \enm$ over multiple runs of the random walk process (again, ignoring single scattering):

*Figure 3: Multiple-scattering Fresnel factors for GGX.*

You can see the whole process in action in this WebGL demo, which reproduces Figure 2 as well as the decomposition shown in Figure 3.

Note: there’s quite a bit going on in this demo, but I didn’t want to get into the weeds in this post, so I’ll write up some separate notes at a later date.

Losing Weights

Although it may be hard to tell from Figure 3, there is some residual energy in the bottom right corner³ of $\wc$, and this is true for a few higher orders as well. This implies that we’d need two or three textures to store all of the factors for fully accurate results, as well as a fair amount of maths in the shader. It would be nice to slim this down!

Fortunately, since most of the energy is in the lower orders, the source polynomials for $\Ems$ can be accurately refitted to lower-order cubic curves. For instance, here’s a plot of $\Ems$ for the bottom right corner, together with the refit:

This means that, in practice, we only need a single four-channel texture and three MAD instructions (when Equation $\ref{eq:ms_sum}$ is written in Horner form) to reconstruct $\Ems$ in the shader:

1vec3 EFms(vec3 F0, float mu, float rough)
2{
3    vec4 w = EFmsFactors(mu, rough);
4    return w.x + F0*(w.y + F0*(w.z + w.w*F0));
5}

I did this refitting in Mathematica using MiniMaxApproximation to minimise relative error, which proved important for accurately reproducing the lower end of the curves.

You can see the end result for yourself in this second WebGL demo, which compares Imageworks along with the improvements discussed so far (left half), against the Heitz reference (right half).

Note: I realise that not everyone has ready access to Mathematica, so I’m planning to rewrite the refitting step in either Python or C++ before publishing the complete end-to-end process on GitHub.

Bonus (Reading) Exercise

Up until now, I’ve been probing Imageworks’ energy-preservation solution, which approximates microfacet multiple scattering with a diffuse-like lobe. However, there are other potentially valid options here. For instance, concurrent to Imageworks’ investigations, Emmanuel Turquin developed a different solution while at Industrial Light & Magic, which instead approximates the multiple scattering via a rescaled single scattering lobe. This approach has since been adopted within Unity’s High Definition Render Pipeline⁴ as well Google’s Filament.

With permission, I have the privilege of hosting a technical report from Manu that goes into all of the details of his approach and various shortcuts, along with some comparisons to Imageworks and Heitz. That said, in conclusion he also states:

“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.”

In the next post, I hope to go some way towards providing this more detailed analysis. In the meantime, you should read Manu’s TR!

For an example of this, see the Image-Based Lighting section of [Karis 2013]. ↩︎
Since the focus is still on conductors, I’m ignoring refraction and transmission for now. ↩︎
Corresponding to low view angles and high roughness values. ↩︎
See the Multiple Scattering GGX section of [Lagarde and Golubev 2018]. ↩︎