Light Transport and Stochastic Q-Processes

This is the first in a series of tutorials about the mathematical underpinnings of light transport. This tutorial covers, among other things:

• The derivation of the radiative transport differential equation.
• The corresponding Lagrangian stochastic $Q$-process.
• Scattering statistics and mean-free path.

This tutorial is aimed at students at the graduate level or advanced undergraduates.

The Photon

Let's first start by thinking about the simplest object involved in light transfer: the photon. If you're interested in what a photon actually is - an excitation of a quantum field - you'll have to learn quantum field theory, but for our purposes, imagine that a photon is just a discrete particle of light travelling through space. A photon has a few elementary properties, including:

• Position: At any given time $t$, every photon has a position $X_t \in \mathbb{R}^3$.

• Direction: At any given time $t$, every photon has a direction $\omega_t \in \mathbb{S}^2$ (here $\mathbb{S}^2 \subset \mathbb{R}^3$ denotes the 2-sphere of points unit $\ell_2$ distance away from the origin). While direction may vary, the speed of a photon in free space is always constant as $c \approx 3 \cdot 10^8 \, m/s$.

• Wavelength: A photon has a specific wavelength $\lambda>0$. The electromagnetic field carrying a photon usually has a periodicity and $\lambda$ is the length of a single period.

• Frequency: A photon has a frequency $\nu$. The frequency measures the time periodicity at which the underlying electromagnetic field - locally, the oscillations of electromagnetic field returns to its original configuration every $1/\nu$ seconds. The frequency and wavelength are related by the speed of the photon,

• Energy: A photon carries with it a certain amount of energy. This energy comes from the light source that generate the photon via spontaneous or stimulated emission. The Planck-Einstein relation postulates that the energy a photon carries is directly porpotional to its frequency. The constant of porportionality is called the Planck constant $h$:

• Momentum: As postulated by de Broglie, the momentum $p$ of a particle is inversely porpotional to its wavelength, i.e.,

For our purposes the important parts of a photon are really only its position, direction, and wavelength. In the coming sections, we will therefore focus on $(X_t,\omega_t,\lambda)$.

Photon Motion in Free Space

Photons don't just sit there, they tend to move around. In fact, they move around at the speed of light (which is pretty fast). Therefore, if we're describing the motion of a single photon $(X_t,\omega_t,\lambda)$ in free space, it should satisfy the basic equation of motion:

In linear optics, photons tend not to interact with each other, so if we have a collection of photons $(X_t^{(i)},\omega_t^{(i)},\lambda^{(i)})$, the equations are motions are not coupled in any way:

Note that to complete the equations, we need to specify the directions $\omega_t^{(i)}$. In free space, once a photon has been released, it will retain its direction and speed as long as it doesn't hit anything. Therefore, the free space equations for photon motion are exceedingly simple, as the directions $\omega_t^{(i)}$ are all constants $\omega^{(i)}$:

Starting from initial positions $X_0^{(i)}$, these ODEs are very easy to solve:

Note that this analysis is in terms of individual non-interacting particles $X_t^{(i)}$. This is usually refered to as a Lagrangian formulation. Most of the time however, the motion of individual particles is not the end goal. Instead, we are often interested the distribution of particles $\rho_t$. This distribution $\rho_t(x, \omega)$ is specified as a function of both a position in space $x \in \mathbb{R}^3$ and a direction $\omega \in \mathbb{S}^2$. Note that the evolution of $\rho_t(x,ω)$ is not tied to the viewpoint of the individual particles, but rather to the viewpoint of individual locations in free space. This distribution-based view is referred to as an Eulerian formulation.

How does one go from the Lagrangian formulation to an Eulerian formulation? One way is via a technique known as characteristic analysis. Note that the equations of motion are completely reversible - if I specify the position and orientation of a photon at time $t$, you can calculate backwards the position and orientation at time $0$. Therefore, all particles follow non-intersecting characteristic paths $(X_t, \omega_t)$ through the phase space $\mathbb{R}^3 \times \mathbb{S}^2$. This means that the distribution of particles along these paths should remain constant over time. This gives the equation:

Using the chain rule,

In the particular case of free space, $\omega_t$ is constant, therefore, using the equation (6) for $dX_t/dt$ from above, we get

Rearranged this gives

If we fix $\omega$, you might recognize this as the transport equation from partial differential equations.

Scattering in Homogenous Media

Of course, the above equations are very simplistic. They simply say that photons move in straight lines when in free space. Though this is useful, often-times photons are moving through a medium where they constantly collide with other particles and are scattered and absorbed throughout. The different ways materials scatter light are wholy reasonable for the different ways they look. For example, in the above figure, light enters the lion and then scatters around inside it until it exists and hits the viewer's eye. This gives the impression of soft lighting and partial translucency (see how the tail lets some light through with its red-orange highlights). So, if we are going to talk about light transport with an eye towards rendering things like the statue above, we are going to need to discuss how light changes direction.

Unfortunately, this process is not deterministic. Or at least, it is not deterministic from the macroscopic viewpoint. Two light particles can hit the same surface and be reflected in completely different directions. We model this by saying that a light particle hitting a material (or fluid) at location $x$ in direction $\omega_i$ has a certain probability of being redirected in direction ωoωo. It also has a certain probability of being absorbed. In the most general case, we want to think about rates rather than concrete probabilities. When a photon enters a scattering medium like the one pictured above, the photon will move a certain distance until it impacts something in the medium that redirects it. You should think of the contents of the medium as being random, and hence the distance travelled before a photon impacts an obstacle will be random as well. The average distance travelled before a photon hits an obstacle (i.e. scatters) is called the mean free path. The higher the mean free path, the more the photon can travel through the medium unobstructed and therefore the clearer the medium.

Let us take this information and think about the parameters of a particle $(X_t,\omega_t)$ travels through a homogenous scattering medium. Occasionally $\omega_t$ will suddenly change randomly when the particle hits something. But how do we characterize this randomness?

We can model this scattering process by using a continuous $Q$-process on the sphere $\mathbb{S}^2$. Recall that a $Q$-process $W_t$ on $\mathbb{S}^2$ is a continous time random process whose value changes abruptly at random times $T_1,T_2,...,T_n,...$. In particular, a $Q$-process on $\mathbb{S}^2$ is characterized by a generator operator $\mathcal{L}$ which satisfies

where

Here $d\Omega$ denotes the surface area measure on $\mathbb{S}^2$. We'll talk more about $\mathcal{L}$ later, but for now, let's focus on $\sigma(\omega_i,\omega_o)$. You should think of the function $\sigma(\omega_i,\omega_o)$ as the probability transition rates to go from a state $x \rightarrow y$. That is, $W_t$ satisfies:

In other words, after a very small amount of time $dt$, $W_t$ has a probability of about $\sigma(\omega_i,\omega_o) \, d\Omega dt$ to jump to a small solid angle $d\Omega$ around $\omega_o$. This function $\sigma(\omega_i,\omega_o)$ is called the phase function and determines the rate at which a photon travelling in direction $\omega_i$ is redirected to direction $\omega_o$ after it hits an obstacle.

How do we go from a Lagrangian formulation of a $Q$-process to an Eulerian formulation? Well let us ignore the position of our particles (the medium is homogenous, so position doesn't really matter when it comes to rotation) and simply examine the probability distribution of the directions $\omega^{(i)}$ of our particles. We write the distribution of directions of our particles as $\rho_t(\omega)$. You should think of this function as a representing the fact that there are $\rho_t(\omega)$ photons anywhere in the medium that are in direction $\omega$.

Let us compute how the density $\rho_t(\omega)$ changes over time at a given orientation $\omega$. From the definition of the $Q$-process above, we have that

This expression is saying that the number of particles in direction $\omega$ at time $t+dt$ is the number of particles in direction $\omega$ at time $t$, plus the number of particles that are redirected from any other direction $\omega_i$ into direction $\omega$ over the time $dt$, minus the number of particles that are redirected from direction $\omega$ to any other direction $\omega_o$ over the time $dt$. Note that we can take $\rho_t(\omega)$ out of the integral above to get

Rearranging the above and taking $dt\rightarrow0$ gives us

This often appreviated by writing

where $\mathcal{L}^*$ is the adjoint of $\mathcal{L}$ (i.e., $\mathcal{L}^*(\omega_i, \omega_o) = \mathcal{L}(\omega_o, \omega_i)$). This equation is called the forward equation for a $Q$-process. Here, one says that the operator $\mathcal{L}^*$ generates the time evolution of $\rho_t$. This is very similar to deriving the forward equations for a random walk by using Ito's lemma.

The Statistics of Scattering Events

The next question one might ask is for the actual statistics of the scattering times $T_1, T_2, ..., T_n$. These are relatively straightforward. Note that the $Q$-process is Markovian and hence the intervals $\Delta T_i =T_i-T_{i-1}$ are independent and identically distributed. This makes sense, because once the particle has scattered, it doesn't care about anything other than its current position and orientation. In other words, the path a particle has taken to get to a certain point is irrelevant, it is only the current state that matters. To derive statistics for the $\Delta T_i$, we need to consider the distribution $\rho_t^{(0)}(\omega)$. $\rho_t^{(0)}(\omega)$ is defined as the distribution of particles that have not yet scattered. At any given time the probability that a particle scatters in the next $dt$ seconds is given by:

Let $N[t,s]$ be the number of scattering events between times $t$ and $s$. We can determine the distribution of the scattering time $T$ as follows (we simply write $T$ here for a random variable with the distribution of $\Delta T_i$):

For disjoint intervals $A$ and $B$, a property of the $Q$-process is that $N(A)$ and $N(B)$ are independent. Therefore, we have that

We now use the fact that $\mathbb{P}(N[t,t+dt]=1)=\mu_s(\omega)dt+o(dt)$ and $\mathbb{P}(N[t,t+dt]=0)=1-\mu_s(\omega)dt+o(dt)$ to get:

Recall that

Therefore, we get

This is the exponential distribution with rate $\lambda=\mu_s(\omega)$, hence the times between scattering events have distribution

A nice outcome of this fact is that we can easily compute the mean scattering time to be:

The mean free path is therefore

If the medium is isotropic (i.e., the function $\sigma(\omega_i, \omega_o)$ is rotationally invariant), then $\mu_s(\omega)$ does not depend on $\omega$ and this can just be written as

Full Eulerian Formulation

One can obtain a full Eulerian formulation for the evolution of $\rho_t(x,\omega)$ as follows. Let us consider the number of particles in a set $A$ with direction $\omega$, i.e.,

The change in the number of particles in $A$ moving in direction $\omega$ is given by the particle flux through the boundary $\partial A$ as well as the particles inside $A$ that scatter either from some direction $\omega_i$ to $\omega$ or scatter from direction $\omega$ to some direction $\omega_o$. We can write this as

The $cdt \int_{\partial A} \rho_t(x,\omega) (\omega \cdot dn)$ represents the sum of all particles that move out of the region $A$ over a time interval $dt$ and the $dt \int_{A} \mathcal{L}_\omega^∗ \rho_t(x,\omega)\, dx\, dt$ represents the change in particle number from scattering inside the region $A$ (as we derived two sections ago). Stokes' Theorem can be applied to the term over $\partial A$ to get

Since this is true for any $A$, it must be the case that

Now rearranging and letting $dt \rightarrow 0$,

Voila. This is the Eulerian formulation of light transport in a homogenous medium with phase function $\sigma(\omega_i, \omega_o)$. Writing it out in full gives:

For comparison, the Langrangian form of these equations is

where $W_t$ is a $Q$-process on $\mathbb{S}^2$ with generator $\mathcal{L}$.

Isotropic Scattering

The simplest example of a scattering medium is one that is fully isotropic, i.e., regardless of the incomming and outgoing directions $\omega_i$ and $\omega_o$, the rate of scattering is the same. That is,

Usually this is referred to as isotropic scattering. In this case, the scattering coefficient $\mu_s$ does not depend on $\omega$ and is easily calculated to be

The equations in the previous section now take on the form

Here $\bar{\rho}_t$ is shorthand for the average photon density at $x$ over all directions $\omega$,

Scattering in Heterogenous and Absorbing Media

In the previous examples, we assumed that the scattering medium was the same everywhere, i.e., that the phase function $\sigma(\omega_i,\omega_o)$ did not depend on $x$. In this section, we see what happens when the phase function $\sigma(\omega_i,\omega_o)$ does depend on $x$. This is relatively common, many objects are not entirely uniform. For example, your hand scatters light, but your hand is not uniformly made of the same material. There's skin, muscle, bone, blood, etc. And each of these materials have different properties. How do we model this? Well, now the statistics of the direction $\omega_t$ of a particle will depend heavily on that particle's position and path.

If this is the case, then $W_t$ is now a $Q$-process whose phase function $\sigma(\omega_i,\omega_o)$ changes depending on the position XtXt of the particle. We will write the new phase function as $\sigma(x;\omega_i,\omega_o)$. This means that the generator of the process now depends explicitly on the value of $x$. Hence, we have

where once again,

In addition to this, we also give a particle a certain probability of being absorbed at any point in time. This probability is quantified by the absorption coefficient $\mu_a(x,\omega)$ which is defined in a similar way to the phase function:

What this means is that the $Q$-process $W_t$ now has a chance of jump to a dead state from which it can never return. Therefore, the number of particles in a region $A$ now changes over a small period dtdt according to:

The terms with a $dt$ factor represent, in order of appearence:

1. The net number of particles that leave the region $A$.
2. The number of particles that are scattered in $A$ from some direction $\omega_i$ to $\omega$.
3. The number of particles that are scattered in $A$ from $\omega$ to some other direction $\omega_o$.
4. The number of particles absorbed in $A$.

Repeating the procedure of the previous section will give us:

Alternatively, we can just write this as:

Where the generator $\mathcal{L}_{x,\omega}$ of the full joint process $(X_t,W_t)$ is given by:

We have therefore derived the full radiative transport equation.

Statistics of Scattering Events

Let us assume for now that the scattering medium is non-absorbing, i.e., that $\mu_a=0$ everywhere. Just like in the homogenous case, every particle will have a series of scattering times $T_1,T_2,...,T_n,...$. The difference now is that these scattering times depend explicitly on the path the photon is travelling along.

We can calculate the statistics of scattering events by using the same trick as last time. Let $N[t,s]$ be the number of scattering events in the time interval $[t,s]$. Suppose that the particle has direction ωω and is traveling along the ray $X_t(0)=X_0+c\omega t$. The we have that:

When $x \ll 1$, we know that $1-x \sim \exp⁡(-x)$. Therefore,

We then have:

As $dt \rightarrow 0$ this converges to

Therefore, if one wishes to simulate a photon scattering in a heterogenous medium, one might sample scattering times from the above distribution.

Statistics of Absorption Events

The above section only covers scattering events. However, a photon may experience either a scattering or an absorption event while it is traversing through a medium. Recall, that photons travelling in direction $\omega$ at a position $x$ are absorbed at a rate $\mu_a(\omega)$. Combining this with scattering events, at every time step of $dt$, a photon has a probability $\mu_a(\omega)\, dt$ of being absorbed and a $\mu_s(\omega)\, dt$ probability of being scattered. Let $A$ be the event where the photon is absorbed. Therefore, one can follow the derivation of the previous section to obtain:

Therefore, the probability that the photon is absorbed is given by:

Similarly, the probability that the photon is scattered is given by:

If the media is homogenous, that is if $\mu_s(x,\omega)=\mu_s(\omega)$ and $\mu_a(x,\omega)= \mu_a(\omega)$, this reduces to:

Emissive Media

In addition to absorbing and scattering light, a general medium can also emit photons. To quantify this emission, we introduce an emission rate $j(x,\omega)$. The quantity of light emitted at $x$ in direction $\omega$ over a short period $dt$ is given by

Therefore, it is relatively simple to see that we should update the Eulerian equations by adding a $j(x,\omega)$ term to the right hand side:

Or written out in full:

Lights are examples of emissive materials.

Steady-State Radiative Transfer

Often times, we aren't interested in the evolution of the distribution $\rho_t(x,\omega)$ over a short time scale, but rather, we're interested in the average distribution over a very long time. Our eyes act like integrators, so if we want to talk about rendering, we care about the time averaged distribution of photons,

This time-averaged distribution should represent a sort of "steady state" of the system when a given $j(x,\omega,t)$ is applied as an input. Computing this also requires time integrating out fluctuations in the emissive term $j(x,\omega,t)$,

At equilibrium, we expect that the distribution of light no longer changes, so an easy way to compute it is to set

This can be written out as:

These are the elliptic equations for the steady state of a radiative transfer system. Note that the equations are linear in $\rho$. Therefore, we know that the function ρρ must be linear in $\bar{j}$ ,. This linearity can be written in terms of a scattering operator $\mathcal{S}$ that acts on the function $\bar{j}$. This function corresponds to applying the Green's function $\mathcal{G}$ of the elliptic equation above:

Steady-State Transfer with Boundary Conditions

Often times we're not interested in scattering on the whole of $\mathbb{R}^3$, but rather on some restricted domain $D$. We can easily imagine the following question: given the geometry and makeup of someone's face, suppose we put a flashlight under it and want to know what the face will look like when after the light propagates through it.

From such a problem to be well posed, we need to specify an ingoing light at the boundary. We will denote this as $f(x,\omega)$ (for "flashlight"). $f(x,\omega)$ is defined on the boundary $x \in \partial D$ for ingoing directions $\omega \cdot n < 0$. Furthermore, we may have emissive materials inside the domain $D$, which we will have to model with the function $\bar{j}(x,\omega)$ defined on the interior of $D$. Therefore, we have

And we have boundary conditions

In the same way as before, the output $\rho$ is linear in the inputs $f$ and $\bar{j}$. It can therefore be written as:

where $\mathcal{S}_{\partial D}$ and $\mathcal{S}_{\text{int}(D)}$ are scattering operators (from boundary $\partial D$ to $D$ and interior $\text{int}(D)$ to $D$ respectively). They also have corresponding Green's functions:

The Green's function $\mathcal{G}_{\partial D}(x,\omega, x', \omega')$ when the output is restricted to the boundary $x \in \partial D$ and $\omega' \cdot n(x)<0$ is called the bidirectional scattering distribution function (BSDF) of the object $D$.

Radiant Flux, Irradiance, and Radiance

There are some fundamental quantities that often get tossed around in light transport. It is good to put everything together with standard temonology. Let us go through all the necessary terms one-by-one:

1. Radiant Flux (across a surface $S$): The radiant flux through $S$ is given by the energy that moves through $S$ per unit time. Recall that the energy per particle is given by

Therefore, for monochromatic light with wavelength $\lambda$, the radiance flux across $S$ is given by:

2. Irradiance: The irradiance at a location $x$ on a surface $S$ is a measure of the energy recieved/emitted by the surface per unit time per unit area. It is written as:

3. Radiance: The radiance at a location $x$ and in direction $\omega$ on a surface is a measure of the energy recieved/emitted by the surface per unit time per unit area per solid angle. It is often written:

It is common to write $\omega\cdot n(x)$ as $\cos \theta$ where $\theta$ is the angle between the direction of the photon $\omega$ and the normal $n(x)$ of the surface,