Fog rendering has a very rich history in video game graphics. Early games relied on the fixed pipeline with its built-in fog functions. Back in those days graphics cards were simply accelerators of predetermined set of graphics processing functions. Some people might even remember comparisons of linear and exponential fogs. Of course, both don’t hold much physical meaning. Early games, like Turok as many remember, even overused the effect to hide the very limited draw distance. You might see that effect referred to as “distance fog”. With the advent of the programmable pipeline many more types of fog came into use. Another kind of fog that came into use is sometimes called “height fog”. We speak a lot about physically based rendering in video game graphics, but somehow forget what is the actual phenomena leading to those effects. Of course none of those physically exist. They are just approximations of very specific atmospheric conditions. In reality, what we have is a participating medium. A volume where light bounces around many times until it reaches the eye. Many physical phenomena can be successfully modeled in that framework: fog, clouds, atmosphere, smokes (particle systems), liquids, dielectric solids, skin, cloth. But here I will talk only about fog and atmosphere because both can be modeled quite successfully with analytic media with a caveat that atmosphere is in practice almost a spherical shell around the Earth surface, so flat approximations need corrections for that fact. Before we step into the depths of modeling fogs, I will make a short introduction of some concepts for completeness.

The more general equation in rendering

Many people in graphics when asked what is the most important equation will often cite The Rendering Equation by Kajiya from 1986 [4]. Will that equation allow us to model fogs and smokes? Unfortunately, no. It is a boundary condition of a much more complex equation. We have to go back in time to discover a slightly better appoximation that balances the radiance passing through a cross section. Chandrasekhar [3] defined it as a differential equation, but lets just express it as simple terms,

$$\begin{eqnarray} (\omega \cdot \nabla) L(\mathbf{x}, \omega) &= -&\mathrm{Absoption} -\mathrm{Out\mbox{-}scattering} + \\\
&& \mathrm{In\mbox{-}scattering} +\mathrm{Emission}. \end{eqnarray}$$

So given the direction \(\omega\) of light flowing through the cross section, certain amount of light is absorbed, other scatters in and out, and emitted in the direction. Well, let’s see it in its mathematical glory for completeness,

$$\begin{eqnarray} (\omega \cdot \nabla) L(\mathbf{x}, \omega) &= -&\sigma_a(\mathbf{x}) L(\mathbf{x}, \omega) - \sigma_s(\mathbf{x}) L(\mathbf{x}, \omega) + \\\
&& \sigma_s(\mathbf{x}) L_s(\mathbf{x}, \omega) + \sigma_a(\mathbf{x}) L_e(\mathbf{x}, \omega). \end{eqnarray}$$

So we have absorption governed by the absorption coefficient and \(\sigma_a\) and in- and out-scattering by the scattering coefficient \(\sigma_s\), but we can slightly shorten that equation using the extinction coefficient \(\sigma_t = \sigma_a + \sigma_s\),

$$ (\omega \cdot \nabla) L(\mathbf{x}, \omega) = -\sigma_t(\mathbf{x}) L(\mathbf{x}, \omega) + \sigma_s(\mathbf{x}) L_s(\mathbf{x}, \omega) + \sigma_a(\mathbf{x}) L_e(\mathbf{x}, \omega). $$

By integrating both sides along a ray, we end up with the equation that is solved by path tracers in graphics engines,

$$\begin{eqnarray} L(\mathbf{x}, \omega) = \int_0^\infty \mathrm{d}y \, e^{-\int_0^y \mathrm{d}s \, \sigma_t(\mathbf{x} - s\omega)} (&\sigma_s(\mathbf{x}) L_s(\mathbf{x} - y \omega, \omega) + \\\
&\sigma_a(\mathbf{x}) L_e(\mathbf{x} - y \omega, \omega)). \end{eqnarray}$$

Writing the transmittance term in every equation will be highly inconvenient, so we can substitute it with the following expression,

$$ T(\mathbf{x}, y) = e^{-\int_0^y \mathrm{d}s \, \sigma_t(\mathbf{x} - s\omega)}, $$

which leads to the much shorter equation,

$$\begin{eqnarray} L(\mathbf{x}, \omega) = \int_0^\infty \mathrm{d}y \, T(\mathbf{x}, y) \, (&\sigma_s(\mathbf{x}) L_s(\mathbf{x} - y \omega, \omega) + \\\
&\sigma_a(\mathbf{x}) L_e(\mathbf{x} - y \omega, \omega)). \end{eqnarray}$$

The only thing left is what do we do when light hits a surface. That’s where the rendering equation or any other equation plugs into this equation. For convenience let’s call it \(L_r(\mathbf{x}, \omega)\) and add it to the equation above by introducing a boundary at distance \(d\),

$$\begin{eqnarray} L(\mathbf{x}, \omega) &=& T(\mathbf{x}, d) \, L_r(\mathbf{x}, \omega&) + \hspace{5.2cm} \\\
&& \int_0^d \mathrm{d}y \, T(\mathbf{x}, y) (&\sigma_s(\mathbf{x}) L_s(\mathbf{x} - y \omega, \omega) + \\\
&&&\sigma_a(\mathbf{x}) L_e(\mathbf{x} - y \omega, \omega)). \end{eqnarray}$$

That equation might be familiar to some of you, because that’s the main reason why you need to do premul alpha for particles. The transmittance up to the surface is \(T(\mathbf{x}, d)\) and the remaining terms are the Destination buffer (rendered with the rendering equation) and the Source buffer containing the computed scattering,

$$\begin{eqnarray} L(\mathbf{x}, \omega) &=& \mathrm{Alpha} \times \mathrm{Destination} + \mathrm{Source}. \end{eqnarray}$$

When you look at some derivations of BRDFs used in the rendering equation you might discover some similarities with that equation, and certain frameworks for layering materials are using various approximations for combining contributions of thin layers. Does that equation tell the entire truth? Not really, it is also a convenient approximation. You might notice that it is based on geometrical optics, so it doesn’t really care about wave optical phenomena. For example, the refraction of light through the atmosphere leads to the Sun appearing slightly higher than it should. But I will leave this for a future discussion. Another simplification is that the locations of particles inside the medium are assumed to be uncorrelated. There are more general equations that can even cover that case [5]. Finally, everything remains in relatively steady state, so there are no significant time correlations.

Solving the equation

The above equation can be relatively easily solved for media that is separable in multiple homogeneous volumes by doing basic Monte Carlo. The premise is that we are going to pick random distances and sample the integral and compensate for the rest of the integral by dividing out the probability of picking a specific distance,

$$\begin{eqnarray} L(\mathbf{x}, \omega) &=& T(\mathbf{x}, d) \, L_r(\mathbf{x}, \omega) & + \hspace{5.6cm} \\\
&& \frac1N \sum_{i=0}^N \, \frac{T(\mathbf{x}, y_i)}{p(y_i)} (&\sigma_s(\mathbf{x}) L_s(\mathbf{x} - y_i \omega, \omega) + \\\
&&&\sigma_a(\mathbf{x}) L_e(\mathbf{x} - y_i \omega, \omega)). \end{eqnarray}$$

Sometimes the integral is solved by breaking it up into multiple segments of equal size and evaluating the integral, i.e.

$$\begin{eqnarray} L(\mathbf{x}, \omega) &=& T(\mathbf{x}, d) \, L_r(\mathbf{x}, \omega) & + \hspace{5.6cm} \\\
&& \sum_{i=0}^N \, T(\mathbf{x}, y_i) \Delta y \, (&\sigma_s(\mathbf{x}) L_s(\mathbf{x} - y_i \omega, \omega) + \\\
&&&\sigma_a(\mathbf{x}) L_e(\mathbf{x} - y_i \omega, \omega)), \end{eqnarray}$$

so if we take into account that

$$ \Delta y = \frac{1}{N\frac1d} = \frac1N\frac{1}{p(y_i)}, $$

we arrive to the conclusion that the Monte Carlo integrator is equivalent to straight up discetizing the integral. The equidistant approach is sometimes called ray marching. However, what Monte Carlo integrators give us more than the regular discretization is that due to linearity, you can arbitrary reorder the integrated rectangular sections and if your distribution of samples is uniform you arrive to the same solution in the limit. Another interesting thing to note is that the probability density has units reciprocal to the differential that it replaces. For more complicated spaces you can also show that the Jacobian introduced by the substitution matches the probability density function.

Making your own sampling function

Suppose that you have a probability density function describing the probability of for free-flight (without collision) at a distance t (\(\mathrm{PDF}(t) = P(t < X < t + \mathrm{d}t) / \mathrm{d}t\)), then we first need the integral of the said PDF to derive the probability of free-flight to distance t (\(\mathrm{CDF}(t) = P(X \leq t)\)). From that function we can remap the uniform distribution to distances because given the random value

$$ \xi = \mathrm{CDF}(t). $$

We are actually trying to compute the distance

$$ t = \mathrm{CDF}^{-1}(\xi), $$

so in practice we are turning the diagram of our function sideways to remap from uniform values in range \([0, 1]\) to distances in range \([0, d]\).

Single scattering of collimated light in various media with analytically defined densities

Today, I am going to discuss just analytically defined media or in other words I am going to discuss various solutions of that part of the equation,

$$ L^*(\mathbf{x}, \omega) = \sum_{i=0}^N \, \frac{T(\mathbf{x}, y_i)}{p(y_i)} \sigma_s(\mathbf{x}) L_s(\mathbf{x} - y_i \omega, \omega). $$

Note that I ignore emission inside the medium because I am going to limit discussion to smoke, dust and fog. Fires and lightning are going to be ignored for simplicity.

Uniform density medium

Uniform density or constant density media are the simplest case that we can discuss. The uniform density leads to the fact that sampling either according to transmittance or the product of density and transmittance will yield the same results.

Probability density function (transmittance or product of transmittance and density)

So let’s show that what I just said is the case. We can make a PDF out of any function as long as we normalize it over the entire domain \(\mathbb{D}\) that is covered by the function, so having our completely arbitrary function g(x)

$$ PDF(x) = \frac{g(x)}{\int_\mathbb{D} g(x) \, \mathrm{d}x}. $$

That works for most bounded functions. So if we sample according to transmittance, our domain is defined as \(\mathbb{D} \in (0, +\infty)\), and thus the probability density function,

$$ PDF(x) = \frac{e^{-\sigma \, x}}{\int_{0}^{+\infty} e^{-\sigma \, x} \, \mathrm{d}x} = \frac{e^{-\sigma \, x}}{-\frac{1}{\sigma} \left. e^{-\sigma \, x} \right|_0^{+\infty}} = \sigma e^{-\sigma \, x}, $$

so we end up with the same solution, if have started directly from \(\sigma e^{-\sigma \, x}\) because the constant terms \(\sigma\) would have canceled out. The solution of the integral is, of course, one of those trivial table integrals.

Cumulative distribution function and its inverse

Relatively straightforward,

$$ \mathrm{CDF}(x) = \int_0^x \! \mathrm{PDF}(t) \, \mathrm{d}t = \int_0^x \! \sigma e^{-\sigma \, t} \, \mathrm{d}t = \left. -e^{-\sigma \, t} \right|_0^x = 1 - e^{-\sigma \, x} $$

The inverse is relatively straightform, reorder and apply logarithm on both sides,

$$ \begin{align} \xi &= 1 - e^{-\sigma \, x} \\\
1 - \xi &= e^{-\sigma \, x} \\\
\log(1 - \xi) &= -\sigma \, x \\\
x &= -\frac{\log(1 - \xi)}{\sigma} \end{align} $$

Bounding the function

Suppose that we don’t want to sample the entire domain. For most functions derived by the inversion method it is possible to compute a new PDF, CDF and inverted them by simply diving out the CDF up to our boundary. Actually, that’s just a variant of the equation that we used to derived the PDF earlier. So a bounded PDF will be,

$$ \mathrm{PDF}_\mathrm{bounded}(x, d) = \frac{\mathrm{PDF}(x)}{\mathrm{CDF}(d)}. $$

The same goes for the CDF because it is an integral of the PDF,

$$ \mathrm{CDF}_\mathrm{bounded}(x, d) = \frac{\mathrm{CDF}(x)}{\mathrm{CDF}(d)}. $$

Therefore, sampling gets modified relatively easily by rescaling,

$$ t = \mathrm{CDF}^{-1}(\mathrm{CDF}(d) \cdot \xi) = (\mathrm{CDF}^{-1})_\mathrm{bounded}(\xi, d), $$

For convenience, I will write all of them,

$$\begin{align} PDF_\mathrm{bounded}(x, d) = \frac{\sigma e^{-\sigma \, x}}{1 - e^{-\sigma \, d}} \\\
CDF_\mathrm{bounded}(x, d) = \frac{1 - e^{-\sigma \, x}}{1 - e^{-\sigma \, d}} \\
x = -\frac{\log(1 - (1 - e^{-\sigma \, d}) \, \xi)}{\sigma} \end{align}$$

In-scattering

The equation that was expressed by \(L^*(\mathbf{x}, \omega)\) is technically the in-scattered radiance weighted by the transmittance along the current ray, so actually it turns out that we can compute this integral analytically for any linear (bounded) section,

$$\begin{align} &L = e^{-(t_S + d_S) \, \sigma} \int_0^{t_B - t_S} \! e^{-\sigma \, c \, t} \, \mathrm{d}t \nonumber\\
&\hspace{10pt}= e^{-(t_S + d_S) \, \sigma} \, (1 - e^{-\sigma\,(t_B - d_S) \, c}) = E_l \, (1 - E_u), \end{align}$$

where $d_S$ is the base of the trapezoid segment at start, $d_E$ is the base at the end, $t_S$ is the start along the view vector, while $t_B$ is the end of the segment, i.e. the end of its leg. And there is a constant slope parameter c,

$$ c = 1 + \frac{d_E - d_S}{t_E - t_S}. $$

Linear density decay medium

That computation looks like it should be simple, but it turns out to be completely cursed and having many corner cases. We will assume that it linearly decays along the Y-axis and it is basically constant along X- and Z-axis.

Optical thickness

We will start with the most simple thing – the argument in the exponential. The extent of the layer is

$$ h = h_u - h_l, $$

where $h_u$ is the upper bound and $h_l$ is the lower bound. If we have the lowest point inside the medium where the computation starts $y_S$ and the distance traveled through the medium $t$, and the slope of the view vector is $v_y$, then the optical thickness is

$$ \tau = \sigma \, \frac t h \left(h_u - y_S - \frac 1 2 \, v_y \, t\right) $$

Transmittance

Pretty much straightforward when using the result in the previous section,

$$ \mathcal{T} = e^{-\tau} $$

In-scattering

And that’s where things become pretty much cursed. What is immediately noticeable about the optical thickness is that it is quadratic, so we will have a second order in the exponential. In other words we will end up with the Error function in the final solution. It is even worse, it turns complex, so you will need the Imaginary Error Function, which is very unstable to compute using a straightforward polynomial expansion. So you might need to first learn about the Dawson function and its approximations, but, anyway, here is the ugly solution of in-scattering,

$$\begin{align} l_d &= \frac{l_E - l_S}{t} \\
l_0 &= t_S - \frac 1 2 l_y l_S \\
l_1 &= -v_y - \frac 1 2 l_y l_d \\
t_2 &= - \sigma \, \frac{\frac 1 2 v_y + l_d \, l_1}{h} \\
t_1 &= \sigma \, \frac{h_S + l_S \, l_1 + l_d \, l_0}{h} \\
t_0 &= l_0 \, l_S \frac{\sigma}{h} \\
a_0 &= \frac{t_1}{2 \sqrt{t_2}} \\
a_1 &= \frac{t_1 + 2 \, t_2 \, t}{2 \sqrt{t_2}} \\
x_0 &= \begin{cases}\mathrm{erf}(a_0) & t_2 > 0 \\ \mathrm{erfi}(a_0) & \mathrm{otherwise}\end{cases} \\
x_1 &= \begin{cases}\mathrm{erf}(a_1) & t_2 > 0 \\ \mathrm{erfi}(a_1) & \mathrm{otherwise}\end{cases} \\
A &= \sqrt{\pi} \, \left(2 \, t_2 \, \frac{h_S \sigma}{h}\right) \\
B &= A \, e^{\frac{t_1^2}{4 \, t_2 \sqrt{|t_2|}} } \\
C &= A \, e^{(t_1 + 2*t_2)} \\
D &= e^{-(t_1 + t_2 \, t) \, t} \\
E &= \frac{\sigma t_S}{h} \\
L &= (1 - D) \, E + \frac{D \, C - B}{4 \, t_2 \sqrt{|t_2|}}) e^{-t_0}. \end{align}$$

We have for inconvenience also distances to the exit along the light ray at the start of the ray $l_S$ and end of the ray $l_E$ and light slope $l_y$. As it is immediately obvious, it is quite involved, but it is even more cursed, since it has singularities when $t_2 = 0$ and $t_1 = 0$ that have to be handled in a very different manner.

Sampling primary rays

In-scattering is quite awkward to invert, but it is possible to invert it for primary rays, sort of, if we are ok with transcedental functions. Usually, they can be approximated relatively well with polynomials or rational functions, if necessary. Either, polyfit or a Pade approximation might be fine for some of those using a direct solve of the mapping of the target value and CDF, i.e. solving $y - \mathrm{CDF}(x) = 0$ to find $x = \mathrm{CDF}^{-1}(y)$. So the inverse CDF is

$$\begin{align} t_2 &= - \frac 1 2 v_y \frac{\sigma}{h} \\
t_1 &= (h_u - h_S) \frac{\sigma}{h} \\
A &= \frac{1}{2 \sqrt{|t_2|}} \\
B &= \frac{A}{\sqrt{|t_2|}} \\
a_l &= \mathrm{erf}(t_1 \, A) \\
a_u &= \mathrm{erf}((t_1 + 2 \, t_2 \, t) \, A) \\
t_\mathrm{sample} &= \begin{cases}t_S + \frac{1}{\sqrt{|t_2|}} \mathrm{erf}^{-1}(a_l + (a_u - a_l) \xi) - B & t_2 > 0 \\ t_S + \frac{1}{\sqrt{|t_2|}} \mathrm{erfi}^{-1}(a_l + (a_u - a_l) \xi) - B & \mathrm{otherwise}\end{cases}. \end{align}$$

Exponential density decay medium

That kind of medium is often used to approximate the Earth atmosphere assuming that it is relatively flat around the camera, but it sort of fails during sunset, so it might require a correction curve that refits it to the values of the spherical atmosphere. What is important is to define the base height $h_0$ of the medium where the density is at its highest - equal to $\alpha$, and the decay factor $\beta$. In other words, density is defined as

$$\begin{align} \alpha e^{-\beta (h - h_0)} \end{align}$$

Optical thickness

The integral over density in that case is

$$\begin{align} \tau(t) = \frac{\alpha}{\beta} (1 - e^{-\beta \, t}), \end{align}$$

if we assume that light travels from the base towards space. Some surveys explain this case [1], but that’s not particularly useful when you render things, so we need to define it at an arbitrary point, which can be surprisingly done with the same equations and doing just a few transformations. We need to redefine the decay factor as

$$ \beta'(h) = \beta \, e^{-\beta \, h}, $$

where $h$ it the starting height. Then you need to re-weight the maximum density with the slope of the view vector,

$$ \alpha'(v_y) = \alpha \, v_y. $$

We can also compute the optical thickness to a distant light source in space because the exponetially-decaying density asymptotically approaches a limit,

$$ \tau_\mathrm{light} = \frac{\alpha'}{l_y \beta} e^{-\beta' \, t} $$

Uniform in-scattering PDF/CDF

When the medium has a boundary that is parallel to the view vector, the secondary ray doesn’t play a role and scattering is approximated with the primary ray. Actually, a lot of samplers function like that, so I will provide the solution for convenience. The PDF is

$$ PDF(t) = \frac{1}{1 - e^{-\frac{\alpha}{\beta}}} \alpha e^{-\beta t - \tau(t)}. $$

Note, that there is a normalization factor in front. Some papers forget about it and then you can’t sample from a regular uniform distribution in range [0, 1]. The CDF is just the integral, which happens to be related to the transmittance,

$$ CDF(t) = 1 - e^{-\tau(t)}. $$

Anyway, soon we will discover that having a varying distance to a collimated source is also fine.

In-scattering in the general case

The integrand is relatively straightforward. We have the optical thickness along the camera,

$$ \tau_\mathrm{camera} = \frac{\alpha'}{\beta'} (1 - e^{-\beta' \, t}), $$

and along the light,

$$ \tau_\mathrm{light} = \frac{\alpha'}{l_y \, \beta} e^{-\beta' \, t}. $$

Now the total optical thickness is just the sum,

$$ \tau = \tau_\mathrm{camera} + \tau_\mathrm{light}. $$

Transmittance is thus just the negative exponential,

$$ \mathcal{T} = e^{-\tau}. $$

And we finally need to add the density in front to get just the integrand of in-scattered light,

$$ L' = \alpha \, e^{-\beta' \, t} \, \mathcal{T}. $$

Surprisingly, integrating this integrand is possible by adding some extra factors and reduces to

$$ L(t) = \frac{1}{1 - \frac{v_y}{l_y}} (e^{-\frac{\alpha}{l_y \, \beta}} - e^{-\tau(t)}). $$

Sampling

Now we just need to turn the integral of in-scattering into a CDF and we can sample it, which also happens to be possible. We first need to normalize it, as explained earlier using the formula from the previous section,

$$ \mathrm{CDF} = \frac{L(t)}{L(t_\mathrm{max})}. $$

We can furthermore, substitute some of the constants,

$$\begin{align} A &= \frac{1}{1 - \frac{v_y}{l_y}} \frac{1}{L(t_\mathrm{max})} \\
B &= e^{-\frac{\alpha}{l_y \, \beta}} \\
\mathrm{CDF} &= (B - e^{-\tau(t)}) \, A. \end{align}$$

Then we need to substitute back the optical thickness and substitute the constants,

$$\begin{align} \mathrm{CDF} &= (B - e^{-(\frac{\alpha'}{\beta'} (1 - e^{-\beta' \, t}) + \frac{\alpha'}{l_y \, \beta} e^{-\beta' \, t})}) \, A \\
\mathrm{CDF} &= (B - e^{-(\frac{\alpha'}{\beta'} + (\frac{\alpha'}{l_y \, \beta} - \frac{\alpha'}{\beta'}) e^{-\beta' \, t})}) \, A \\
D &= \frac{\alpha'}{\beta'} \\
E &= \frac{\alpha'}{l_y \, \beta} - \frac{\alpha'}{\beta'} \\
\mathrm{CDF} &= (B - e^{-(D + E \, e^{-\beta' \, t})}) \, A \\
\end{align}$$

which happens to be relatively straightforward to invert,

$$\begin{align} \xi &= (B - e^{-(D + E \, e^{-\beta' \, t})}) \, A \\
\frac{\xi}{A} &= B - e^{-(D + E \, e^{-\beta' \, t})} \\
e^{-(D + E \, e^{-\beta' \, t})} &= B - \frac{\xi}{A} \\
-(D + E \, e^{-\beta' \, t}) &= \log{(B - \frac{\xi}{A})} \\
e^{-\beta' \, t} &= -\frac{\log{(B - \frac{\xi}{A}) - D}}{E} \\
e^{-\beta' \, t} &= \log{\left(-\frac{\log{(B - \frac{\xi}{A}) - D}}{E}\right)} \\
t &= -\frac{1}{\beta'} \log{\left(-\frac{\log{(B - \frac{\xi}{A}) - D}}{E}\right)} \end{align}$$

Where does that lead us?

Having those estimators is useful if you want to try out to importance sample the entire medium or pre-warp the domain before applying an even more sophisticated statistical frameworks. Actually, I have a relatively recent publication that I will present at i3D that will discuss some of those frameworks and you can read my paper [2] at this website.

I hope that this information will be useful for someone and it complements some of my earlier writings.

References

[1] Jan Novák, Iliyan Georgiev, Johannes Hanika, Wojciech Jarosz. Monte Carlo Methods for Volumetric Light Transport Simulation. Computer Graphics Forum (presented at EUROGRAPHICS 2018 - State of the Art Rep.), vol. 37, no. 2 Delft, Netherlands, April 16-20, 2018 link

[2] Zdravko Velinov and Kenny Mitchell. Collimated Whole Volume Light Scattering in Homogeneous Finite Media. IEEE Transactions on Visualization and Computer Graphics. (to appear) link

[3] S. Chandrasekhar, Radiative Transfer, ser. Dover Books on Intermediate and Advanced Mathematics. Dover Publications, 1960.

[4] KAJIYA J. T.: The rendering equation. In Proc. 13th Ann. Conf. on Computer Graphics and Interactive Techniques (1986), SIGGRAPH ’86, ACM, pp. 143–150 link

[5] Eugene D’Eon A Hitchhiker’s Guide to Multiple Scattering. link