10  Subsurface and Volume Scattering

This chapter studies subsurface scattering and volume scattering. While superficially different, they involve the same forms of light-matter interaction and are modeled in the same way. We start with an example to build some useful intuitions, and then we will get into the weeds of modeling. In our modeling, we start from modeling local events (absorbing and scattering photons by particles), from which we will build a general framework, called the Radiative Transfer Equation and its variant Volume Rendering Equation (VRE), to reason about subsurface and volume scattering globally. Finally, we will connect VRE to (neural) radiance-field rendering, a modern iteration of image-based rendering (Section 8.6) that uses VRE to parameterize the image formation process.

\[ \def\oi{{\omega_i}} \def\os{{\omega_s}} \def\Oi{{\Omega_i}} \def\Os{{\Omega_s}} \def\d{{\text{d}}} \def\D{{\Delta}} \def\do{{\d\omega}} \def\Do{{\Delta\omega}} \def\doi{{\d\omega_i}} \def\dos{{\d\omega_s}} \def\Doi{{\D\omega_i}} \def\Dos{{\D\omega_s}} \def\H{{\mathbf{H}}} \newcommand{\cM}{\mathcal{M}} \newcommand{\cL}{\mathcal{L}} \]

10.1 An Informal Discussion to Build Intuition

Once inside the material (through surface refraction), a photon roams about until it meets a particle. The interactions between photons and particles are governed by the subsurface scattering (SSS) or volume scattering processes. As noted before, photon emission, absorption, and scattering all take place during the SSS/volume scattering processes, not just scattering, even though the names suggest otherwise. We will generally ignore emission in our discussion unless otherwise noted, but just note that emission does happen and is correlated with absorption, since emission is the result of absorbed photons having (e.g., chemical) reactions with the particles.

Also a reminder that SSS and volume scattering are governed by exactly the same principles, because they are exactly the same thing. In computer vision and graphics literature they might be used to refer to superficially different phenomena. Volume scattering is concerned with materials that can be modeled as a volume of particles, like fog, clouds, and smoke; they are given the name participating media in computer graphics. SSS is, instead, more commonly used to refer to solids where subsurface-scattered photons contribute to their observed colors.

Subsurface scattering is so termed to distinguish itself from surface scattering, but what is beneath the surface is nothing more than a volume of particles. In fact, what is above the surface is also a volume of particles. Looking at Figure 10.1, the air, Material₁, and Material₂ can all be thought of as participating media. We usually model the air as a vacuum so photons traverse in straight lines undisturbedly, but if we were to be exact, we would want to model the particles in the air, which becomes a participating medium. So “above-surface scattering” is as different from as surface scattering as is subsurface scattering.

Figure 10.1: At the Air-Material₁ interface, photons are either reflected directly back or penetrate into the material through refraction. The refracted photons interact with the material particles through the volume scattering processes, where some photons are absorbed and others penetrate into Material₂. For someone observing from the outside, a portion of the photons would eventually leave the material composite altogether and re-enter the air. Some of these leaving photons are called the back-scattered photons that contribute to the apparent surface reflectance; others transmit through the materials and contribute to the apparent transmittance of the material composite.

10.1.1 General Intuitions

We will use Figure 10.1 as a running example to discuss the life of photons inside the material. At the Air-Material 1 interface, photons are either reflected directly back or penetrate into the material through refraction. When a refracted photon meets a particle, the particle might absorb the photon or scatter it away. If absorbed, the photon is “dead” and can be removed from the discussion. If scattered, the photon might appear to change its direction and continue to travel on a straight line until it meets another particle, so in principle a photon can be scattered multiple times.

There are three fates a photon eventually has to accept: 1) it might be absorbed along the way, 2) it might re-emerge from Material 1 back to the air, or 3) it might emerge to the air from the bottom of Material 2. Absorption is easy to understand: a photon has a certain probability of being absorbed when it meets a particle, so the longer it travels, the more likely it will be absorbed. Let’s examine the other two cases where a photon escapes the media.

• After multiple scattering, some of the initial photons that enter Material 1 from the air will reach the Material 1-air boundary again, but this time from the material side. At that point, the photons necessarily go through another round of reflection-refraction governed by the surface scattering processes. The refracted photons will re-emerge from Material 1.

  This is called back-scattering, because these photons are scattered back to where they come. As a consequence, when we observe the material from the same side of the illumination, the lights that enter our eye come from two sources: the initial surface scattering and the back-scattering.

• Some photons might leave Material 1 from the other side and enter Material 2, in which photons go through the same volume scattering processes, where some are absorbed, some can be turned back to Material 1, and some, critically, can hit the Air-Material 2 interface. Just like what happens at the Air-Material 1 interface, some of the photons will eventually emerge from Material 2. These photons essentially survive the absorption of all the particles in the media. When you observe the material from the opposite side of the illumination, it is these transmitted photons that dictate the color of the material.

Sometimes people will also say, “sub-surface scattering is caused by photons exiting at a point different from the incident point”. It points to the fact that a photon can re-emerge anywhere from the material after SSS, whereas surface scattering is modeled to be taking place only at the incident point (although we will see later that this is just a useful macroscopic abstraction or, rather, modeling strategy).

10.1.2 Transparent vs. Opaque vs. Translucent Materials

We often hear materials being described as opaque, translucent, and transparent. We can now more scientifically approach these terms given the intuitions we have built so far.

Transparent Materials

Transparent materials either scatter light predominantly in forward directions or they scatter very little light (other than surface scattering). As a result, most photons traveling through the material are either absorbed or go through without changing much of their the directions. So if you hold a transparent material against a light source, you can clearly see through the material and see the light on the other side. This does not mean transparent materials always have the same color as the light source — absorption could be wavelength-selective. An example is aqueous/dye solutions where dye molecules are very small (\(\sim nm\) range) and, thus, scatter little light so they look transparent, but depending on the absorption spectrum (which depends on how the dye molecules interact with molecules in the solvent), most dye solutions are not colorless.

Opaque Materials

In many materials, photons arriving at the material surface are either reflected right away at the surface or, for those that do penetrate into the materials, are all absorbed by the subsurface particles. Examples include conductors like metals, whose subsurface absorption is very strong, or sufficiently thick dielectrics. These materials are opaque in two senses. First, their transmittance is practically 0. Because of strong absorption, no photon re-emerges at the other side of the material. If you hold, say, a brick (dielectric) against a light bulb, the brick would completely block the light. Second, their reflectance is independent of the substrate or the material beneath them, so they completely hide the color of the substrate. Painters know that if they want to cover a layer in their painting, they will need to apply a very thick layer of paint on top.

Translucent Materials

Translucent materials such as jade, wax, and human skin are neither opaque nor transparent. If you hold wax against a light bulb, the wax will not completely block the light, so you will see some light, but you will not be able to see clearly the other side through the wax, since photons from the light bulb are very much volume-scattered after passing through the wax. Clearly modeling SSS is critical for accurately estimating the color of translucent materials. In fact, in graphics literature we sometimes see things like “modeling translucent material must consider sub-surface scattering”. In this sense, we might be tempted to classify participating media as translucent materials, because their colors certainly very much depend on volume scattering. While it is technically correct, people rarely do that, perhaps just because of the weirdness of calling, say, smokes, a material rather than a medium?

It is not true that SSS is important only in modeling translucency. Modeling SSS can be important for opaque materials. Consider the wax case: what if we make the wax very thick? The thick wax will eventually become opaque in that it will completely hide the material behind it. But that does not mean volume scattering does not matter here; the back-scattered photons do contribute to the apparent color of a thick wax.

Oil Painting Example

To put things together, consider a painting. One way paintings are characterized is by how they were painted, and we might see things like ``oil on canvas’’. Oil means the paint is oil paint, where paint pigments are dispersed into (usually linseed) oil, which is usually called the binder or the vehicle. Canvas is the substrate, which is nothing more than another material that is right beneath the painting.

The oil itself is somewhat transparent, especially when you just apply a thin layer on the canvas. But with the paint pigments, the entire oil paint becomes a translucent material. When photons leave the oil paints, they immediately interact with the canvas. If the paint layer is thick enough, virtually no photon can ever reach the canvas. But if the paint is relatively thin, the property of the substrate will contribute to the overall color of the paint. For instance, if the canvas is white-ish, a good percentage of the photons will be reflected back. The same paint would look much darker if the canvas is black, which absorbs a lot of photons.

10.1.3 Equilibrium

We can view the light-material interaction as a dynamical system under an equilibrium. To appreciate this, consider again Figure 10.1. Some photons entering Material 1 are back-scattered and hit the Air-Material 1 interface and some of those photons will re-enter Material 1 through internal reflection. Those photons will then go through multiple scattering, and as a result some will be back-scattered again and hit the Air-Material 1 interface. The cycle goes on. The secondary back-scattering is weaker in power than the first back-scattering, and the third-order back-scattering is even weaker, and so on. So eventually you can imagine that the total number of photons back-scattered at the surface will reach a constant.

In fact, this sort of dynamics takes place everywhere inside the material along every direction. If you pick a point \(p\) in the material (or at the surface) and a direction \(\omega\) starting at the point, the radiance at (\(p, \omega\)) is a constant under equilibrium. In other words, the spatial radiance distribution (a.k.a., the light field) is not changing over time.

The equilibrium is reached almost instantaneously, since light propagates incredibly fast. So the equilibrium discussion is probably of no practical impact in modeling or actual measurement, but it is still important to keep this in mind. The (spectral) reflectance/BRDF modeling/measurement is done assuming equilibrium, and later when we model volume scattering, we will set up the differential equations under the equilibrium assumption, too.

10.2 Absorption

We will focus on modeling absorption in this chapter, and the way we build the models is fundamental to how scattering will be dealt with later.

10.2.1 A Simple Case: Collimated Illumination on Uniform Medium

Imagine that a beam of light hits a volume of particles. The light is collimated in that all photons travel along the same direction. We take a slice of the material perpendicular to the incident direction. The slice is so thin that no particles in that material cover each other from the direction of the incident light. This is shown in Figure 10.2 (a). We also, for now, assume that the medium is uniform in that the number concentration \(c\) (i.e., the number of particles per unit volume) of each slice is exactly the same.

Figure 10.2: Conceptual model to help reason about photon absorption. (a): the setting for calculating the radiance reduction over a very thin slice. (b): the radiance reduction over a finite length is calculated by accumulating the radiance reduction over infinitely many thin slices.

Say the slice has a depth of \(\D s\) and a geometrical cross-sectional area of \(E\). All the particles have the same geometrical cross-sectional area of \(\epsilon_g\). In the simplest model, a photon is absorbed whenever it hits a particle. In reality, the chance of absorption can be higher or lower. The effective area available for absorption is:

\[ \begin{align} \epsilon = \epsilon_g Q_a, \end{align} \]

where \(Q_a\) is called the absorption efficiency and is usually smaller than 1 for molecules (which have small \(\epsilon_g\)) and greater than 1 for large particles (whose \(\epsilon_g\) can be large). \(Q_a\) is wavelength dependent, so we should have written it as \(Q_a(\lambda)\), but we will omit the wavelength in our notations for simplicity’s sake. In physics, \(\epsilon\) is called the absorption cross section of the particle; it characterizes the intrinsic capability of a particle to absorb photons. Mind the subtle but important difference between the geometrical cross-sectional area and the cross section of a particle.

The question we are interested in is, if the incident radiance is \(L\), what is the radiance leaving the slice \(L+\D L\)? By convention, \(\D L\) is defined as the exitant radiance minus the incident radiance and, in this case, has to be negative. The percentage of photons that are absorbed by this slice of particles (\(-\frac{\D L}{L}\)) is equivalent to the cross-sectional area of the slice that is covered by the total cross sections of the particles:

\[ \begin{align} -\frac{\D L}{L} &= \frac{cE\D s\epsilon}{E}, \end{align} \tag{10.1}\]

where \(c\) is the particle concentration of the slice, and \(E\D s\) is the total volume of the slice. So \(cE\D s\) is the number of particles in this thin slice, and \(cE\D s\epsilon\) is the total cross section of all the particles. Given the assumption that no particles are covering each other, \(\frac{cE\D l\epsilon}{E}\) is then the percentage of the thin slice’s cross-sectional area that is available for photon absorption and, thus, the percentage of the incident photons that are absorbed. The negative sign on the left-hand side of Equation 10.1 signals the fact that \(\D L\) is negative.

We rewrite Equation 10.1 as Equation 10.2: \[ \begin{align} \frac{\D L}{\D s} &= -c\epsilon L = -\sigma_a L, \end{align} \tag{10.2}\]

which shows that the amount of photon absorption per unit length (\(\frac{\D L}{\D s}\)) is proportional to the current amount of photons up to a scaling factor \(c\epsilon\). In the computer graphics literature, \(c\epsilon\) is called the absorption coefficient, denoted \(\sigma_a\).

Bouguer-Beer-Lambert’s Law

When \(\D s\) approaches infinity, we can rewrite Equation 10.2 as a differential equation:

\[ \begin{align} \frac{\d L}{\d s} = \lim_{\D s \rightarrow 0}\frac{\D L}{\D s} = -\sigma_a L \end{align} \tag{10.3}\]

This equation is a classic case of exponential decay, and its solution is given by:

\[ \begin{align} L(s) = L_0 e^{-\sigma_a s}. \end{align} \tag{10.4}\]

where \(L_0 = L(0)\) is the initial radiance of the light before interacting with the particles, as visualized in Figure 10.2 (b), \(L(s)\) denotes the radiance at a particular length \(s\).

Equation 10.4 allows us to calculate the remaining radiance after the light travels a length \(s\). Equation 10.4 is called the Bouguer-Beer-Lambert’s law (BBL), which is a geometrical optics’ simplification of the electromagnetic theory of light-matter interaction where the matter is purely absorptive (Mayerhöfer, Pahlow, and Popp 2020).

An Alternative Derivation

An equivalent way of deriving the BBL law is the following. We divide the entire volume (with a total length of \(s\)) into \(N\) thin slices, each with a length of \(\D s\). After the first slice, the surviving portion of the initial radiance is \(L = L_0(1-\sigma_a \D s)\), so after going through all the \(N\) slices, the remaining radiance is given by:

\[ \begin{align} L_N = L_0(1-\sigma_a \D s)^N = L_0(1-\sigma_a \frac{s}{N})^N. \end{align} \tag{10.5}\]

Now when \(\D s\) becomes infinitesimally small, \(N\) approaches infinity, so the limit of the remaining radiance as a function of the total length \(s\) is given in Equation 10.6, which is the same as Equation 10.4.

\[ \begin{align} L(s) = \lim_{N \rightarrow \infty}L_0(1-\sigma_a \frac{s}{N})^N = L_0 e^{-\sigma_a s}. \end{align} \tag{10.6}\]

10.2.2 Absorption Coefficient

The absorption coefficient is an important measure of the medium’s ability to absorb photons. It has a unit of \(\text{m}^\text{-1}\), which means it is not bound by 0 and 1. One way to interpret the absorption coefficient is to observe that \(\sigma_a \d s = \d L/L\), which is the fraction of the radiance absorbed or the probability of light absorption by an infinitesimal slice. So \(\sigma_a = (\d L/L)/\d s\) can be interpreted as the probability density of photon absorption, i.e., the probability of absorption per unit length traveled:

\[ \begin{align} \sigma_a = \lim_{\D s \rightarrow 0} \frac{\D L}{L}/\D s = \frac{\d L}{L \d s}. \end{align} \]

Like any density measure, absorption coefficient is most useful when it is integrated: when we integrate \(\sigma_a\) over the length that light travels, we get the fraction/percentage of the light absorbed. One can also show that \(1/\sigma_a\) is the expected value of the distance a photon can travel before being absorbed (Bohren and Clothiaux 2006, Chpt. 5.1.3); this quantity is given the name mean free path (l). To derive \(l\), observe that the probability that a photon is absorbed after traveling a distance \(s\) is \(1-e^{-\sigma_as}\). So the probability density of absorption as a function of the distance \(s\) is:

\[ \begin{align} f(s) = \frac{\text{d}(1-e^{-\sigma_a s})}{\text{d}s} = \sigma_ae^{-\sigma_a s}. \end{align} \]

So the expected value of \(s\), which we can interpret as the distance a photon can travel on average before being absorbed, is:

\[ \begin{align} l = \int_0^\infty sf(s)\text{d}s = 1/\sigma_a. \end{align} \tag{10.7}\]

10.2.3 A Few Important Quantities

We can now define a few other commonly used quantities (omitting the wavelength dependence for simplicity). The transmittance \(T\) of a volume with a total thickness of \(s\) is defined as the percentage of the transmitted/unabsorbed photons after traveling the length of \(s\) (Equation 10.8): \[ \begin{align} T = \frac{L(s)}{L_0} = e^{-\sigma_a s}. \end{align} \tag{10.8}\]

The absorbance \(A\) is the product of \(\sigma_a s\): \[ \begin{align} A = -\ln T = \ln\frac{L(s)}{L_0} = \sigma_a s. \end{align} \tag{10.9}\]

The absorptance \(a\) of a volume is defined as the percentage of the absorbed photons by the volume, which relates to \(T\) and \(A\) by Equation 10.9: \[ \begin{align} a = 1 - T = 1 - e^{-A}. \end{align} \tag{10.10}\]

We have seen these definitions in Section 3.2.1. One very nice thing about the absorbance \(A\) is that it is approximately equivalent to absorptance \(a\) when \(A\) is small (which would be true when, e.g., the length \(s\) is very small, as is the case when discussing how a photoreceptor absorbs photons when illuminated transversely).

Another nice thing about absorbance is that absorbances add, because absorption coefficients add. Imagine you have \(n\) kinds of particles mixed up in a medium, each with a different absorption coefficient \(\sigma_a^i\). The overall absorbance of the medium is the sum of the individual absorbance \(A^i\) derived as if the medium is made up of only one kind of particles. That is: \[ \begin{align} A = \sum_i^n A^i = s\sum_i^n\sigma_a^i = s\sum_i^n c^i\epsilon^i, \end{align} \tag{10.11}\]

where \(c^i\) and \(\epsilon^i\) are the concentration and absorption cross section of the \(i^{th}\) particles. Specifically, \(c^i\) is defined as: \[ \begin{align} c^i = \frac{n_i}{V}, \end{align} \]

where \(n_i\) is the number of the \(i^{th}\) kind of particles in the material, and \(V\) is the material volume.

This is not a surprising result. As long as particles in a thin slice of this new heterogeneous medium do not cover each other, we can easily extend Equation 10.1 and the rest of the derivation to consider multiple kinds of particles; eventually Equation 10.11 would be a natural conclusion. We will omit the derivation here for simplicity sake.

Equation 10.11 is a nice conclusion to have, because usually we are dealing with hybrid media. For instance, paint is a mixture of binder particles and pigment particles, and a mist is a mixture of water droplets and air particles. If we do not want to model individual matters, we can use a single absorption coefficient to describe the aggregate behavior of the mixture. That absorption coefficient does have a physical meaning: it is the concentration-weighted sum of the individual absorption coefficients.

There are a bunch of other quantities defined in the literature. The state of the definitions is a bit of a mess, largely because different communities use different definitions.

• In visual neuroscience people sometimes use a quantity called specific absorbance (see, e.g., Bowmaker and Dartnall (1980)), which is the absorbance per unit length \(\frac{A}{s}\). Whenever you see a quantity that starts with the word “specific”, chances are that the quantity is defined per unit length. You can see that specific absorbance is actually just our absorption coefficient.

• In scientific communities, especially chemistry and spectroscopy, people define \(\epsilon\), rather than \(c\epsilon\), to be the absorption coefficient. You can see the appeal of doing that — \(\epsilon\) is a more fundamental measure of a medium’s ability to absorb photons, independent of the particle concentration \(c\) (and certainly independent of the traversal length \(s\)).

• The absorbance defined in Equation 10.9 is technically called the Naperian absorbance, because we take the natural logarithm of \(T\). Sometimes people also use the decadic absorbance, which is defined as \(-\log T\). This quantity is also called the optical density.

• Finally, the number concentration \(c\) here is defined in terms of the absolute quantity per unit volume, but sometimes people want to define \(c\) as the molar concentration, which is the number of moles per unit volume. If so, all other derived quantities are then prefixed with “molar”. Next time when you see something like the molar decadic absorption coefficient, you know what it is!

The annoying thing is that people do not always tell you which definition they use. The plea I have to you is to be specific about which definition you use in your writing and tell me when I am being vague!

10.2.4 General Case

So far we have assumed that the absorption coefficient \(\sigma_a = c\epsilon\) is a constant regardless of the position \(p\) in the medium and along any direction \(\omega\). The former property assumes that the medium is uniform, and the latter property is called isotropic¹ in that the medium’s ability to absorb photons is independent of the light direction.

Both assumptions are problematic in practice. The concentration can change spatially and should be denoted \(c(p)\), where \(p\) is an arbitrary position in space. \(\epsilon\) can also change with \(p\) and, more importantly, change with the direction of light incidence \(\omega\). For instance, the particles might not be spherical, so their geometrical cross-sectional area and, thus, the cross section \(\epsilon\) available for absorbing photons can depend on \(\omega\). As a result, the absorption coefficient should generally be denoted \(\sigma_a(p, \omega)\).

Figure 10.3: A conceptual model to help reason about photon absorption in the general case, where the absorption coefficient can vary spatially and directionally. The medium is divided into many tiny elemental volumes, each of which is so small that particles do not cover each from any direction.

Effectively, our conceptual model, shown in Figure 10.3, has to be changed to one where the entire body of particles is divided into many equally-sized volumes (with a length \(\D s\) and an area \(\D A\)), each of which is so small that particles do not cover each other from any direction. The radiance reduction per unit length in a small volume is then expressed as:

\[ \begin{align} \frac{\D L(p, \omega)}{\D s} = -\sigma_a(p, \omega)L. \end{align} \tag{10.12}\]

Given this model, we can calculate the exitant radiance after light travels a length \(s\) through the medium:

\[ \begin{align} L(p+s\omega, \omega) = L(p, \omega) e^{-\int_0^s \sigma_a(p+t\omega, \omega) \d t}, \end{align} \tag{10.13}\]

where \(\omega\) is the (unit) direction of the incident radiance, \(L(p, \omega)\) is the incident radiance, and \(L(p+s\omega, \omega)\) is the exitance radiance (radiance toward \(\omega\) leaving the entire medium after traveling \(s\)).

You would notice that for a beam with an oblique incident direction, the distance traveled, say \(\D s'\), can be different (longer or shorter than) from \(\D s\). Our model can account for this by folding the factor \(\D s'/\D s\) specific to a particular direction \(\omega'\) into the absorption coefficient \(\sigma_a(p, \omega')\). Note that the \(\D s'/\D s\) factor should be the average for all the incident photons with the same direction \(\omega'\) across the entire \(\D A\).

10.2.5 Nature and Applicability of the Model

The absorption model (the BBL law) derived before (Equation 10.4 and Equation 10.13) is a continuous one, but it is derived based on modeling discrete particles and events. It is another example of the modeling methodology discussed on Section 9.5.1.

Equation 10.13 seems to suggest that absorption coefficient \(\sigma_a(p, \omega)\) is continuously defined at any position \(p\) in the medium along any direction \(\omega\). It is not true. For starters, concentration \(c\) is not continuous. Rather, it exhibits the triphasic profile shown in Figure 9.2. As we keep shrinking the size of the volume to the molecular scale, eventually the concentration depends on whether the tiny volume contains any molecules or not, so it becomes wildly discontinuous, not to mention the headache of dealing with a partial molecule in a volume — should it be counted or not? In general, the absorption coefficient can be an arbitrary discontinuous function that is not integrable.

What about Equation 10.4 where the absorption coefficient is uniform so we do not have to take the integral? Well, that is a lie too: concentration is not continuous, so it cannot be uniform everywhere, and, by extension, the absorption coefficient cannot be a constant everywhere either. So Equation 10.3 is technically wrong when we let \(\D s \rightarrow 0\) (i.e., \(N \rightarrow \infty\)), which is necessary for us to construct the differential equation (or take the limit in Equation 10.6). For Equation 10.3 to be true, the concentration/absorption coefficient must be a constant everywhere, which can be true only if the volume is continuous.

What has to happen is that the limit of \(\D s\) cannot be literally 0 and the limit of \(N\) cannot be infinity. What we do is to keep reducing \(\D s\) to the point where the concentration (and thus absorption coefficient) is insensitive to slight perturbation of \(\D s\) (i.e., operating in the stable range in Figure 9.2), and call it the concentration/absorption coefficient of that specific \(\D s\). And we repeat this for all the \(\D s\). This certainly applies to the general-case models in Equation 10.12 and Equation 10.13, where we iterate over not the thin slices \(\D s\) but all the tiny volumes (\(\D A \times \D s\)). So all the integral symbols are secretly summing over an extremely fine-grained grid.

How big of an error are we introducing here? Technically, we should sum all \(N\) slices across the total traversal length \(s\) in Equation 10.5. If we assume \(\D s\) to be very small (even though not infinitesimal) compared to \(s\), \(N\) would be large, so taking the integration (equivalent to letting \(N \rightarrow \infty\)) would be very close to summing over \(N\). Similarly, the integral in should have been a summation of the concentration in each of the \(N\) slices. If you want to be pedantic, however, the integration there is exact: we can model \(c\) as a piece-wise function, where the value at each piece is the concentration of the corresponding volume. Integrating over a piece-wise function is the same as summing all the pieces. Only the exponential expression in Equation 10.13 is inexact.

The discontinuity of the medium is, of course, orthogonal to the discontinuity and non-uniformity in the light field itself. For instance, the fact that we use \(\frac{cE\D l\epsilon}{E}\) as the percentage of photon absorption in Equation 10.1 (and implicitly in Equation 10.4) assumes that the irradiance of the incident illumination is continuous and uniform in the small volume. This is technically not true because photons are discrete packets of energy. But in practice this is not a concern because we can assume that there is an enormous amount of photons incident on the small volume, and these photons are randomly distributed.

In essence, we are using the aggregated behavior of many photons to model the behavior of a small volume. This is similar to the microfacet models, where we use the aggregated behavior of many microfacets to statistically model the behavior of a small macro-surface.

This sort of modeling strategy is a weird case where the discrete model provides the “ground truth”, which is approximated by a continuous model. I say ground truth — to the extent that the geometrical optics can approximate the electromagnetic theory of light-matter interaction. The BBL law fails when the wave nature of photons has to be considered (Mayerhöfer, Pahlow, and Popp 2020).

10.3 Scattering

Scattering is much more difficult to reason about than absorption, primarily because a scattered photon is not “dead” and continues to participate in light-matter interaction. The way to study scattering is to first understand the behavior of a single scattering event and then consider the overall behavior of a large of collection of particles.

This section focuses a single scattering event (Section 10.3.2), and the next section discusses the general case where a large collection of particles interacts with photons. Before all these, though, it is useful to first build some intuitions as to why there is a distinction between a single scattering event and scattering by a particle collection and explicitly lay out the assumptions made for the rest of our discussions (Section 10.3.1).

10.3.1 Scattering by a Particle vs. a Collection of Particles

In geometric optics terms, scattering can be thought of as an event that takes place between a photon and a particle. In the real world, however, objects and media are usually made of a large collection of particles, which introduces two complications: multiple scattering and interference.

Multiple Scattering

First, it is possible that a scattered photon, after traveling a certain distance, meets another particle and gets scattered again. This makes it considerably more difficult to analyze the effect of scattering by a medium than does the scattering of a single particle.

Figure 10.4: Left: The atmospheric scattering of the sand coming from the Sahara during Harmattan glows in the sun and gives a hazy view of the remote mountains. Nigeria’s National Mosque is in the foreground; from . Right: illustrations of the glow of the sun and the haze. Both are due to scattering, and the difference is purely visual but not fundamental.

Look at Figure 10.4 (left) taken during Harmattan, where the atmosphere is full of sand and dust blown from the Sahara. The large collection of particles in the atmosphere scatters light, glowing the sun and giving the remote mountains a hazy view. The right panel illustrates the scattering events that give rise to the glow and the haze.

Without scattering, sunlight enters the eye directly. With scattering, some photons from the sun are first knocked out of the view and could potentially be then scattered again back to the eye. Some photons that enter the eye might even come from nearby objects other than the sun. The scattering creates a glow around the sun, and, for the observer, the sun appears larger than it actually is. The hazy view of the mountains is created by exactly the same scattering processes. The photons that enter the eyes are mixed up from different parts of the mountains and from other objects. The mountains appear hazy rather than glowing as the sun does simply because the sun has a higher brightness contrast against the background than does a region on the mountain. So the distinction between “glow” and “haze” is nothing more than a visual difference at a superficial level rather than anything deeper in physics.

You can see why multiple scattering by large collections of particles poses challenges to our analysis. If a photon is scattered once in the medium, the only effect of scattering would be to knock photons out of our line of sight, and thus, remote objects would only look dimmer rather than hazy. In this case, scattering would function exactly like absorption, for modeling purposes at least. With multiple scattering, we have to track not only photons that are scattered out but also photons that are scattered into the rays that enter our eyes². This is a daunting task considering that we are usually dealing with millions of particles and billions of photons, if not more.

If you want to be absolutely pedantic, we can distinguish the following cases:

• a single scattering event, where a photon meets a particle and is scattered away;
• single scattering, where a photon is scattered once by a medium (a large collection of photons), which is under
  • a collimated illumination, so the radiance of a ray can only be weakened because photons are scattered to other directions,
  • an arbitrary illumination, so the radiance of a ray can be both weakened and augmented (by photons scattered from other directions);
• multiple scattering, where a photon is scattered multiple times by a medium, so the radiance of a ray can both be weakened and augmented.

Section 10.3.2 studies Case 1 and Case 2(a) together, because the latter is the statistical consequence of the former. Section 10.4 studies Case 2(b) and Case 3 together because they have the same observable effects and, thus, are modeled in the same way.

Interference and Coherence

Second, when there is a large collection of particles, the scattered radiation fields of individual particles can interfere with each other. The exact impact of interference can only be calculated by considering the wave nature of the light. But to the first order, the inference depends on how densely packed the particles are.

In fact, the specular surface scattering we discussed in Section 9.4 is just a macroscopic approximation of the microscopic volume scattering where particles interfere non-randomly. In a mirror or a glass of water, the particles/molecules are very densely packed to the point that the distance between two particles is smaller than the wavelength of the light. As a result, the scattering is coherent, which gives rise to the illusion of a specular surface. One can show that the Fresnel equations are the solution to the Maxwell’s equations when surface particles are densely packed.

Why would the particle density matter? If particles are very close to each other, their radiation fields are close too, so the interference is stronger and cannot be ignored. More importantly, when particles are close to each other, their spatial positions can no longer be treated as random, so the interference can become coherent. Imagine you drop particles into a vast empty space; the particle sizes are much smaller relative to the space, so their spatial distribution can be roughly described as random. But if the particles are very densely packed, where the next particle can be is very much restrained, so their positions are highly correlated, leading to coherent scattering.

We will generally assume incoherent scattering unless otherwise noted, where individual scattering events interfere each other in random ways, so we are spared of the complication of thinking of the wave nature of the photons. Under this assumption, the total power scattered by a collection of particles is the same as the sum of the power scattered by the individual particles. This happens when the particles are sufficient sufficiently distant (separated by more than multiple wavelengths) and their spatial arrangements are uncorrelated.

10.3.2 A Single Scattering Event

Scattering Efficiency and Coefficient

Intuitively, scattering has a similar effect as absorption: it weakens the radiance by taking photons away from a beam of light. The difference is that scattered photons are not dead; they are re-directed to other directions. We can define two important quantities, one to characterize a particle’s ability to scatter photons and the other to characterize a medium’s ability to scatter photons.

Similar to the situation in absorption, the intrinsic capability of a particle to scatter photons is defined by the particle’s scattering cross section \(\epsilon_s\), which itself is the product of the geometrical cross-sectional area of the particle \(\epsilon_g\) and the scattering efficiency \(Q_s\). We can then define the scattering coefficient \(\sigma_s\) of a medium (a large collection of particles), which characterizes the ability of the medium to scatter photons away from its incident radiance. \(\sigma_s\) is the product of the particle concentration of the medium \(c\) and the particle’s scattering cross section \(\epsilon_s\). Again, \(\sigma_s\) has a unit \(\text{m}^\text{-1}\) and is not bound by 0 and 1; it is best interpreted as the probability density (i.e., probability per unit length) of light being scattered away. Of course, both the scattering efficiency and scattering coefficient can vary spatially, angularly, and spectrally.

The effects of scattering and absorption add up, because they both weaken a radiance. We can extend Equation 10.13 to consider scattering (again omitting the wavelength from the equations):

\[ \begin{align} L(p+s\omega, \omega) = L(p, \omega) e^{-\int_0^s (\sigma_a(p+t\omega, \omega) + \sigma_s(p+t\omega, \omega)) \d t}, \end{align} \tag{10.14}\]

where \(\sigma_s(p, \omega)\) is the scattering coefficient at \(p\) toward the direction \(\omega\).

Rearranging the terms, we can express the transmittance between \(p\) and \(p+s\omega\) along the direction \(\omega\), denoted \(T(p \rightarrow p+s\omega)\): \[ \begin{align} T(p \rightarrow p+s\omega) = \frac{L(p+s\omega, \omega)}{L(p, \omega)} = e^{-\int_0^s (\sigma_a(p+t\omega, \omega) + \sigma_s(p+t\omega, \omega)) \d t}. \end{align} \tag{10.15}\]

Equation 10.14 can be derived using the same idea as that used for deriving the absorption equation in Section 10.2.1 by modeling a thin layer \(\D x\) — with an additional assumption that \(\D x\) is so thin that a photon is scattered at most once before leaving \(\D x\). Therefore, scattering by a single particle has the same effect as absorption: they both take the photon out of the radiance, and that is why the absorption equation (Equation 10.13) can be directly extended here.

Think of the applicability of Equation 10.14: it says that the radiance of a ray can only be weakened. If the incident light has only one direction (e.g., a collimated beam), Equation 10.14 is true when a photon is scattered at most once in the medium. This is because a scattered photon will not have a chance to get back to the ray. If the incident light is not mono-directional, e.g., diffuse illumination, Equation 10.14 in general does not apply — even if we consider only single scattering. This is because photons originally not along the direction \(\omega\) can be scattered toward it through just one single scattering event. We can see how limited Equation 10.14 is: it applies only when the illumination is collimated and we assume only single scattering. We will relax this constraint later.

Just like the absorption case (Equation 10.11), if a medium is mixed with different particles, each with a different scattering coefficient, the overall scattering coefficient is the sum of the individual scattering coefficients as if the medium is made up of a particular kind of particles.

The sum of the scattering coefficient and absorption coefficient is called the extinction coefficient or attenuation coefficient, denoted \(\sigma_t(p, \omega)\):

\[ \begin{align} \sigma_t(p, \omega) = \sigma_a(p, \omega) + \sigma_s(p, \omega). \end{align} \]

The ratio between the scattering coefficient and the attenuation coefficient is called the single-scattering albedo of the medium:

\[ \begin{align} \rho = \frac{\sigma_s(p, \omega)}{\sigma_t(p, \omega)}. \end{align} \]

This albedo can be seen as the volumetric counterpart of the surface albedo discussed in Equation 9.11. The two forms of albedo have the same physical meaning: the fraction of the incident energy that is scattered away (i.e., not absorbed). A dark medium (e.g., smoke) has a lower albedo, and a bright medium (e.g., mist) has a higher albedo.

The sum of the scattering and absorption cross sections is called the extinction cross section or attenuation cross section, denoted \(\epsilon_t = \epsilon_a + \epsilon_s\). And of course \(1/\sigma_t\) is the mean free path in a medium where both absorption and scattering take place, i.e., the mean distance a photon can travel without being absorbed or scattered away.

10.3.3 Scattering Direction Distribution: Phase Function

While the scattering efficiency (coefficient) characterizes how well a particle (medium) is able to scatter photons, it tells us nothing about the direction of scattering. The direction of a single scattering event is characterized by the phase function \(f_p(p, \os, \oi)\), which can be interpreted as the probability density function that a photon incident from a direction \(\oi\) is scattered toward a direction \(\os\). We will omit \(p\) and write the phase function as \(f_p(\os, \oi)\) when the discussion is unconcerned of \(p\).

\(f_p(\os, \oi)\) is defined as the fraction of the irradiance incident from an infinitesimal solid angle \(\doi\) that is scattered toward an infinitesimal solid angle \(\dos\) per unit solid angle:

\[ \begin{align} f_p(\os, \oi) = \lim_{\Dos \rightarrow 0} \lim_{\Doi \rightarrow 0} \frac{\D E_o(\os)}{\D E_i(\oi)} / \Dos = \frac{\d^2 E_o(\os)}{\d E_i(\oi)\dos} = \frac{\d^2 E_o(\os)}{L(\oi) \doi \dos}. \end{align} \tag{10.16}\]

\(\D E_i(\oi)\) is the incident irradiance over a small solid angle \(\Doi\) and scatters in all directions. \(\D E_o(\oi)\) is the outgoing irradiance over a small solid angle \(\Dos\), so \(\frac{\D E_o(\os)}{\D E_i(\oi)}\) is the fraction of the photons incident from \(\Doi\) that are scattered over \(\Dos\) or, alternatively, the probability that a photon incident from \(\Doi\) is scattered toward \(\Dos\); this ratio/fraction is clearly a value between 0 and 1. Dividing that fraction by \(\Dos\) gets us the probability per unit solid angle. When both the incident solid angle \(\Doi\) and the outgoing solid angle \(\Dos\) approach 0, the fraction can be interpreted as the directional-directional reflectance (Section 9.2), and the probability per solid angle within \(\Dos\) becomes the probability density toward \(\os\).

Like all density functions, the meaning of a phase function is most clear when it is integrated to compute some other quantity. Integrating Equation 10.16 over all the outgoing directions \(\os\):

\[ \begin{align} \d E_o &= \int^{\Omega = 4\pi} f_p(\os, \oi)L(\oi)\doi\dos \label{eq:energy_con_phase_func_1a} \\ &= L(\oi)\doi\int^{\Omega = 4\pi} f_p(\os, \oi)\dos \label{eq:energy_con_phase_func_1b} \\ &= \d E_i\int^{\Omega = 4\pi} f_p(\os, \oi)\dos \label{eq:energy_con_phase_func_1c}. \end{align} \tag{10.17}\]

To interpret this integration, consider a point that receives an incident radiance of \(L(\oi)\) over an infinitesimal solid angle \(\doi\). The point receives a total irradiance of \(\d E_i = L(\oi)\doi\), which is scattered in all directions. The density of the irradiance scattered toward a particular direction \(\os\) is \(f_p(\os, \oi)L(\oi)\doi\)³, which when multiplied by \(\dos\) gives us the actual irradiance scattered over a small solid angle \(\dos\) around \(\os\). Integrating all outgoing directions over the entire sphere (\(4\pi\)) we have Equation 10.17.

Now, of course, some of the photons in \(\d E_i\) might not be scattered; they could be absorbed, or they could simply not hit the cross section of any particle. So technically \(\d E_o \leq \d E_i\) in Equation 10.17, just like how energy conservation is expressed in surface scattering in Equation 9.8. By the convention in the volume scattering literature, however, the phase function is defined such that \(\d E_i\) refers to only the portion of the incident irradiance that does get scattered. Therefore, \(\d E_o = \d E_i\), so we have:

\[ \begin{align} \int^{\Omega = 4\pi} f_p(\os, \oi) \dos = \int^{\Omega = 4\pi} f_p(\os, \oi) \doi = 1. \end{align} \tag{10.18}\]

That is, the phase function integrates to 1; the second integral can be derived using the Helmholtz reciprocity (since we are still dealing with geometrical optics):

\[ \begin{align} f_p(\oi, \os) = f_p(\os, \oi). \end{align} \]

One way to interpret the fact that the phase function integrates to 1 is that the phase function is the conditional probability density function of scattering: given that a photon is scattered, what is the probability (density) of scattering to a particular direction?

Phase Function vs. BRDF

The phase function can be seen as the volumetric counterpart (in the sense that we are talking about volume scattering) of BRDF (Section 9.1) — with two differences. First, the definition of the BRDF accounts for absorption, so the BRDF integrates to at most 1, whereas the integral of the phase function is normalized to 1. This difference in definition is born purely of convention.

The second difference is more fundamental. There is no \(\cos\theta\) term when using the phase function; see, e.g., Equation 10.17, unlike how the BRDF is used to turn irradiance into radiance (e.g., Equation 9.8). In fact, from Equation 10.17 we can see that given a radiance \(L(\oi)\) and a solid angle \(\doi\), the irradiance is simply \(L(\oi)\doi\) rather than \(L(\oi)\cos\theta_i\doi\). Didn’t we say that there is a cosine fall-off between radiance and irradiance (Section 8.3.4)?

One intuition that might help is that in volume scattering we are dealing with points, which can receive flux from the entire sphere and have no definition of a normal (because points are dimensionless and shapeless) or, perhaps more conveniently, have a “flexible” normal that changes with the illumination direction and is always facing directly at the illumination. Entertain this thought experiment. We set up a small surface detector at a point and measure the power of the detector; if the incident light is parallel to the surface, the detector would receive no power, but would you say that the point does not receive any light and that the radiation field has no power? Of course not.

The fact that a parallel surface would receive no photons absolutely does not mean the illumination has no power; the radiation field is the same whether it is illuminating a surface or illuminating a point. But if we are modeling a surface, we want our model to say that the power received by the surface is 0, because it matches our phenomenological observation (that a detector arranged that way would receive no recording); when we are modeling a point in volume scattering, we want the point to receive a power as if the point has a “normal” that is directly facing the illumination because, again, this matches our phenomenological observation.

Ultimately, the difference is a conscious choice of modeling strategy even though the underlying physics is exactly the same. That is why models based on BRDF and phase function are phenomenological models. If you deal with electromagnetic theories and QED, you would not have to have this distinction between modeling surface and volume scattering.

With the understanding that there is no cosine fall-off in volume scattering, Equation 10.16 can be re-written as:

\[ \begin{align} f_p(\os, \oi) = \frac{\d^2 E_o(\os)}{\d E_i(\oi)\dos} = \frac{\d}{\d E_i(\oi)}\frac{\d E_o(\os)}{\dos} = \frac{\d L_o(\os)}{\d E_i(\oi)} = \frac{\d L_o(\os)}{L_i(\oi)\doi}, \end{align} \tag{10.19}\]

where \(\d L_o(\os)\) is the infinitesimal outgoing radiance toward \(\os\). In this sense, the phase function operates in exactly the same way as the BRDF (Equation 9.33): they both operate on irradiance and turn infinitesimal irradiance into infinitesimal radiance.

Isotropic Medium and Isotropic Scatters

Given the normalization in the phase function, the scattering efficiency should actually be parameterized as \(\bar{Q_s}(p, \os, \oi)\):

\[ \begin{align} \bar{Q_s}(p, \os, \oi) = Q_s(p, \oi) f_p(p, \os, \oi), \end{align} \tag{10.20}\]

where \(Q_s(p, \oi)\) should be be interpreted as the total scattering efficiency at \(p\) over all outgoing directions for a given incident direction \(\oi\). Similarly, the scattering coefficient would be expressed as:

\[ \begin{align} \bar{\sigma_s}(p, \os, \oi) = \sigma_s(p, \oi) f_p(p, \os, \oi), \end{align} \tag{10.21}\]

where \(\sigma_s(p, \oi)\) is interpreted as the total scattering coefficient at \(p\) over all outgoing directions for a given incident direction \(\oi\).

Figure 10.5: Visualizations of common phase functions; adapted from Novák et al. (2018), where \(\omega\) and \(\bar\omega\) are the incident direction \(\oi\) and the outgoing direction \(\os\) in our notation, respectively. For an isotropic medium, the phase function depends on only the angle \(\theta\) subtended by \(\omega\) and \(\omega'\) and is axially symmetric about \(\omega\).

Figure 10.5 visualizes a few common phase functions. While \(f_p(\cdot)\) is technically a 4D function parameterized by \(\os\) and \(\oi\), the phase function of many natural media is 1D and depends only on the angle \(\theta\) subtended by \(\os\) and \(\oi\). In Figure 10.5, the distance of a point on the contour to the center represents the magnitude of the phase function at that particular \(\theta\) Consider under what conditions this simplification can be true:

• First, it says that the phase function does not depend on the absolute incident direction \(\oi\) but the relative angle between \(\oi\) and \(\os\). To get a visual intuition, see Figure 10.6; if the phase function is invariant to the photon incident direction \(\oi\), we can, without losing any generality, assign \(\oi\) to the \(z\)-axis; the scattered direction \(\os\) is parameterized by \(\theta\) and \(\phi\).
• Second, it also says the phase function depends on only \(\theta\) but not \(\phi\). That is, the phase function is rotationally symmetric about the incident direction \(\oi\). So \(f_p(\os, \oi) = f_p(\os', \oi) \neq f_p(\os'', \oi)\).

Figure 10.6: The phase function of a spherical particle is 1) invariant to the incident direction \(\oi\), which, without losing generality, is taken to be the \(z\)-axis here, and 2) also invariant to the azimuthal angle \(\phi\) of the outgoing direction \(\os\) but depends on the polar angle \(\theta\). So \(f_p(\os, \oi) = f_p(\os', \oi) \neq f_p(\os'', \oi)\). Media consisting of such particles are called isotropic media, but it does not mean the particle itself is an isotropic scatterer, which does not exist, but if it did, its phase function would be a constant (invariant to both \(\theta\) and \(\phi\)).

Intuitively, the phase function has the following two properties:

• If you fix the incident direction, no matter how you rotate the particle, the phase function distribution is the same. Alternatively, if you change the incident direction, the phase function distribution moves along with the incident direction.
• Given an incident direction, the phase function distribution is axially symmetric about the incident direction.

The two conditions above are met only when the medium consists of randomly distributed spherically symmetric particles⁴, in which case 1) there is no reason to think that any incident direction is special, so the phase function certainly is invariant to \(\oi\), and 2) there is no reason to think \(\os\) and \(\os'\) are any different since one should not expect the scattering behavior to change if we rotate the sphere about the incident direction (\(z\)-axis).

A medium consisting of spherically symmetric particles is called a symmetric or an isotropic medium. Usually when we refer to an isotropic medium, not only is the phase function but also the total scattering coefficient \(\sigma_s(p, \oi)\) (Equation 10.21) rotationally invariant to the incident direction⁵.

As we said earlier, “isotropic” is an unbelievably overloaded term. People also call a particle an isotropic scatterer if its phase function is a constant, i.e., invariant to \(\os\); such a phase function is sometimes called an isotropic phase function. An isotropic scatterer does not exist; it is a purely theoretical construction, but if it existed, its phase function would take the value of \(\frac{1}{4\pi}\) given Equation 10.18, as shown in the first graph in Figure 10.5.

10.3.4 Common Models and General “Rules”

There are many factors that determine the exact scattering efficiency and scattering direction, which can be calculated by solving the Maxwell’s equations. We will talk about a few common models here; we focus on the intuitions while omitting the exact mathematical expressions, which can be found in standard texts. From the models, we can identify a few general “rules” or, rather, approximations under certain assumptions.

The main theory or model for a single scattering event is called the Mie scattering theory, which, strictly speaking, applies only when the particle is spherical (Sharma 2003, Chpt. 3.5.2; Bohren and Clothiaux 2006, Chpt 3.5; Melbourne 2004, Chpt. 3). Mie scattering is not somehow a different scattering process from any other scattering, and the Mie theory is nothing more than the solution to the Maxwell’s equations under certain conditions⁶.

The Mie theory predicts that the overall scattering efficiency \(Q_s\) is:

\[ \begin{align} Q_s &= \frac{8}{3}\gamma^4\big(\frac{m^2-1}{m^2+2}\big)\big[1+\frac{6}{5}\big(\frac{m^2-1}{m^2+2}\big)\gamma^2 + \cdots \big]\\ \gamma &= \frac{r}{\lambda_m}, \\ m &= \frac{n}{n_m}, \end{align} \tag{10.22}\]

where \(m\) is the the relative refractive index between the particle and the medium surrounding the particle, and \(\gamma\) is the ratio between the particle radius \(r\) and the incident light wavelength in the surrounding medium \(\lambda_m\). The notion of surrounding media might come across as a little surprising: doesn’t the material consist merely of its particles? Hardly. For instance, in paints, pigments are surrounded by binders (e.g., linseed oil in oil paints, egg yolk in tempera paints, and beeswax in encaustic paints) and usually some amount of water (except oil paints). When paint dries, some water might be evaporated, leaving pockets of air, which also contributes to the surrounding media.

We can draw a few general conclusions from the model.

Small-Particle (Rayleigh) Scattering

For small particles where \(\gamma \ll 1\) (generally when the radius is ten times smaller than the wavelength of the incident light), only the first term in Equation 10.22’s bracket matters, so the scattering efficiency is inversely proportional to \(\lambda_m^{-4}\). The inverse proportionality to \(\lambda_m^{-4}\) largely (but apparently not entirely) explains why the sky is blue and why the sun is red (Bohren and Clothiaux 2006, Chpt. 8.1). Why? First, recognize that individual molecules, such as air molecules, are usually sub-nm in size, so they scatter in this small-particle regime. Short wavelength lights from the sun are scattered by the atmospheric molecules more toward the sky and eventually enter your eyes, so if you look at the sky (against the sun) it would appear blue; when you look at the sun directly, the photons entering your eyes are mostly those unscattered ones that transmit directly through the atmosphere, and they are mostly longer-wavelength photons.

By then water molecules are also similarly small, so why would water look so different from the air? It is because water molecules are very densely packed, so their scatterings are coherent. In fact, the end result of such coherent scatterings by a collection of water molecules is that water appears specular.

The photopigments in a photoreceptor are very small in size compared to the wavelengths of visible light (each rhodopsin has a cross-section area of about (Milo and Phillips 2015, p. 144)), so they almost do not scatter lights at all, only absorption. That is why we could use microspectrophotometry (MSP) to measure a photoreceptor’s (transverse) absorption rate (Section 3.2.1): MSP measures the amount of light transmitted through a photoreceptor, and if there is little scattering, then all the photons that are not measured must be absorbed by the photoreceptor.

Scattering in the small-particle regime is also called Rayleigh scattering, which, again, is not somehow a fundamentally different scattering process, and the Rayleigh scattering theory⁷ is nothing more than a special case of the Mie scattering theory (Sharma 2003, Chpt. 3.5.1; Bohren and Clothiaux 2006. Chpt. 3.2).

The phase function in the Rayleigh regime is proportional to \(1+\cos^2\theta\), so the backward and forward scatterings are roughly equally probable. Taking the phase function into account, the scattering efficiency (in the form defined in Equation 10.20) in Rayleigh scattering is proportional to:

\[ \begin{align} Q_s \propto (\frac{r}{\lambda_m})^4 \frac{m^2-1}{m^2+2} (1 + \cos^2\theta). \end{align} \]

Impact of Particle Size

When the particle size increases, the scattering efficiency increases, initially very quickly, but eventually saturates. In fact, the Mie theory predicts that when the particle size is much larger than the wavelength (e.g., more than 100 times larger), the scattering efficiency approaches a constant 2 regardless of \(m\) and \(\lambda_m\) (Johnsen 2012, fig. 5.4). This is evident in Figure 10.7, which shows the scattering efficiency of a kind of particle as a function of particle radius (\(x\)-axis) under different \(m\) (different curves); the incident light wavelength is 500 nm.

Figure 10.7: Scattering efficiency as a function of particle radius under different relative refractive index \(m\); the incident light wavelength is 500~nm. From Johnsen (2012, fig. 5.4).

The particle size also affects the phase function. As we have discussed above, small particles in the Rayleigh regime tend to scatter photons equally in the forward and backward directions, while large particles primarily scatter photons in the forward directions. The last two graphs in Figure 10.5 show the phase functions predicted by the Mie scattering theory under different particle sizes (both are larger than that in the Rayleigh regime). The forward fraction increases as the particle size increases.

Consider the scenario in Figure 10.1, where Material₁ sits on top of Material₂, and our goal is to hide Material 2 so that the color of Material 1 is dependent only on the illumination (not the property of Material 2). There is an interesting trade-off between scattering efficiency and scattering direction here. If we want Material 1 to hide Material 2, we want the particles in Material 1 to scatter a lot of light (high scattering efficiency) backwards. If the scattering efficiency is low (so photons march on and are hindered only by absorption) or the scattering is heavy in the forward directions, photons penetrate through Material 1 and reach Material 2, which would then contribute to the overall color.

Now, to scatter a lot of light, we need the particles to be large, but then the scattering will be mostly in the forward directions. So there exists a sweet spot of the particle size that provides the highest “hiding power” for a material per unit volume. If we work out the math, we will see that the sweet spot falls roughly in the visible wavelength range. That is why most paint pigments have a diameter between 100 nm and (Bruce MacEvoy 2015). Of course, no matter how poor the hiding power is for a particular paint, if you apply enough of it, it will eventually hide whatever is behind it. Dye pigments are rather small in size (nm range), so they scatter few photons and that is why dye solutions look relatively transparent.

Impact of Refractive Index

Figure 10.7 also shows the impact of the relative refractive index \(m\) (between the particle and the surrounding media) on scattering efficiency. Generally, the scattering efficiency increases with \(m\) at all particle sizes until when the particles are so large that the scattering efficiency becomes a constant. This is supported by Equation 10.22, too (\(\frac{m^2-1}{m^2+2}\) monotonically increases and has a limit of 1).

For large particles, while \(m\) does not affect the scattering efficiency, it influences the scattering directions. When \(m\) is small, the scattering tends to be more forward, whereas when \(m\) is large, the scattering tends to be toward large angles (i.e., more photons will be back-scattered). This is why wet objects look darker (recall the unpleasant experience of accidentally spilling water on your pants). In dry paints, the medium surrounding the textile particles is air, and in wet paints it is water. \(m\) becomes smaller when the material is wet (i.e., the relative refractive-index difference becomes smaller between the textile particles and water), so most of the scattering will be forward, increasing the traversal length of photons and essentially giving photons more opportunities to be absorbed.

Aspherical Particles

What if the particle is not spherical? The Mie theory does not apply. Analytical or even numerical solutions to the Maxwell’s equations would be difficult, so perhaps a better approach is just to parameterize a model and fit it with the experimental data.

One popular one-parameter parameterization of the phase function is the Henyey–Greenstein phase function (Pharr, Jakob, and Humphreys 2018, Chpt. 11.3.1; Bohren and Clothiaux 2006, Chpt. 6.3.2), which takes the form:

\[ \begin{align} p(\theta) = \frac{1}{4\pi} \frac{1-g^2}{(1+g^2-2g\cos\theta)^{3/2}}, \end{align} \]

where \(g\) is the free parameter and is usually called the asymmetry parameter.

We hasten to emphasize that the Henyey–Greenstein function has absolutely zero physical meaning; it is designed for fitting experimental phase function data, so in the modern deep learning era, you might as well try a deep neural network. The second graph in Figure 10.5 shows one instantiation of the Henyey–Greenstein function.

10.4 Radiative Transfer Equation and Volume Rendering

So far we have assumed that a ray can only be attenuated, which can happen only when the illumination is collimated and we assume single scattering. Under this assumption, Equation 10.14 allows us to calculate any radiance in the medium (by weakening the initial radiance). General media are much more complicated: illumination can be from anywhere, and multiple scattering must be accounted for. As a result, external photons can be scattered into a ray of interest, as we have intuitively discussed in Section 10.3.1.

In the realm of geometric optics and radiometry, the general way to model lights going through a material/medium amounts to solving the so-called Radiative Transfer Equation (RTE), whose modern version was established by Chandrasekhar (1960)⁸, and Kajiya and Von Herzen (1984) was the first to introduce it to computer graphics. The RTE provides a mathematical way to express an arbitrary radiance in a medium.

10.4.1 Radiative Transfer Equation

The basic idea is to set up a differential equation to describe the (rate) of the radiance change. Given an incident radiance \(L(p, \os)\), we are interested in \(L(p+\D s \os, \os)\), the radiance after the ray has gone a small distance \(\D s\). The radiance can be:

• attenuated by the medium because of absorption;
• attenuated by the medium because photons are scattered out into other directions; this is called out-scattering in graphics;
• augmented by photons that are scattered into the ray direction from all other directions — because of multiple scattering⁹; this is called in-scattering in graphics;
• augmented because particles can emit photons.

The attenuation (reduction) of the radiance over \(\D s\) is:

\[ \begin{align} -L(p, \os) \sigma_t(p, \os) \D s. \end{align} \tag{10.23}\]

The radiance augmentation due to in-scattering is given by:

\[ \begin{align} \int^{\Omega = 4\pi} f_p(p, \os, \oi) \sigma_s(p, \os) \D s L(\oi) \doi = \sigma_s(p, \os) \D s \int^{\Omega = 4\pi} L(p, \oi) f_p(p, \os, \oi) \doi. \end{align} \tag{10.24}\]

The way to interpret Equation 10.27 is the following. \(L(p, \oi)\) is the incident radiance from a direction \(\oi\), \(L(p, \oi) \doi\) is the irradiance received from \(\doi\), of which \(\sigma_s(p, \oi)\D s L(p, \os) \doi\) is the irradiance scattered in all directions after traveling a distance \(\D s\). That portion of the scattered irradiance is multiplied by \(f_p(\os, \oi)\) to give us the radiance toward \(\os\) (see Equation 10.19). We then integrate over the entire sphere, accounting for the fact that lights can come from anywhere over the space, to obtain the total augmented radiance toward \(\os\).

If we consider emission, the total radiance augmentation is:

\[ \begin{align} \sigma_a(p, \os) \D s L_e(p, \os) + \sigma_s(p, \os) \D s \int^{\Omega = 4\pi} L(p, \oi) f_p(p, \os, \oi) \doi, \end{align} \tag{10.25}\]

where \(L_e(p, \os)\) is the emitted radiance at \(p\) toward \(\os\), so the first term represents the total emission over \(\D s\). If we let:

\[ \begin{align} L_s(p, \os) = \sigma_a(p, \os) L_e(p, \os) + \sigma_s(p, \os) \int^{\Omega = 4\pi} L(p, \oi) f_p(p, \os, \oi) \doi, \end{align} \tag{10.26}\]

the total augmentation can be simplified to:

\[ \begin{align} L_s(p, \os) \D s, \end{align} \tag{10.27}\]

where the \(L_s\) term is sometimes called the source term or source function in computer graphics, because it is the source of power at \(p\).

Combining Equation 10.23 and Equation 10.27, the net radiance change is:

\[ \begin{align} \D L(p, \os) &= L(p + \D s \os, \os) - L(p, \os) \\ &= -L(p, \os) \sigma_t(p, \os) \D s + L_s(p, \os) \D s. \end{align} \tag{10.28}\]

As \(\D s\) approaches 0, we get (assuming \(\os\) is a unit vector as in Equation 10.13 and Equation 10.14):

\[ \begin{align} \os \cdot \nabla_p L(p, \os) &= \frac{\d L(p, \os)}{\d s} = \lim_{\D s \rightarrow 0} \frac{L(p + \D s \os, \os) - L(p, \os)}{\D s} \nonumber \\ &= -\sigma_t(p, \os) L(p, \os) + L_s(p, \os), \end{align} \tag{10.29}\]

where \(\nabla_p\) denotes the gradient of \(L\) with respect to \(p\), and \(\os \cdot \nabla_p\) denotes the directional derivative, which is used because technically \(p\) and \(\os\) are both defined in a three-dimensional space, so what we are really calculating is the rate of radiance change at \(p\) along \(\os\).

Equation 10.29 is the RTE, which is an integro-differential equation, because it is a differential equation with an integral embedded. The RTE has an intuitive interpretation: if we think of radiance as the power of a ray, as a ray propagates, its power is attenuated by the medium but also augmented by “stray photons” from other rays. The latter is given by \(L_s(p, \os)\), which can be thought of as the augmentation of the radiance per unit length.

The RTE describes the rate of change of an arbitrary radiance \(L(p, \os)\). But our ultimate goal is to calculate the radiance itself? Generally the RTE has no analytical solution. There are two strategies to solve it. First, we can derive analytical solutions under certain certain assumptions and simplifications.

• For instance, the integral in Equation 10.29 can be approximated by a summation along \(N\) directions; then we can turn Equation 10.29 into a system of \(N\) differential equations to be solved. This is sometimes called the N-flux theory. We will omit a formal treatment but refer you to Bohren and Clothiaux (2006, Chpt. 6.1), Volz and Simon (2001, Chpt. 3.1.2), and Klein (2010, Chpt. 5.5) for details. You might have heard of the famous Kubelka-Munk model (Kubelka and Munk 1931b, 1931a; Kubelka 1948) widely used in modeling the color of pigment mixture; it is essentially a special case of the N-flux theory where \(N=2\), which we will discuss in Section 10.5.

• Another assumption people make is to assume that volume scattering is isotropic and can be approximated as a diffusion process. This is called the diffusion approximation (Ishimaru 1977; Ishimaru et al. 1978), which is widely used in both scientific modeling (Farrell, Patterson, and Wilson 1992; Eason et al. 1978; Schweiger et al. 1995; Boas et al. 2001) and in rendering (Stam 1995, Chpt. 7; Jensen et al. 2001; Dong et al. 2013); see Bohren and Clothiaux (2006, Chpt. 6.2) for a theoretical treatment.

The second approach deserves its own section.

10.4.2 Volume Rendering Equation

The second approach, which is particularly popular in computer graphics, is to first turn the RTE into a purely integral equation and then numerically (rather than analytically) estimate the integral using Monte Carlo integration, very similar to how the rendering equation is dealt with for surface scattering (Section 9.3).

Figure 10.8: (a): Illustration of the continuous VRE (Equation 10.30). (b): Illustration of a discrete VRE (Equation 10.31), where the integral in the continuous VRE is replaced by a summation between \(p_0\) and \(p\) at an interval of \(\D s\); \(t_i\) is the total transmittance between \(p_i\) and \(p_{i+1}\); \(L_i\) is a shorthand for \(L_s(p_i, \os)\), the source term of at \(p_i\) toward \(\os\).

The way to think of this is that in order to calculate any given radiance \(L(p, \os)\), we need to integrate all the changes along the direction \(\os\) up until \(p\). Where do we start the integration? We can start anywhere. Figure 10.8 (a) visualizes the integration process. Let’s say we want to start from a point \(p_0\), whose initial radiance toward \(\os\) is \(L_0(p_0, \os)\). Let \(p = p_0 + s\os\), where \(\os\) is a unit vector and \(s\) is the distance between \(p_0\) and \(p\). An arbitrary point \(p'\) between \(p_0\) and \(p\) would then be \(p' = p_0 + s' \os\)¹⁰.

Now we need to integrate from \(p_0\) to \(p\) by running \(s'\) from 0 to \(s\). Observe that the RTE is a form of a non-homogeneous linear differential equation, whose solution is firmly established in calculus. Without going through the derivations, its solution is:

\[ \begin{align} L(p, \os) = T(p_0 \rightarrow p)L_0(p_0, \os) + \int_{0}^{s} T(p' \rightarrow p) L_s(p', \os)\d s', \end{align} \tag{10.30}\]

where \(T(p_0 \rightarrow p)\) is the transmittance between \(p_0\) and \(p\) along \(\os\), and \(T(p' \rightarrow p)\) is the transmittance between \(p'\) and \(p\) along \(\os\). Recall the definition of transmittance in Equation 10.15: it is the remaining fraction of the radiance after attenuation by the medium after traveling the distance between two points. In our case here:

\[ \begin{align} T(p' \rightarrow p) = \frac{L(p+s\os, \os)}{L(p+s'\os, \os)} = e^{-\int_{s'}^s \sigma_t(p+t\omega, \omega) \d t}, \\ T(p_0 \rightarrow p) = \frac{L(p+s\os, \os)}{L(p, \os)} = e^{-\int_{0}^s \sigma_t(p+t\omega, \omega) \d t}, \end{align} \]

The integral equation Equation 10.30 in the graphics literature is called the volume rendering equation (VRE) or the volumetric light transport equation — the counterpart of the surface LTE (Section 9.3). Looking at the visualization in Figure 10.8 (a), the VRE has an intuitive interpretation: the radiance at \(p\) along \(\os\) is the the contribution of \(p_0\) plus and contribution of every single point between \(p_0\) and \(p\).

• The contribution of \(p_0\) is given by its initial radiance \(L_0\) weakened by the transmittance between \(p_0\) and \(p\);
• Why would a point \(p'\) between \(p_0\) and \(p\) make any contribution? It is because of the source term (Equation 10.26): \(p'\) might emit lights, and some of the in-scattered photons at \(p'\) will be scattered toward \(\os\). The contribution of \(p'\) is thus given by the source term \(L_s\) weakened by the transmittance between \(p'\) and \(p\).

The form of the VRE might appear to suggest that it is enough to accumulate along only the direct path between \(p_0\) and \(p\), which is surprising given that there are infinitely many scattering paths between \(p_0\) and \(p\) (due to multiple scattering). For instance, it appears that we consider only the outgoing radiance toward \(\os\) from \(p_0\), but \(p_0\) might have outgoing radiances over other directions, which might eventually contribute to \(L(p, \os)\) through multiple scattering. Are we ignoring them?

The answer is that the VRE implicitly accounts for all the potential paths between \(p_0\) and \(p\) because of the \(L_s\) term, which expands to . That is, every time we accumulate the contribution of a point between \(p_0\) and \(p\), we have to consider the in-scattering from all the directions at that point. Another way to interpret this is to observe that the radiance term \(L\) appears on both sides of the equation. Therefore, the VRE must be solved recursively by evaluating it everywhere in space.

Does this remind you of the rendering equation (Equation 9.18)? Indeed, the VRE can be thought of as the volumetric counterpart of the rendering equation. Similarly, we can use Monte Carlo integration to estimate it, just like how the rendering equation is dealt with — with an extra complication: the VRE has two integrals: the outer integral runs from \(p_0\) to \(p\) and, for any intermediate point \(p'\), there is an inner integral that runs from \(p'\) to \(p\) to evaluate the transmittance \(T(p' \rightarrow p)\). Therefore, we have to sample both integrands.

Similar to the situation of the rendering equation, sampling recursively would exponentially increase the number of rays to be tracked. Put it another way, since there are infinitely many paths from which a ray gains its energy due to multiple scattering, we have to integrate infinitely many paths. Again, a common solution is path tracing, for which Pharr, Jakob, and Humphreys (2023, Chpt. 14) is a great reference.

A simplification that is commonly used is to assume that there is only single scattering directly from the light source. In this way, the \(L_s\) term does not have to integrate infinitely many incident rays over the sphere but only a fixed amount of rays emitted from the light source non-recursively. This strategy is sometimes called local illumination in volume rendering, as opposed to global illumination, where one needs to consider all the possible paths of light transport. The distinction is similar to that in modeling surface scattering (Section 9.3).

10.4.3 Discrete VRE and Scientific Volume Visualization

Sometimes the VRE takes the following discrete form:

\[ \begin{align} L = \sum_{i=0}^{N-1}\big(L_i\D s\prod_{j=i+1}^{N-1}t_j\big). \end{align} \tag{10.31}\]

Equation 10.31 is the discrete version of Equation 10.30: the former turns the two integrals in the latter (both the outer integral and the inner one carried by \(T(\cdot)\)) to discrete summations using the Riemann sum over \(N\) discrete points along the ray between \(p_0\) and \(p\) at an interval of \(\D s = \frac{s}{N}\).

The notations are slightly different; Figure 10.8 (b) visualizes how this discrete VRE is expressed with the new notations.

• \(L\) is \(L(p, \os)\), the quantity to be calculated;
• \(L_i\) is a shorthand for \(L_s(p_i, \os)\), i.e., the source term (Equation 10.26) for the \(i^{th}\) point between \(p_0\) and \(p\) toward \(\os\); by definition, \(p_0\) is the \(0^{th}\) point (so \(L_0\) is the initial radiance \(L_0(p_0, \os)\) in Equation 10.30) and \(p\) is the \(N^{th}\) point;
• \(t_i\) (or more explicitly \(t(p_{i} \rightarrow p_{i+1})\)) represents the total transmittance between the \(i^{th}\) and the \((i+1)^{th}\) point and is given by \(e^{-\sigma_t(p_i, \os) \D s}\) (notice the integral in continuous transmittance Equation 10.15 is gone, because we assume the transmittance between two adjacent points is a constant in the Reimann sum);
• \(\alpha_i\) is the opacity between the \(i^{th}\) and the \((i+1)^{th}\) point, which is defined as the residual of the transmittance between the two points: \(1-t_i\).

See Max (1995, Sect. 4) or Kaufman and Mueller (2003, Sect. 6.1) for a relatively straightforward derivation of Equation 10.31, but hopefully this form of the VRE is equally intuitive to interpret from Figure 10.8 (b). It is nothing more than accumulating the contribution of each point along the ray, but now we also need to accumulate the attenuation along the way just because of how opacity is defined by convention (per step), hence the product of a sequence of the opacity residuals.

We can also re-express Equation 10.31 using opacity rather than transmittance: \[ \begin{align} L &= \sum_{i=0}^{N-1}\big(L_i\D s\prod_{j=i+1}^{N-1}(1-\alpha_j)\big) \\ &= L_{N-1} \D s + L_{N-2} \D s(1-\alpha_{N-1}) + L_{N-3} \D s(1-\alpha_{N-1})(1-\alpha_{N-2}) +~\cdots~ \\ &~~~+ L_1\D s\prod_{j=2}^{N-1}(1-\alpha_j)+ L_0\D s\prod_{j=1}^{N-1}(1-\alpha_j). \end{align} \tag{10.32}\]

The discrete VRE is usually used in the scientific visualization literature, where people are interested in visualizing data obtained from, e.g., computer tomography (CT) scans or magnetic resonance imaging (MRI). There, it is the relative color that people usually care about, not the physical quantity such as the radiance, so people sometimes lump \(L_i\D s\) together as \(C_i\) and call it the “color” of the \(i^{th}\) point. The VRE is then written as:

\[ \begin{align} C = \sum_{i=0}^{N-1}\big(C_i\prod_{j=i+1}^{N-1}(1-\alpha_j)\big). \end{align} \tag{10.33}\]

The \(C\) terms are defined in a three-dimensional RGB space, and Equation 10.33 is evaluated for the three channels separately, similar to how Equation 10.31 and Equation 10.30 are meant to be evaluated for each wavelength independently. Since color is a linear projection from the spectral radiance, the so-calculated \(C\) (all three channels) is indeed proportional to the true color, although in visualization one usually does not care about the true colors anyway (see Section 10.4.3.3).

Equation 10.33 is also called the back-to-front compositing formula in volume rendering, since it starts from \(p_0\), the farthest point on the ray to \(p\). We can easily turn the order around to start from \(p\) and end at \(p_0\) in a front-to-back fashion (\(C_{N-1}\) now corresponds to \(p_0\)):

\[ \begin{align} C = \sum_{i=0}^{N-1}\big(C_i\prod_{j=0}^{i-1}t_j\big). \end{align} \tag{10.34}\]

While theoretically equivalent, the latter is better in practice because it allows us to opportunistically terminate the integration early when, for instance, the accumulated opacity is high enough (transmittance is low enough), at which point integrating further makes little numerical contribution to the result.

Another Discrete Form of VRE

A perhaps more common way to express the discrete VRE is to approximate the transmittance \(t\) using the first two terms of its Taylor series expansion and further assume that the medium has a low albedo, i.e., \(\sigma_t \approx \sigma_a\) and \(\sigma_s \approx 0\) (that is, the medium emits and absorbs only); we have:

\[ \begin{align} & 1 - \alpha_i = t_i = t(p_{i} \rightarrow p_{i+1}) = e^{-\sigma_t(p_i, \os) \D s} = 1 - \sigma_t(p_i, \os) \D s + \frac{(\sigma_t(p_i, \os) \D s)^2}{2} - \cdots \\ \approx & 1 - \sigma_t(p_i, \os) \D s \\ \approx & 1 - \sigma_a(p_i, \os) \D s. \end{align} \]

Therefore:

\[ \begin{align} \alpha_i \approx \sigma_a(p_i, \os) \D s. \end{align} \tag{10.35}\]

Now, observe that the \(L_i\) term in Equation 10.31 is the source term in Equation 10.26, which under the low albedo assumption has only the emission term, so:

\[ \begin{align} \label{eq:vre_3} L &= \sum_{i=0}^{N-1}\big(L_i\D s\prod_{j=i+1}^{N-1}(1-\alpha_j)\big) \\ &= \sum_{i=0}^{N-1}\big(\sigma_a(p_i, \os) L_e(p_i, \os)\D s\prod_{j=i+1}^{N-1}(1-\alpha_j)\big). \end{align} \tag{10.36}\]

Now plug in Equation 10.35, we have:

\[ \begin{align} L = \sum_{i=0}^{N-1}\big(L_e(p_i, \os)\alpha_i\prod_{j=i+1}^{N-1}(1-\alpha_j)\big). \end{align} \tag{10.37}\]

If we let \(C_i = L_e(p_i, \os)\), the discrete VRE is then expressed as (Levoy 1988):

\[ \begin{align} C = \sum_{i=0}^{N-1}\big(C_i\alpha_i\prod_{j=i+1}^{N-1}(1-\alpha_j)\big). \end{align} \tag{10.38}\]

This can be interpreted as a form of alpha blending (Smith 1995), a typical trick in graphics to render transparent materials. It makes sense for our discrete VRE to reduce to alpha blending: our derivation assumes that the volume does not scatter lights, so translucent materials become transparent.

Again, Equation 10.38 is the back-to-front equation, and the front-to-back counterpart looks like:

\[ \begin{align} C = \sum_{i=0}^{N-1}\big(C_i\alpha_i\prod_{j=0}^{i-1}t_j\big). \end{align} \tag{10.39}\]

If you compare the two discrete forms in Equation 10.33 and Equation 10.38, it would appear that the two are not mutually consistent! Of course we know why: 1) Equation 10.38 applies two further approximations (low albedo and Taylor series expansion) and 2) the two \(C\) terms in the two equations refer to different physical quantities (compare Equation 10.32 with Equation 10.37).

The Second Form is More Flexible

What is the benefit of this new discrete form, comparing Equation 10.34 and Equation 10.39? Both equations can be interpreted as a form of weighted sum, where \(C_i\) is weighted by a weight \(w_i\), which is \(\prod_{j=0}^{i-1}t_j\) in the first case and \(\alpha_i\prod_{j=0}^{i-1}t_j\) in the second case. The most obvious difference is that the weights in the first case are correlated but less so in the second case. The weights are strictly decreasing as \(i\) increases in the first case, since \(t_i < 1\).

In the second case, the weights are technically independent. One way to understand this is to observe, in the second case, that \(w_0 = \alpha_0\) and \(w_{i+1} = w_i \frac{\alpha_{i+1}(1-\alpha_i)}{\alpha_i}\), so there is generally a unique assignment of the \(\alpha\) values for a given weight combination. This “flexibility” will come in handy when we can manually assign (Section 10.4.3.3) or learn the weights (Section 10.4.4). Note, however, that if we impose the constraint that \(\alpha\in[0, 1]\), we are effectively constraining the weights too.

Visualization is Not (Necessarily) Physically-Based Rendering!

These discrete VRE forms might give you the false impression that we have avoided the need to integrate infinitely many paths, because, computationally, the evaluation of the VRE comes down to a single-path summation along the ray trajectory. Not really. Calculating the \(C_i\) terms in the new formulations still requires recursion if the results are meant to be physically accurate. Of course we can sidestep this by, e.g., applying the local-illumination approximation, as mentioned before, to avoid recursion.

Scientific visualization offers another opportunity: we can simply assign values to the \(C\)s and even the \(\alpha\)s without regard to physics. The goal of visualization is to discover/highlight interesting objects and structures while de-emphasizing things that are irrelevant. So the actual colors are not as important, which gives us great latitude to determine VRE parameters.

Figure 10.9: Comparing volume rendering for visualization and rendering. (a): two examples of scientific visualization (of CT data) using volume rendering; from Sjschen (2025) and MathiasRav (2009). (b): photorealistic volume rendering (a scene from Disney’s Moana, 2016); from Gnash (2017).

Figure 10.9 compares volume-rendered images for scientific visualization (a) and for photorealistic rendering (b). In the case of visualization, the data were from CT scans. In both scans, the outer surface is not transparent but is rendered so just because we are interested in seeing the inner structures that are otherwise occluded. The user makes an executive call to assign a very low transparency to the bones in the knee model 0 but a very high transparency value to the skin and other tissues: this is not physically accurate but a good choice for this particular visualization. Photorealistic rendering, in contrast, has to be physically based and does not usually have this flexibility. See figures in Wrenninge and Bin Zafar (2011) and Fong et al. (2017) for more examples.

There is volume rendering software that would allow the users to make such an assignment depending on what the user wants to highlight and visualize. With certain constraints and heuristics, one can also procedurally assign the \(\alpha\) and \(C\) values from the raw measurement data, usually a density field (see below) acquired from whatever measurement device is used (e.g., CT scanners or MRI machines), using what is called the transfer functions in the literature¹¹. Making an assignment usually is tied to a classification problem: voxels/points of different classes should have different assignments.

Density Fields

Physically speaking, the medium in RTE/VRE is parameterized by its absorption and scattering coefficients, which are a product of cross section and concentration, which is sometimes also called the density. In physically-based volume rendering, this is indeed how the density is used from the very beginning (Kajiya and Von Herzen 1984; Blinn 1982).

In visualization where being physically accurate or photorealistic is unimportant, the notion of an attenuation coefficient loses its physical meaning; it is just a number that controls how the brightness of a point weakens. People simply call the attenuation coefficient the density (Kaufman and Mueller 2003), presumably because, intuitively, if the particle density is high the color should be dimmer. If you want to be pedantic, you might say that the attenuation coefficient depends not only on the density/concentration but also on the cross section (Equation 10.2), so how can we do that? Remember in visualization one gets to make an executive call and assign the density value (and thus control \(\alpha\)), so it does not really matter if the value itself means the physical quantity of density/concentration. This is apparent in early work that uses volume rendering for scientific visualization (Sabella 1988; Williams and Max 1992), where attenuation coefficients are nowhere to be found.

For this reason, the raw volume data obtained from raw measurement device for scientific (medical) visualization are most often called the density field, even though what is being measured is almost certainly not the density field but a field of optical properties that are related to, but certainly do not equate, density. For instance, the raw data you get from a CT scanner is actually a grid of attenuation coefficients (Bharath 2009, Chpt. 3).

10.4.4 Discrete VRE in (Neural) Radiance-Field Rendering

There is another field where the discrete VRE (especially our second form Equation 10.38) is becoming incredibly popular: (neural) radiance-field rendering. The two most representative examples are NeRF (Mildenhall et al. 2021) and 3D Gaussian Splatting (3DGS) (Kerbl et al. 2023). They are fundamentally image-based rendering or light-field rendering (Section 8.6), where they sample the light field of the scene by taking a set of images at different camera poses and learn to reconstruct the underlying light field, from which they can then re-sample a new camera pose and, thus, render the corresponding image as if it was taken at that new camera pose. This is re-branded in the modern era as “novel view synthesis”.

We will assume that you have read the two papers above so that we can focus on interpreting these radiance-field methods within the fundamental framework of physically-based rendering. Such a re-interpretation would allow us to better understand where these methods come from, how they have introduced new tweaks to physically-based rendering, and what their limitations might be.

The first thing to notice is that NeRF and 3DGS, being essentially a sampling and reconstruction method, use the discrete form of the VRE (mostly Equation 10.38) as the reconstruction function. This means they assume that the radiance is calculated by a single-path summation along the ray trajectory without the need for path tracing and solving the actual RTE/VRE. The reason they can reduce infinite paths to a single-path evaluation is very similar to that in scientific visualization, except now instead of assigning the “color” values \(C\) and opacity values \(\alpha\) (or equivalently the density field as discussed above), they train a neural network to directly learn these values.

VRE for Surface Rendering?

It is interesting to observe that both NeRF and 3DGS (and the vast majority of their later developments) use the discrete VRE form in Equation 10.38 as their forward model and can evidently do a very good job at rendering opaque surfaces. Is this surprising? Isn’t Equation 10.38 designed to render transparent materials/volumes (alpha blending)?

For opaque surfaces, the “ground truth” is the surface rendering equation (Section 9.3), which can be seen as a form of weighted sum, where the BRDF \(f_r(p, \os, \oi)\) weighs the incident light irradiance \(L(p, \oi) \cos\theta_i \text{d}\oi\). In theory, the weights are independent of each other: the value of \(f_r(p, \os, \oi)\) at different \(\oi\) can technically be completely uncorrelated. But in reality they are most likely somewhat correlated for real materials: the appearance of a material does not change dramatically when the incident light direction changes slightly. For volumes/translucent materials, the “ground truth” is the VRE, which you could also say is a weighted sum, although the weights are constrained if the \(\alpha\) values are constrained (e.g., between 0 and 1), which is the case in NeRF and 3DGS training (Section 10.4.3.2).

So effectively, when rendering opaque surfaces, we are using a form of (theoretically) constrained weighted sum to approximate another (practically) constrained weighted sum, and we hope that we can learn the approximation from a large amount of offline samples. The learned parameters (color and opacity of each point) should not be interpreted literally in the physical sense. One advantage of this parameterization is that it could be used to render volumes or translucent materials if needed, where the ground truth is the VRE, in which case the learned parameters might be more amenable to physical interpretations.

Volume Graphics vs. Point-Based Graphics

Related to volume graphics, there is also a subtly different branch of graphics called point-based graphics (PBG) (Levoy and Whitted 1985; Gross and Pfister 2011). The boundary is somewhat blurred, but given the way the two terms are usually used, we can observe a few similarities and distinctions. Both volume graphics and PBG use discrete points as the rendering primitives (as opposed to continuous surfaces such as a mesh), although the input points in volume graphics are usually placed on uniform grids (Engel et al. 2006, Chpt. 1.5.2) whereas points in PBG can be spatially arbitrary.

Traditionally, PBG is almost exclusively used for photorealistic rendering of surfaces. In fact, the points used in PBG are usually acquired as samples on continuous surfaces (Gross and Pfister 2011, Chpt. 3). PBG usually uses object-order rendering through splatting, although ray casting is used too, but RTE/VRE is not involved in the rendering process (Gross and Pfister 2011, Chpt. 6).

In contrast, the use of volume graphics is much broader. Volume rendering can be used for photorealistic rendering of participating media and translucent surfaces (by solving the RTE/VRE), or it can be used for non-photorealistic data visualization (by evaluating the single-path, discrete VRE), at which point whether the object to be rendered is called a participating medium, a translucent surface, or anything else is irrelevant, because visualization does not care much about being physically accurate.

3DGS is a somewhat interesting case. It is largely a form of PBG because the rendering primitives are unaligned surface samples, and its splatting technique (which we will discuss shortly) resembles that developed in the PBG literature (Gross and Pfister 2011, Chpt. 6.1). However, 3DGS does use the discrete VRE as the forward model. Again, as discussed just above, VRE is just a way for 3DGS to parameterize its forward mode, so the comparison with traditional volume graphics and PBG should not be taken literally.

Splatting is Signal Filtering

Splatting, initially proposed by Westover (1990) for visualizing volume data, is a common rendering technique used in PBG and 3DGS-family models. We discuss what splatting is, why it works, and how it is used in NeRF and 3DGS. We start by asking: how can we render continuous surfaces from discrete points? If we directly project the points to the sensor plane, we obviously will get holes, as shown in Figure 10.10 (a). This is, of course, not an issue if the rendering primitives are meshes (or procedurally-generated surfaces).

Figure 10.10: (a): directly projecting discrete points to the image plane would create holes in the rendered image. (b): in splatting, each point is associated with a splat or a footprint function, which can distribute the color of the point to a spatial region on the image plane. (c): splatting essentially allows signal interpolation, which amounts to first reconstructing the underlying signal from the samples (with potential anti-aliasing filtering) followed by re-sampling at new, desired positions.

The key is to realize that “meshless” does not mean surfaceless: the fact that we do not have a mesh as the rendering primitives does not mean the surface does not exist. Recall that the points used by PBG are actually samples on the surface. To render an image pixel is essentially to estimate the color of a surface point that projects to the pixel (ignoring supersampling for now). From a signal processing perspective, this is a classic problem of signal filtering: reconstruction and resampling.

That is, ideally what we need to do is to reconstruct the underlying signal, i.e., the color distribution of the continuous surface¹² (combined with an anti-aliasing pre-filter¹³) and then resample the reconstructed/filtered signal at positions corresponding to pixels in an image. This is shown in Figure 10.10 (c). %This amounts to applying a single composite filter \(F\) at new sampling positions, and The name of the game is to design proper filters. The issue of signal sampling, reconstruction, and resampling is absolutely fundamental to all forms of photorealistic rendering and not limited to PBG; Pharr, Jakob, and Humphreys (2023, Chpt. 8) and Glassner (1995, Unit II) are great references.

Another way to think of this is that the color of a surface point is very likely related to its nearby points that have been sampled as part of the rendering input, so one straightforward thing to do is to interpolate from those samples to calculate colors of new surface points. Signal interpolation is essentially signal filtering/convolution.

There is one catch. In classic signal sampling theories (think of the Nyquist-Shannon sampling theorem), samples are uniformly taken and, as a result, we can use a single reconstruction filter. But in PBG the surface samples are non-uniformly taken, so a single reconstruction filter would not work. Instead, we need a different filter for each point. The filter in the PBG parlance is called a splat, or a footprint function; it is associated with each point (surface sample) and essentially distributes the point color to a local region, enabling signal interpolation. This is shown in Figure 10.10 (b). The exact forms of the footprint functions would determine the exact forms of the signal filters. Gaussian filters are particularly common, and Gaussian splatting is a splatting method that uses Gaussian filters (Greene and Heckbert 1986; Heckbert 1989; Zwicker et al. 2001b).

From a 3D modeling perspective, instead of having a continuous mesh, the scene is now represented by a set of discrete points, each of which is represented by a 2D Gaussian distribution, which is called a surface element or a surfel (Pfister et al. 2000). Each surfel is then projected to the screen space to generate a splat as the corresponding reconstruction kernel of that surface point; the reconstruction kernel, after projection, is another Gaussian filter (which can be cascaded with an anti-aliasing pre-filter). The color of each screen space point is then calculated by summing over all the splats (each of which is, of course, scaled by the color of the sample), essentially taking a weighted sum of the colors of the neighboring surface samples (see Zwicker et al. (2001b, fig. 3) for a visualization) or, in signal processing parlance, resampling the reconstructed signal with the reconstruction kernels.

This rendering process is traditionally called surface splatting (Zwicker et al. 2001b). Yifan et al. (2019) is an early attempt to learn the surfels through a differential surface splatting process. Surface splatting is a reasonable rendering model in PBG: the “ground truth” in rendering surfaces is the rendering equation, which can also be interpreted as a weighted sum.

One can also apply the same splatting idea to volume rendering. In this case, each point in the scene is represented by a 3D Gaussian distribution, which is projected to a 2D splat in the screen space. The color of a pixel is then calculated through, critically, alpha blending the corresponding splats, not weighted sum. This is called volume splatting (Zwicker et al. 2001a). 3DGS can be seen as a differential variant of traditional volume splatting, even though it is also very effective in rendering opaque surfaces (and we have discussed the reason on Section 10.4.4.1).

10.4.5 Integrating Surface Scattering with Volume Scattering

The rendering equation governs the surface scattering or light transport in space, and the RTE/VRE governs the volume/subsurface scattering or light transport in a medium. Both processes can be involved in a real-life scene. For instance, the appearance of a translucent material like a paint or a wax is a combination of both forms of scattering/light transport (Figure 10.1). Another example would be rendering smoke against a wall.

Conceptually nothing new needs to be introduced to deal with the two forms of light transport together. Say we have an opaque surface (a wall) located with a volume (smoke) in the scene. If we want to calculate the radiance of a ray leaving a point on the wall, we would evaluate the rendering equation there, and for each incident ray, we might have to evaluate the VRE since that ray might come from the volume. In practice it amounts to extending the path tracing algorithm to account for the fact that a path might go through a volume and bounce off between surface points. See Pharr, Jakob, and Humphreys (2023, Chpt. 14.2) and Fong et al. (2017, Sect. 3) for detailed discussions.

Another approach, which is perhaps more common when dealing with translucent materials (whose appearance, of course, depends on both the surface and subsurface scattering), is through a phenomenological model based on the notion of Bidirectional Scattering Surface Reflectance Distribution Function (BSSRDF) (Nicodemus et al. 1977). The BSSRDF is parameterized as \(f_s(p_s, \os, p_i, \oi)\), describing the infinitesimal outgoing radiance at \(p_s\) toward \(\os\) given the infinitesimal power incident on \(p_i\) from the direction \(\oi\):

\[ \begin{align} f_s(p_s, \os, p_i, \oi) = \frac{\text{d}L(p_s, \os)}{\text{d}\Phi(p_i, \oi)}. \end{align} \tag{10.40}\]

BSSRDF can be seen as an extension of BRDF in that it considers the possibility that the radiance of a ray leaving \(p_s\) could be influenced by a ray incident on another point \(p_i\) due to SSS/volume scattering. Given the BSSRDF, the rendering equation can be generalized to:

\[ \begin{align} L(p_o, \os) = \int^A\int^{\Omega=2\pi} f_s(p_s, \os, p_i, \oi) L(p_i, \oi) \cos\theta_i \doi \d A, \end{align} \tag{10.41}\]

where \(L(p_o, \os)\) is the outgoing radiance at \(p_o\) toward \(\os\), \(L(p_i, \oi)\) is the incident radiance at \(p_i\) from \(\oi\), \(A\) in the outer integral is the surface area that is under illumination, and \(\Omega=2\pi\) means that each surface point receives illumination from the entire hemisphere.

We can again use path tracing and Monte Carlo integration to evaluate Equation 10.41 if we know the BSSRDF, which can, again, either be analytically derived given certain constraints and assumptions or measured (Frisvad et al. 2020). To analytically derive it, one has to consider the fact that the transfer of energy from an incident ray to an outgoing ray is the consequence of a cascade of three factors: two surface scattering (refraction) factors, one entering the material surface \(p_i\) from \(\oi\) and the other leaving the material surface at \(p_i\) toward \(p_o\), and a volume scattering factor that accounts for the subsurface scattering between the incident ray at \(p_i\) and the exiting ray at \(p_o\) (Pharr, Jakob, and Humphreys 2018, Chpt. 11.4). If all three factors have an analytical form, the final BSSRDF has an analytical form too. This is the approach that, for instance, Jensen et al. (2001) takes.

10.5 The Kubelka–Munk Model

As discussed in Section 10.4.1, the general RTE is difficult to solve, and there are two general strategies. Section 10.4.1 discusses one strategy that numerically approximates the solution using Monte Carlo methods. Another common strategy is to make some simplified assumptions and/or apply additional constraints, which would allow us to derive analytical solutions. This section discusses perhaps the most aggressive form of simplification that, nevertheless, is very widely used, especially in the printing, painting, and dye industry (and to some extent in graphics).

Kubelka and Munk, two Czechoslovakian chemists, built a phenomenological model that estimates the spectral reflectance/transmittance of a material (Kubelka and Munk 1931b, 1931a; Kubelka 1948). Their model cares only about the hemispherical-hemispherical reflectance (Section 9.2), i.e., the ratio of total flux scattered upward to the hemisphere to the total toward flux incident from the hemisphere. The K-M model also considers that the material under modeling is in immediate contact with a Lambertian substrate that reflects light uniformly over all directions. This is a common scenario in scenarios such as paints on a canvas, dyes on textiles, printer inks on paper, coatings on films, etc.

We will go through an excruciatingly long derivation, but the intuitions behind the derivation are exactly the same as that of the RTE, because the K-M model is a simplification of the RTE. The derivation allows us to make clear what simplifications we have made and, thus, when the model is and is not applicable.

10.5.1 Deriving the Model

Figure 10.11 illustrates the setup that we will use to derive the model. The first important assumption that the K-M model makes is that the every point on the material surface receives exactly the same irradiance and that the material itself is homogeneous, consisting of particles that are statistically randomly distributed and oriented. As a result, there is no difference between different positions at the same depth anywhere inside the material. Of course the radiation fields at different depths are different; at the very least, the deeper you go, the fewer photons there are due to absorption. Another way to think of this is that we are intentionally limiting our consideration to only a small area where there is no spatial difference at the same depth, so we can analyze information only at different depths rather than different positions at the same depth.

Figure 10.11: The setup of deriving the K-M model. We focus on the change of downward irradiance (\(E_{\downarrow}\)) and upward irradiance (\(E_{\uparrow}\)) as they travel a thin layer \(\D x\). We make the assumption that the material is an isotropic medium and is homogeneous spatially.

We focus on a very thin layer \(\D x\); since there are no differences between horizontal positions, we can arbitrarily pick a point at depth \(0 \leq x \leq X\) for analysis, where \(X\) is the total depth (the material surface has \(x=0\) and the material bottom has \(x=X\)). The point at \(x\) receives photons from all directions. Let’s say all the downward photons (going to the lower hemisphere) have a total irradiance of \(E_{\downarrow}(x)\) and all the upward photons (going to the upper hemisphere) have a total irradiance of \(E_{\uparrow}(x)\). The (hemispherical-hemispherical) reflectance at \(x\) is then the ratio of the two:

\[ \begin{align} R(x) = \frac{E_{\uparrow}(x)}{E_{\downarrow}(x)}, \end{align} \]

and the reflectance at the material surface (which is what the K-M model is interested in calculating) is \(R(0)\).

How do we express \(E_{\downarrow}(x)\)? As usual, we set a differential equation to describe the change of \(E_{\downarrow}(x)\). Consider the following conservation of energy when \(E_{\downarrow}(x)\) goes through the thin layer \(\D x\):

\[ \begin{align} E_{\downarrow}(x + \D x) - E_{\downarrow}(x) = & - \int^{S^2_+}\sigma_a(x, \omega) \D x L(x, \omega) \do \label{eq:km_down_1} \\ & - \int^{S^2_+}\int^{S^2_+}\sigma_s(x, \omega) f_p(x, \omega', \omega) \D x L(x, \omega) \do' \do \label{eq:km_down_2} \\ & + \int^{S^2_-}\int^{S^2_-}\sigma_s(x + \D x, \omega) f_p(x + \D x, \omega', \omega) \D x L(x + \D x, \omega) \do' \do, \label{eq:km_down_3} \end{align} \tag{10.42}\]

where \(S^2_+\) and \(S^2_-\) represent the upper and lower hemispheres, respectively; \(L(x, \omega)\) is the radiance coming from direction \(\omega\) incident on a point of depth \(x\), \(f_p(x, \omega', \omega)\) is the phase function between the incident direction \(\omega\) and outgoing direction \(\omega'\), and \(\sigma_s(x, \omega)\) is the scattering coefficient at \(x\) for an incident direction \(\omega\) (Section 10.3.2).

The interpretation of this equation is exactly the same as that of the RTE. The equation essentially says that the downward irradiance after traveling \(\D x\) is (omitting emission):

• reduced by absorption (the first negative term in Equation 10.42), which is the absorption component in the RTE;
• reduced by upward scattering (the second negative term in Equation 10.42), which is the out-scattering component in the RTE;
• increased by the downward-scattering of upward irradiance reaching the bottom of the thin layer (the positive term in Equation 10.42), which is the in-scattering component in RTE.

Now let’s take a closer look at the absorption term in Equation 10.42 and re-write it:

\[ \begin{align} \int^{S^2_+}\sigma_a(x, \omega) \D x L(x, \omega) \do = \sigma_a(x, \bar{\omega}) \D x \int^{S^2_+} L(x, \omega) \do. \end{align} \tag{10.43}\]

Equation 10.43 is derived based on the mean-value theorem¹⁴ in integral calculus, which says that there exists some value of \(\bar{\omega}\in S^2_+\) that would allow us to take the \(\sigma_a\) term out of the integral. Once we do that, the integral on the right-hand side of Equation 10.43 is simply the total downward irradiance, which we denote \(E_{\downarrow}(x)\).

If we further assume that the material absorption is isotropic, then \(\sigma_a(x, \bar{\omega})\) is independent of \(\bar{\omega}\) and is simply a function of \(x\), so we write it as \(K_{\downarrow}(x)\), which allows us to re-write Equation 10.43 as:

\[ \begin{align} \int^{S^2_+}\sigma_a(x, \omega) \D x L(x, \omega) \do = K_{\downarrow}(x) \D x E_{\downarrow}(x). \end{align} \tag{10.44}\]

\(K_{\downarrow}(x)\) is a phenomenological coefficient, but it is related to the fundamental absorption coefficient \(\sigma_a\) and, thus, carries physical meanings. This is clear from Equation 10.45:

\[ \begin{align} K_{\downarrow}(x) = \frac{\int^{S^2_+}\sigma_a(x, \omega) L(x, \omega) \do}{\D x E_{\downarrow}(x)}. \end{align} \tag{10.45}\]

\(K_{\downarrow}(x)\), by definition, represents the fraction of downward irradiance that is absorbed per unit length. Therefore, \(K_{\downarrow}(x)\) can be intuitively interpreted as the effective downward absorption coefficient at depth \(x\). Note that while the mean-value theorem tells us that some \(\bar{\omega}\in S^2_+\) exists, it does not tell us what its values is. In reality, \(K_{\downarrow}(x)\) will very much depend on \(L(x, \omega)\).

Similarly, we can rewrite the upward scattering term in Equation 10.42 by first rearranging the terms:

\[ \begin{align} &\int^{S^2_+}\int^{S^2_+}\sigma_s(x, \omega) f_p(x, \omega', \omega) \D x L(x, \omega) \do' \do \\ = &\int^{S^2_+}\Big(\int^{S^2_+}\sigma_s(x, \omega) f_p(x, \omega', \omega) \do'\Big) \D x L(x, \omega) \do. \end{align} \tag{10.46}\]

We then invoke the mean-value theorem again, which says that there exists \(\bar{\omega}\in S^2_+\) that allows us to write the upward scattering term as:

\[ \begin{align} (\int^{S^2_+}\sigma_s(x, \bar{\omega}) f_p(x, \omega', \bar{\omega}) \do' \D x) \int^{S^2_+}L(x, \omega) \do. \end{align} \tag{10.47}\]

We use \(S_{\downarrow\uparrow}(x)\) to denote the first integral in Equation 10.47: \(S_{\downarrow\uparrow}(x) = \int^{S^2_+}\sigma_s(x, \bar{\omega}) f_p(x, \omega', \bar{\omega}) \do'\), which means the upward scattering term can be simplied to:

\[ \begin{align} S_{\downarrow\uparrow}(x) \D x E_{\downarrow}(x). \end{align} \tag{10.48}\]

Combining Equation 10.46 and Equation 10.48, we get: \[ \begin{align} S_{\downarrow\uparrow}(x) = & \frac{\int^{S^2_+}\int^{S^2_+}\sigma_s(x, \omega) f_p(x, \omega', \omega) L(x, \omega) \do' \do}{\D x E_{\downarrow}(x)}. \end{align} \tag{10.49}\]

\(S_{\downarrow\uparrow}(x)\) is again a phenomenological coefficient that is related to the fundamental scattering coefficient. Its physical interpretation is clear from Equation 10.49: it is the fraction of the downward irradiance that is scattered upward per unit length. We can think of it as the effective “upward scattering coefficient of the downward irradiance”.

Finally, we can do the same thing to the downward scattering term in Equation 10.42:

\[ \begin{align} \int^{S^2_-}\int^{S^2_-}\sigma_s(x + \D x, \omega) f_p(x + \D x, \omega', \omega) \D x L(x + \D x, \omega) \do' \do \end{align} \tag{10.50}\]

by first invoking the mean-value theorem and re-express it as:

\[ \begin{align} (\int^{S^2_-}\sigma_s(x + \D x, \hat{\omega}) f_p(x + \D x, \omega', \hat{\omega}) \do' \D x) \int^{S^2_-}L(x + \D x, \omega) \do. \end{align} \tag{10.51}\]

We use \(S_{\uparrow\downarrow}(x)\) to denote the first integral in Equation 10.51. The second integral in Equation 10.51 is essentially the total upward irradiance at the depth \(x+\D x\), which we denote \(E_{\uparrow}(x + \D x)\). Therefore, the downward scattering term becomes:

\[ \begin{align} S_{\uparrow\downarrow}(x + \D x) \D x E_{\uparrow}(x + \D x). \end{align} \tag{10.52}\]

Combining Equation 10.50 and Equation 10.52, we have:

\[ \begin{align} S_{\uparrow\downarrow}(x + \D x) = & \frac{\int^{S^2_-}\int^{S^2_-}\sigma_s(x + \D x, \omega) f_p(x + \D x, \omega', \omega) \D x L(x + \D x, \omega) \do' \do}{\D x E_{\uparrow}(x + \D x)},\\ S_{\uparrow\downarrow}(x) = & \frac{\int^{S^2_-}\int^{S^2_-}\sigma_s(x, \omega) f_p(x, \omega', \omega) \D x L(x, \omega) \do' \do}{\D x E_{\uparrow}(x)}. \end{align} \tag{10.53}\]

Now plug Equation 10.44, Equation 10.48, and Equation 10.52 back into Equation 10.42, we get:

\[ \begin{align} E_{\downarrow}(x + \D x) = E_{\downarrow}(x) - K_{\downarrow}(x) \D x E_{\downarrow}(x) - S_{\downarrow\uparrow}(x) \D x E_{\downarrow}(x) + S_{\uparrow\downarrow}(x + \D x) \D x E_{\uparrow}(x + \D x). \end{align} \]

Rewrite it and take the limit as \(\D x \rightarrow 0\):

\[ \begin{align} \frac{E_{\downarrow}(x + \D x) - E_{\downarrow}(x)}{\D x} &= - K_{\downarrow}(x) E_{\downarrow}(x) - S_{\downarrow\uparrow}(x) E_{\downarrow}(x) + S_{\uparrow\downarrow}(x + \D x) E_{\uparrow}(x + \D x),\\ \frac{\d E_{\downarrow}(x)}{\d x} &= - K_{\downarrow}(x) E_{\downarrow}(x) - S_{\downarrow\uparrow}(x) E_{\downarrow}(x) + S_{\uparrow\downarrow}(x) E_{\uparrow}(x). \end{align} \tag{10.54}\]

Similarly we can express the rate of change of the upward irradiance \(E_{\uparrow}(x)\) based on the same energy conservation constraint:

\[ \begin{align} \frac{\d E_{\uparrow}(x)}{\d x} &= K_{\uparrow}(x) E_{\uparrow}(x) + S_{\uparrow\downarrow}(x) E_{\uparrow}(x) - S_{\downarrow\uparrow}(x) E_{\downarrow}(x), \end{align} \tag{10.55}\]

where the three terms on the right-hand side, again, represent the absorption, out-scattering, and in-scattering in the original RTE.

Now a few more assumptions. If we assume that the material absorption is isotropic (the total absorption per unit length is invariant to incident light direction over the entire sphere), we have \(K_{\uparrow}(x) = K_{\downarrow}(x)\) because:

\[ \begin{align} K_{\uparrow}(x) = \sigma_a(x, \hat{\omega}) = \sigma_a(x, \bar{\omega}) = K_{\downarrow}(x), \end{align} \tag{10.56}\]

for some \(\hat{\omega} \in S^2_-\) and \(\bar{\omega} \in S^2_+\).

If we further assume that 1) the material is an isotropic scattering medium, then \(\sigma_s(x, \omega)\) is invariant to \(\omega\) (the total amount of scattering per unit length does not change with the incident light direction over the entire sphere), and 2) the particles are also isotropic scatters (which theoretically do not exist), then \(f_p(x, \omega', \omega)\) is a constant (\(\frac{1}{4\pi}\)). Therefore, \(S_{\uparrow\downarrow}(x) = S_{\downarrow\uparrow}(x)\), because:

\[ \begin{align} S_{\uparrow\downarrow}(x) &= \int^{S^2_-}\sigma_s(x, \hat{\omega}) f_p(x, \omega', \hat{\omega}) \do' \\ &= \sigma_s(x, \hat{\omega}) \int^{S^2_-}f_p(x, \omega', \hat{\omega}) \do' = \sigma_s(x, \hat{\omega})/2 \\ &= \sigma_s(x, \bar{\omega}) \int^{S^2_-}f_p(x, \omega', \bar{\omega}) \do' = \sigma_s(x, \bar{\omega})/2 \\ &= S_{\downarrow\uparrow}(x), \end{align} \tag{10.57}\]

for some \(\hat{\omega} \in S^2_-\) and \(\bar{\omega} \in S^2_+\).

A few notes on the assumptions here:

• The assumption of isotropic scatters is important; assuming only an isotropic medium is not enough. This is because the integrals in integrate only the upper (or lower) hemisphere rather than the entire sphere, so if the phase function itself is not a constant, the integral result will still depend on \(\bar{\omega}\) or \(\hat{\omega}\). See the distinction between an isotropic medium and an isotropic scatter on Section 10.3.3.2.
• Alternatively, if the particles are not isotropic scatters, \(S_{\uparrow\downarrow}(x) = S_{\downarrow\uparrow}(x)\) can still hold if we assume that the upward and downward irradiance are both diffuse, i.e., the radiance \(L(x, \omega)\) is invariant to \(\omega\). The isotropic medium assumption still needs to hold. This can be proven by going back to the respective definitions of \(S_{\uparrow\downarrow}(x)\) and \(S_{\downarrow\uparrow}(x)\) in Equation 10.49 and Equation 10.53¹⁵. Interestingly, if the particles are isotropic scatters, the outgoing irradiance will necessarily be diffuse.
• Yet another way for \(S_{\uparrow\downarrow}(x) = S_{\downarrow\uparrow}(x)\) to hold is if we are considering a very idealized scenario where photons travel and are scattered only upward or downward (Bohren and Clothiaux 2006, Chpt. 5.2). If the medium is also isotropic (but does not have to consist of isotropic scatters), the fraction of the upward irradiance turned downward would be the same as the fraction of the downward irradiance turned upward, so \(S_{\uparrow\downarrow}(x) = S_{\downarrow\uparrow}(x)\).

Finally, given the assumption that the material is spatially homogeneous, both \(\sigma_a\), \(\sigma_s\), and \(f_p\) are all independent of \(x\). Therefore, we can denote \(K = K_{\uparrow}(x) = K_{\downarrow}(x)\) and \(S = S_{\uparrow\downarrow}(x) = S_{\downarrow\uparrow}(x)\). \(K\) and \(S\) are simply called the absorption and scattering coefficients, respectively, of the medium, but we now should know that they are of the phenomenological nature and, with all the simplifications above, are derived from the fundamental optical absorption/scattering properties of the medium (see Equation 10.56 and Equation 10.57).

Now we combine Equation 10.54 and Equation 10.55 and get to the famous pair of differential equations underlying the K-M model:

\[ \begin{align} \frac{\d E_{\downarrow}(x)}{\d x} &= - (K +S)E_{\downarrow}(x) + S E_{\uparrow}(x), \\ \frac{\d E_{\uparrow}(x)}{\d x} &= (K + S)E_{\uparrow}(x) - S E_{\downarrow}(x). \end{align} \tag{10.58}\]

10.5.2 The Model and Its Interpretation

Equation 10.58 gives us a pair of linear differential equations. What we are interested in solving for is \(R(0) = \frac{E_{\uparrow}(0)}{E_{\downarrow}(0)}\). The boundary condition is that \(R(X) = R_g\), which is the (assumed-to-be) known reflectance of the substrate. Solving the differential equations gives us:

\[ \begin{align} & R(0) = \frac{E_{\uparrow}(0)}{E_{\downarrow}(0)} = \frac{1-R_g[a-b\coth(bSX)]}{a-R_g+b\coth(bSX)},\\ & T(X) = \frac{E_{\downarrow}(X)}{E_{\downarrow}(0)} = \frac{b}{a\sinh(bSX)+b\cosh(bSX)},\\ & a = \frac{K+S}{S},\\ & b = \sqrt{a^2-1}, \end{align} \tag{10.59}\]

where \(R(0)\) is the hemispherical-hemispherical reflectance at the material surface and \(T(X)\) is the hemispherical-hemispherical transmittance at the bottom of the material¹⁶; \(\coth(\cdot)\) is the hyperbolic cotangent function, \(\sinh(\cdot)\) is the hyperbolic sine function, and \(\cosh(\cdot)\) is the hyperbolic cosine function. Note that we omit \(\lambda\) for simplicity but keep in mind that \(R(0), T(X), S, K, R_g, a, b\) are all spectral quantities.

Equation 10.59 is the famous K-M model. Consider the case where the material is purely absorptive and scatters little (e.g., dyes on textiles or fabrics), so \(S\) is close to 0 and the reflectance is simplified to:

\[ \begin{align} R(0) = R_g e^{-2KX} = e^{-KX} R_g e^{-KX}. \end{align} \tag{10.60}\]

We can understand \(R_X\) by decomposing it into the product of three terms, each representing a step in the overall, observed reflection. First, photons go through the material from the top down, being absorbed as they go. The percentage of photons that are still left (i.e., unabsorbed) just before they hit the substrate is \(e^{-KX}\), which is consistent with the Beer-Lambert law. Second, the substrate reflects \(R_g\) amount of light back toward the material. Finally, as the photons make their way back to the surface, they go through another round of absorption governed by the same Beer-Lambert law (Section 10.2.1).

Now consider the case where there is no substrate or when the substrate is a perfect black substrate; in both cases \(R_g = 0\). Reflectance is now:

\[ \begin{align} R_{black} = \frac{1}{a+b\coth(bSX)}. \end{align} \tag{10.61}\]

Finally, consider the case where the material is so thick that no photon reaches the substrate. In this case, the reflectance is not affected by the substrate and can be simplified to (by letting \(X\) approach infinity):

\[ \begin{align} \lim_{X \rightarrow \infty} = R_{\infty} = \frac{1}{a+b} = 1 + \frac{K}{S} - \sqrt{\Big(\frac{K}{S}\Big)^2 + 2 \frac{K}{S}}. \end{align} \tag{10.62}\]

In the painting industry, we say the paint’s “hiding is complete” when the substrate does not influence the material color. We can quantify the “hiding power” of a paint by \(H = \frac{R_{white}}{R_{black}}\), i.e., the ratio of reflectance between when the substrate is black (absorbs everything) and white (reflects everything back). If the hiding is complete, \(H\) would be 1.

You might be wondering how we know \(K\) and \(S\) of a material — we measure them. For instance, observe Equation 10.61 and Equation 10.62; we can measure the reflectance \(R_{black}\) when the substrate is nearly black and the reflectance \(R_{\infty}\) when the hiding is near complete. We can then solve the system of equations to estimate \(K\) and \(S\).

You can see that we are not actually taking a very thin layer of particles, illuminating it, and then measuring how much of the light is scattered vs. absorbed. Instead, we estimate \(K\) and \(S\) macroscopically. We have in mind a model (a set of equations or functions, if you will) that is parameterized by unknown variables. We then probe the model by giving it different inputs and measuring the outputs. From the input-output pairs, we can then estimate the unknown parameters. This is why a model so defined and derived is called a phenomenological model and the parameters (e.g., the \(K\) and \(S\) coefficients) are called the phenomenological parameters — because all we do is to observe the phenomena.

10.5.3 K-M Model for Mixture of Materials

The basic K-M model can be extended to account for materials that consist of a mixture of materials, each with different absorption/scattering behaviors. This is done by first expressing the overall absorption and scattering coefficients of the mixture and then plugging them into the K-M model.

Specifically, if we are mixing \(N\) materials, each with a phenomenological absorption and scattering coefficient \(K_i\) and \(S_i\), the overall absorption and scattering coefficient of the material is:

\[ \begin{align} K &= \sum_{i=1}^N \eta_i K_i, \label{eq:ks_mix_k} \\ S &= \sum_{i=1}^N \eta_i S_i, \label{eq:ks_mix_s} \end{align} \tag{10.63}\]

where \(\eta_i\) is the volume concentration of the \(i^{th}\) material in the overall material mixture and is defined as:

\[ \begin{align} \eta_i = \frac{V_i}{\sum_{j=1}^N V_j}, \end{align} \]

where \(V_j\) is the volume of the \(j^{th}\) material. Once we have the \(K\) and \(S\) coefficients, we simply invoke Equation 10.59 to calculate the reflectance/transmittance of the material mixture.

Why would it make sense to calculate the overall \(K\) and \(S\) using the volume-weighting method in Equation 10.63? Let’s derive it using absorption as an example. The overall absorption coefficient \(K\) of the material mixture is also the optical absorption coefficient \(\sigma_a\) of the mixture (given the set of assumptions we have made so far; see Equation 10.56), which is given by Equation 10.11 as:

\[ \begin{align} K = \sigma_a = \sum_i^N c^i\epsilon^i, \end{align} \]

where \(c^i\) and \(\epsilon^i\) are the number concentration and the absorption cross section of the \(i^{th}\) material in the material mixture. Specifically, \(c^i\) is defined as:

\[ \begin{align} c^i = \frac{n_i}{\sum_{j=1}^N V_j} = \frac{V_i}{\sum_{j=1}^N V_j} \frac{n_i}{V_i}, \end{align} \tag{10.64}\]

where \(n_i\) is the number of particles of the \(i^{th}\) material in the mixture. We assume that when we mix \(N\) materials, their volumes add. That is, the volume of the mixture is the sum of that of the individual materials. This is in general true if different materials do not chemically react and if materials do not dissolve into each other. Therefore:

\[ \begin{align} K = \sum_i^N c^i\epsilon^i \label{eq:k_mix_1} = \sum_i^N \frac{V_i}{\sum_{j=1}^N V_j} \frac{n_i}{V_i} \epsilon^i = \sum_i^N \eta_i \frac{n_i}{V_i} \epsilon^i. \end{align} \tag{10.65}\]

We then, again, use the fact (Equation 10.56) that the phenomenological absorption coefficient of the \(i^{th}\) material \(K_i\) is the same as its optical absorption coefficient \(\sigma_{a, i} = c_i \epsilon_i\), where:

• \(c_i = \frac{n_i}{V_i}\) is the number concentration of the \(i^{th}\) material on its own (note the subtle difference of \(c_i\) here and \(c^i\) in Equation 10.64, which is the number concentration of the \(i^{th}\) material in a mixture), and
• \(\epsilon_i\) is just \(\epsilon^i\): they both refer to absorption cross section of the \(i^{th}\) material, which does not change whether the material is in a mixture or not.

Therefore: \[ \begin{align} K = \sum_i^N \eta_i c_i \epsilon_i = \sum_i^N \eta_i K_i, \end{align} \tag{10.66}\]

which is exactly Equation 10.63.

10.5.4 N-Stream Model

The basic K-M model is called the two-stream or two-flux model, because it considers only the total upward irradiance and total downward irradiance. We basically have divided all the directions possible in the space into only two solid angles, the upper hemisphere and the lower hemisphere. What if we want to know the irradiance at a finer granularity (i.e., over a smaller solid angle)? We divide all the directions into more solid angles and analyze the irradiance change in them using a similar method.

For instance, if we now consider irradiance in four directions: \(E_{\nwarrow}(x)\), \(E_{\nearrow}(x)\), \(E_{\searrow}(x)\), and \(E_{\swarrow}(x)\), we can write the change of the irradiance in \(E_{\searrow}(x)\) as:

\[ \begin{align} \frac{\d E_{\searrow}(x)}{\d x} = & - K_{\searrow}(x) E_{\searrow}(x) \nonumber \\ & + S_{\nwarrow\searrow}(x) E_{\nwarrow}(x) + S_{\nearrow\searrow}(x) E_{\nearrow}(x) + S_{\swarrow\searrow}(x) E_{\swarrow}(x) \nonumber \\ & - S_{\searrow\nwarrow}(x) E_{\searrow}(x) - S_{\searrow\nearrow}(x) E_{\searrow}(x) - S_{\searrow\swarrow}(x) E_{\searrow}(x), \end{align} \tag{10.67}\]

where \(S_{\nwarrow\searrow}\) is the scattering coefficient of from the northwest irradiance to the southeast direction, and so on. We can express the changes of the other three directions similarly. In general, if we divide the space into \(N\) “channels”, each representing a set of directions (a finite solid angle), we can extend the two-flux model to a “N-flux” model.

Figure 10.12: The setup for the N-flux model. All the possible directions (consider all the arrows that can possibly go out from the origin) are divided into “channels” or “streams”, each of which represents a finite solid angle within which an irradiance travels. The N-flux model models the changes of each of these irradiances. Adapted from Li (2003, fig. 4.1).

Figure 10.12 shows the setup for deriving the N-flux model, where all the possible directions (consider all the arrows that can possibly go out from the origin) are divided into “channels” or “streams”, each of which represents a finite solid angle within which an irradiance travels. The figure visualizes eight such channels. The change of each channel is modeled by taking away from each channel photons that are absorbed and scattered to all other channels and by adding photons scattered into the channel from all other channels. In the end, we get \(N\) linear differential equations. Equation 10.67 is one such equation when \(N=4\). The channels are usually rotationally symmetric about the \(z\)-axis, assuming that the material is rotationally symmetric about the \(z\)-axis.

Comparing against the RTE in Equation 10.29, which has an integral term \(L_s\) given in Equation 10.26. What the N-flux model does is to essentially approximate that integral with a finite sum. In general, the larger the \(N\) the better the approximation but also the more computationally intensive to solve. We will omit a formal treatment but refer you to Bohren and Clothiaux (2006, Chpt. 6.1), Volz and Simon (2001, Chpt. 3.1.2), and Klein (2010, Chpt. 5.5) for details.

10.5.5 Correction for Surface Reflection

One thing that is ignored in the K-M model is the surface reflection/refraction yet. Part of the photons will be reflected away at the air/material interface before they enter the material. Similarly, when we consider the transmittance, we have ignored the reflection/refraction at the other side of the material. So the reflectance and transmittance calculated by the K-M model are defined at the point when the photons are just about to leave the material.

To account for the surface phenomena, we can apply what is called the Saunderson correction, derived by Saunderson (1942). See Sharma (2003, Chpt. 3.6.3) for a derivation, but briefly, if the illumination is diffuse, the corrected surface reflectance is:

\[ \begin{align} R = r_s + \frac{(1-r_s)(1-r_i)~R(0)}{1-r_i~R(0)}, \end{align} \]

where \(R(0)\) is the reflectance given by the K-M model, \(r_s\) is the fraction of incident irradiance scattered by the air-material surface, and \(r_i\) is the fraction of the internal irradiance approaching the air-material interface that is scattered back by the interface. Assuming a smooth surface, both \(r_s\) and \(r_i\) can be calculated by the Fresnel equations given the refractive index of the material (Section 9.4).

Bharath, Anil. 2009. Introductory Medical Imaging. Morgan & Claypool.

Blinn, James F. 1982. “Light Reflection Functions for Simulation of Clouds and Dusty Surfaces.” Acm SIGGRAPH Computer Graphics 16 (3): 21–29.

Boas, David A, Dana H Brooks, Eric L Miller, Charles A DiMarzio, Misha Kilmer, Richard J Gaudette, and Quan Zhang. 2001. “Imaging the Body with Diffuse Optical Tomography.” IEEE Signal Processing Magazine 18 (6): 57–75.

Bohren, Craig F, and Eugene E Clothiaux. 2006. Fundamentals of Atmospheric Radiation: An Introduction with 400 Problems. John Wiley & Sons.

Bowmaker, James K, and HJk Dartnall. 1980. “Visual Pigments of Rods and Cones in a Human Retina.” The Journal of Physiology 298 (1): 501–11.

Bruce MacEvoy. 2015. “The material attributes of paints.” https://www.handprint.com/HP/WCL/pigmt3.html#particlesize.

Chandrasekhar, Subrahmanyan. 1960. Radiative Transfer. Courier Corporation.

Dong, Yue, Stephen Lin, Baining Guo, et al. 2013. Material Appearance Modeling: A Data-Coherent Approach. Springer.

Eason, G, AR Veitch, RM Nisbet, and FW Turnbull. 1978. “The Theory of the Back-Scattering of Light by Blood.” Journal of Physics D: Applied Physics 11 (10): 1463.

Engel, Klaus, Markus Hadwiger, Joe M Kniss, Christof Rezk-Salama, and Daniel Weiskopf. 2006. Real-Time Volume Graphics. A K Peters, Ltd.

Farrell, Thomas J, Michael S Patterson, and Brian Wilson. 1992. “A Diffusion Theory Model of Spatially Resolved, Steady-State Diffuse Reflectance for the Noninvasive Determination of Tissue Optical Properties in Vivo.” Medical Physics 19 (4): 879–88.

Fong, Julian, Magnus Wrenninge, Christopher Kulla, and Ralf Habel. 2017. “Production Volume Rendering: Siggraph 2017 Course.” In ACM SIGGRAPH 2017 Courses, 1–97.

Frisvad, Jeppe Revall, Soeren A Jensen, Jonas Skovlund Madsen, António Correia, Li Yang, Søren Kimmer Schou Gregersen, Youri Meuret, and P-E Hansen. 2020. “Survey of Models for Acquiring the Optical Properties of Translucent Materials.” In Computer Graphics Forum, 39:729–55. 2. Wiley Online Library.

Glassner, Andrew S. 1995. Principles of Digital Image Synthesis. Elsevier.

Gnash. 2017. “Rendering of the Extremely Large Telescope from 2009; CC BY-SA 4.0 license.” https://commons.wikimedia.org/wiki/File:Latest_Rendering_of_the_E-ELT.jpg.

Greene, Ned, and Paul S Heckbert. 1986. “Creating Raster Omnimax Images from Multiple Perspective Views Using the Elliptical Weighted Average Filter.” IEEE Computer Graphics and Applications 6 (6): 21–27.

Gross, Markus, and Hanspeter Pfister. 2011. Point-Based Graphics. Elsevier.

Heckbert, Paul S. 1989. “Fundamentals of Texture Mapping and Image Warping. Master’s Thesis.” University of California, Berkeley.

Ishimaru, Akira. 1977. “Theory and Application of Wave Propagation and Scattering in Random Media.” Proceedings of the IEEE 65 (7): 1030–61.

Ishimaru, Akira et al. 1978. Wave Propagation and Scattering in Random Media. Vol. 2. Academic press New York.

Jensen, Henrik Wann, Stephen R Marschner, Marc Levoy, and Pat Hanrahan. 2001. “A Practical Model for Subsurface Light Transport.” ACM SIGGRAPH Computer Graphics, 511–18.

Johnsen, Sönke. 2012. The Optics of Life: A Biologist’s Guide to Light in Nature. Princeton University Press.

Kajiya, James T, and Brian P Von Herzen. 1984. “Ray Tracing Volume Densities.” ACM SIGGRAPH Computer Graphics 18 (3): 165–74.

Kaufman, Arie, and Klaus Mueller. 2003. “Volume Visualization and Volume Graphics.” Technical Report; Stony Brook University.

Kerbl, Bernhard, Georgios Kopanas, Thomas Leimkühler, and George Drettakis. 2023. “3d Gaussian Splatting for Real-Time Radiance Field Rendering.” ACM Trans. Graph. 42 (4): 139–31.

Klein, GA. 2010. “Industrial Color Physics.” Springer Science+ Business Media.

Kubelka, Paul. 1948. “New Contributions to the Optics of Intensely Light-Scattering Materials. Part i.” Josa 38 (5): 448–57.

Kubelka, Paul, and Franz Munk. 1931a. “An Article on Optics of Paint Layers (Translated by Stephen h. Westin).” Z. Tech. Phys 12 (593-601): 259–74.

———. 1931b. “Ein Beitrag Zur Optik Der Farbanstriche.” Z. Tech. Phys 12:593–601.

Levoy, Marc. 1988. “Display of Surfaces from Volume Data.” IEEE Computer Graphics and Applications 8 (3): 29–37.

Levoy, Marc, and Turner Whitted. 1985. “The Use of Points as a Display Primitive.” Technical Report; University of North Carolina at Chapel Hill.

Li, Yang. 2003. “Ink-Paper Interaction-a Study in Ink-Jet Color Reproduction.” Institute of Technology-Linköpings University. Norrköping, Sweden: UniTryck.

MathiasRav. 2009. “Volume rendering of a mouse skull (CT) using shear warp algorithm; CC BY-SA 3.0 license.” https://en.wikipedia.org/wiki/File:VolRenderShearWarp.gif.

Max, Nelson. 1995. “Optical Models for Direct Volume Rendering.” IEEE Transactions on Visualization and Computer Graphics 1 (2): 99–108.

Mayerhöfer, Thomas G, Susanne Pahlow, and Jürgen Popp. 2020. “The Bouguer-Beer-Lambert Law: Shining Light on the Obscure.” ChemPhysChem 21 (18): 2029–46.

Melbourne, William G. 2004. Radio Occultations Using Earth Satellites: A Wave Theory Treatment. 1st ed. John Wiley & Sons.

Mildenhall, Ben, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. 2021. “Nerf: Representing Scenes as Neural Radiance Fields for View Synthesis.” Communications of the ACM 65 (1): 99–106.

Milo, Ron, and Rob Phillips. 2015. Cell Biology by the Numbers. Garland Science.

Nicodemus, FE, JC Richmond, JJ Hsia, IW Ginsberg, and T Limperis. 1977. Geometrical Considerations and Nomenclature for Reflectance. Vol. 160. US Department of Commerce, National Bureau of Standards Washington, DC, USA.

Novák, Jan, Iliyan Georgiev, Johannes Hanika, Jaroslav Křivánek, and Wojciech Jarosz. 2018. “Monte Carlo Methods for Physically Based Volume Rendering.” In ACM SIGGRAPH 2018 Courses.

Pfister, Hanspeter, Matthias Zwicker, Jeroen Van Baar, and Markus Gross. 2000. “Surfels: Surface Elements as Rendering Primitives.” In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, 335–42.

Pharr, Matt, Wenzel Jakob, and Greg Humphreys. 2018. Physically Based Rendering: From Theory to Implementation. 3rd ed. MIT Press.

———. 2023. Physically Based Rendering: From Theory to Implementation. 4th ed. MIT Press.

Sabella, Paolo. 1988. “A Rendering Algorithm for Visualizing 3D Scalar Fields.” ACM SIGGRAPH Computer Graphics 22 (4): 51–58.

Saunderson, JL. 1942. “Calculation of the Color of Pigmented Plastics.” JOSA 32 (12): 727–36.

Schweiger, Martin, SR Arridge, M Hiraoka, and DT Delpy. 1995. “The Finite Element Method for the Propagation of Light in Scattering Media: Boundary and Source Conditions.” Medical Physics 22 (11): 1779–92.

Sharma, Gaurav. 2003. “Color Fundamentals for Digital Imaging.” In Digital Color Imaging Handbook, edited by Gaurav Sharma, 14–127. CRC Press.

Sjschen. 2025. “Volume rendered CT scan of a forearm with different color schemes for muscle, fat, bone, and blood; released into the public domain by the copyright holder.” https://commons.wikimedia.org/wiki/File:CTWristImage.png.

Smith, Alvy Ray. 1995. “Alpha and the History of Digital Compositing.” Citeseer.

Stam, Jos. 1995. “Multiple Scattering as a Diffusion Process.” In Rendering Techniques’ 95: Proceedings of the Eurographics Workshop in Dublin, Ireland, June 12–14, 1995 6, 41–50. Springer.

Volz, Hans G, and Frederick T Simon. 2001. Industrial Color Testing. Vol. 2. Wiley-VCH New York.

Westover, Lee. 1990. “Footprint Evaluation for Volume Rendering.” In Proceedings of the 17th Annual Conference on Computer Graphics and Interactive Techniques, 367–76.

Williams, Peter L, and Nelson Max. 1992. “A Volume Density Optical Model.” In Proceedings of the 1992 Workshop on Volume Visualization, 61–68.

Wrenninge, Magnus, and Nafees Bin Zafar. 2011. “Production Volume Rendering: Siggraph 2011 Course.” In ACM SIGGRAPH 2011 Courses, 1–71.

Yifan, Wang, Felice Serena, Shihao Wu, Cengiz Öztireli, and Olga Sorkine-Hornung. 2019. “Differentiable Surface Splatting for Point-Based Geometry Processing.” ACM Transactions On Graphics (TOG) 38 (6): 1–14.

Zwicker, Matthias, Hanspeter Pfister, Jeroen Van Baar, and Markus Gross. 2001a. “EWA Volume Splatting.” In Proceedings Visualization, 2001. VIS’01., 29–538. IEEE.

———. 2001b. “Surface Splatting.” In Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques, 371–78.

---

1. “Isotropic” is a very overloaded term; it just means some physical property is invariant when measured from different directions. So depending on what physical property you care about, “isotropic” can mean different things. The property we care about here is a volume’s ability to absorb photons, which is different from our earlier use of isotropy, which is concerned with the ability of a surface to scatter photons.

2. Technically photons from other objects can enter our eye through a single-scattering event; see the discussion at the end.

3. A direction \(\os\) has a solid angle of 0, so its associated irradiance is technically 0, too. What \(f_p(\os, \oi)L(\oi)\doi\) represents is the irradiance per solid angle.

4. or when the medium consists of randomly distributed and oriented spherically asymmetric particles, in which case the medium is statistically spherically symmetric.

5. In theory, it is certainly possible to have a medium whose total scattering coefficient/efficiency varies with the incident direction but not the angular distribution/probability of the scattered photons.

6. The modern form of the solution is summarized, not invented, by Gustav Mie but the solution had been developed by many predecessors such as Ludvig Lorenz.

7. worked out by Lord Rayleigh, who won the Nobel Prize in Physics in 1904

8. Subrahmanyan Chandrasekhar won the Nobel Prize in physics in 1983 (not for the RTE).

9. Technically, even single scattering can lead to augmentation if there is illumination coming from anywhere outside the ray direction.

10. There are two alternative parameterizations, both of which are common in graphics literature. The first (Pharr, Jakob, and Humphreys 2023) is to express \(p_0 = p + s\os\) (\(s\) being positive), but then the initial radiance would have to be expressed as \(L(p_0, -\os)\), since \(\os\) now points from \(p\) to \(p_0\). The other is to express \(p_0 = p - s\os\) (\(s\) again being positive) (Fong et al. 2017); this avoids the need to switch directions but uses a negative sign. It is a matter of taste which one to use, but be alert to the different conventions.

11. Some (color) transfer functions could have physical underpinnings, such as applying a single-scattering shading algorithms (i.e., local illumination); see, e.g., Levoy (1988, Sect. 3) or Max (1995, Sect. 5), but the goal there is not to precisely model physics but for better, subjective visualization.

12. assuming a diffuse surface so we care to reconstruct the color of each point, not the radiance of each ray.

13. The compound filter combining reconstruction and anti-aliasing filters is sometimes also called a resampling filter, because the compound filter is used during resampling to calculate the new sample values.

14. which says that there exists \(c\in [a, b]\) such that \(\int_a^b f(x)g(x)\d x = f(c)\int_a^b g(x)\d x\), if \(g(x)\) is integrable and does not change its sign in \([a, b]\).

15. If \(L(x, \omega)\) is invariant to \(\omega\) and \(f_p(x, \omega', \hat{\omega})\), under the isotropic medium assumption, is reduced to \(f_p(x, \theta)\), where \(\theta\) is the angle subtended by \(\omega\) and \(\omega'\), both \(S_{\uparrow\downarrow}(x)\) and \(S_{\downarrow\uparrow}(x)\) are essentially calculating \(\int^{\pi}\int^{\pi}c~\d \theta\d \theta'\).

16. You can see that \(T(X)\) does not involve \(R_g\): when we calculate the transmittance of the material, we assume that there is no substrate.