Light Transport Equation

Basic Form

Previously, we have seen the reflection equation:

$$L_o(\mathbf{p}, \boldsymbol{\omega}_o) = \int_{H^2} f_{\mathrm{r}} \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right) L_i(\mathbf{p}, \boldsymbol{\omega}_i) \cos \theta_i \, d\Omega_i.$$

If we replace the bidirectional reflectance distribution function $f_{\mathrm{r}} \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right)$, which describes surface reflection, with the bidirectional scattering distribution function $f \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right)$, and add the emitted radiance $L_e(\mathbf{p}, \boldsymbol{\omega}_o)$, we obtain the scattering equation with the emission term:

$$L_o(\mathbf{p}, \boldsymbol{\omega}_o) = L_e(\mathbf{p}, \boldsymbol{\omega}_o) + \int_{S^2} f \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right) L_i(\mathbf{p}, \boldsymbol{\omega}_i) \left|\cos \theta_i\right| d\Omega_i.$$

Notice that in the above equation, $L_e(\mathbf{p}, \boldsymbol{\omega}_o)$ and $f \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right)$ are known quantities. If we can express the function $L_i$ in terms of $L_o$, we obtain an integral equation involving only $L_o$. To achieve this, we introduce the ray-casting function $t(\mathbf{p}, \boldsymbol{\omega})$. Suppose a ray is cast from point $\mathbf{p}$ in direction $\boldsymbol{\omega}$; the point where the ray first intersects a surface in the scene is denoted by $t(\mathbf{p}, \boldsymbol{\omega})$. Thus, we have

$$L_i(\mathbf{p}, \boldsymbol{\omega}_i) = L_o(t(\mathbf{p}, \boldsymbol{\omega}_i), -\boldsymbol{\omega}_i).$$

To handle the special case where the ray does not intersect any surface in the scene, we define $t(\mathbf{p}, \boldsymbol{\omega})$ to take a special value $\Lambda$, and for any $\boldsymbol{\omega}$, we have $L_o(\Lambda, \boldsymbol{\omega}) = 0$.

Using the above relation and abbreviating $L_o(\mathbf{p}, \boldsymbol{\omega}_o)$ as $L(\mathbf{p}, \boldsymbol{\omega}_o)$, the resulting equation is known as the light transport equation:

$$L(\mathbf{p}, \boldsymbol{\omega}_o) = L_e(\mathbf{p}, \boldsymbol{\omega}_o) + \int_{S^2} f \left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i \right) L(t(\mathbf{p}, \boldsymbol{\omega}_i), -\boldsymbol{\omega}_i) \left|\cos \theta_i\right| d\Omega_i.$$

The light transport equation is also known as the rendering equation. It is an integral equation concerning the function $L(\mathbf{p}, \boldsymbol{\omega}_o)$, which describes the radiance at any point on a surface in any direction under equilibrium. In fact, we can also obtain the radiance at any point in the scene in any direction because, using the previously derived relation:

$$L(\mathbf{p}, \boldsymbol{\omega}) = L(t(\mathbf{p}, \boldsymbol{\omega}), -\boldsymbol{\omega}),$$

we can convert the radiance in the scene to the radiance on the surfaces of objects.

Surface Form

In the above light transport equation, the geometric information of objects in the scene is entirely encapsulated in the ray-casting function $t(\mathbf{p}, \boldsymbol{\omega})$. To make this information explicit, we will derive the surface form of the light transport equation.

Previously, when integrating the incident light from various directions at a point, we performed the integration on the unit sphere. Considering that the incident light received at a point actually comes from the emitted light at various points on other surfaces, we can convert the integration on the unit sphere to an integration on the surfaces of other objects. If we want to estimate this integral using the Monte Carlo integration method, our sampling region will shift from the unit sphere to the surfaces of other objects.

Next, we introduce some notation. If a ray from point $\mathbf{p}’$ in direction $\boldsymbol{\omega}$ passes through point $\mathbf{p}$, we denote it by

$$L(\mathbf{p}' \to \mathbf{p}) := L(\mathbf{p}', \boldsymbol{\omega}) = L \left(\mathbf{p}', \frac{\mathbf{p} - \mathbf{p}'}{\|\mathbf{p} - \mathbf{p}'\|} \right).$$

As shown in the figure below,

Figure 1 - Illustration of $f(\mathbf{p}^{\prime\prime}\to\mathbf{p}^\prime\to\mathbf{p})$

If $t(\mathbf{p}^\prime, \boldsymbol{\omega}_o) = \mathbf{p}$ and $t(\mathbf{p}^\prime, \boldsymbol{\omega}_i) = \mathbf{p}^{\prime\prime}$, we denote

$$f(\mathbf{p}^{\prime\prime} \to \mathbf{p}^\prime \to \mathbf{p}) := f\left(\mathbf{p}, \boldsymbol{\omega}_o, \boldsymbol{\omega}_i\right) = f\left(\mathbf{p}^\prime, \frac{\mathbf{p} - \mathbf{p}^\prime}{\Vert \mathbf{p} - \mathbf{p}^\prime \Vert}, \frac{\mathbf{p}^{\prime\prime} - \mathbf{p}^\prime}{\Vert \mathbf{p}^{\prime\prime} - \mathbf{p}^\prime \Vert}\right).$$

Let $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ be the visibility function, which is 1 if there is no occlusion between points $\mathbf{p}$ and $\mathbf{p}^{\prime}$, meaning they are mutually visible, i.e.,

$$t\left(\mathbf{p}, \frac{\mathbf{p}^\prime - \mathbf{p}}{\Vert \mathbf{p}^\prime - \mathbf{p} \Vert}\right) = \mathbf{p}^\prime, \quad t\left(\mathbf{p}^\prime, \frac{\mathbf{p} - \mathbf{p}^\prime}{\Vert \mathbf{p} - \mathbf{p}^\prime \Vert}\right) = \mathbf{p}$$

In this case, $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ equals 1; otherwise, $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right)$ equals 0.

The geometric coupling term is defined as the product of the visibility function and the Jacobian determinant of the integral region transformation, i.e.,

$$G\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right) = V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right) \frac{|\cos \theta| \left|\cos \theta^{\prime}\right|}{\left\|\mathbf{p} - \mathbf{p}^{\prime}\right\|^{2}},$$

where $\theta$ is the angle between the normal vector at point $\mathbf{p}$ and $\mathbf{p}^\prime - \mathbf{p}$, and $\theta^\prime$ is the angle between the normal vector at point $\mathbf{p}^\prime$ and $\mathbf{p} - \mathbf{p}^\prime$.

When $V\left(\mathbf{p} \leftrightarrow \mathbf{p}^{\prime}\right) = 1$, the light transport equation can be rewritten as an integral over the surfaces of all objects in the scene, denoted as $A$:

$$L\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right) = L_{\mathrm{e}}\left(\mathbf{p}^{\prime} \rightarrow \mathbf{p}\right) + \int_{A} f\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime} \rightarrow \mathbf{p}\right) L\left(\mathbf{p}^{\prime \prime} \rightarrow \mathbf{p}^{\prime}\right) G\left(\mathbf{p}^{\prime \prime} \leftrightarrow \mathbf{p}^{\prime}\right) \, dA\left(\mathbf{p}^{\prime \prime}\right).$$

This is the surface form of the light transport equation.

Path Integral Form

The surface form of the light transport equation provides the recursive relationship for $L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right)$ concerning $i$:

$$\begin{aligned} &L\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right) \\ =\;&L_{\mathrm{e}}\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right) + \int_{A} f\left(\mathbf{p}_{i+2} \rightarrow \mathbf{p}_{i+1} \rightarrow \mathbf{p}_i\right) L\left(\mathbf{p}_{i+2} \rightarrow \mathbf{p}_{i+1}\right) G\left(\mathbf{p}_{i+2} \leftrightarrow \mathbf{p}_{i+1}\right) \, dA\left(\mathbf{p}_{i+2}\right). \end{aligned}$$

Starting from $i = 0$, we first express $L\left(\mathbf{p}_2 \rightarrow \mathbf{p}_1\right)$ on the right-hand side of the equation in terms of $L\left(\mathbf{p}_3 \rightarrow \mathbf{p}_2\right)$, then express $L\left(\mathbf{p}_3 \rightarrow \mathbf{p}_2\right)$ on the right-hand side of the equation in terms of $L\left(\mathbf{p}_4 \rightarrow \mathbf{p}_3\right)$, and so on. Eventually, we obtain the following equation:

$$\begin{aligned} L\left(\mathbf{p}_1 \rightarrow \mathbf{p}_0\right) &= L_{\mathrm{e}} \left(\mathbf{p}_1 \rightarrow \mathbf{p}_0\right) \\ &\quad + \int_{A} L_{\mathrm{e}}\left(\mathbf{p}_2 \rightarrow \mathbf{p}_1\right) f\left(\mathbf{p}_2 \rightarrow \mathbf{p}_1 \rightarrow \mathbf{p}_0\right) G\left(\mathbf{p}_2 \leftrightarrow \mathbf{p}_1\right) \, dA\left(\mathbf{p}_2\right) \\ &\quad + \int_{A} \int_{A} L_{\mathrm{e}}\left(\mathbf{p}_3 \rightarrow \mathbf{p}_2\right) f\left(\mathbf{p}_3 \rightarrow \mathbf{p}_2 \rightarrow \mathbf{p}_1\right) G\left(\mathbf{p}_3 \leftrightarrow \mathbf{p}_2\right) \\ &\quad \quad \times f\left(\mathbf{p}_2 \rightarrow \mathbf{p}_1 \rightarrow \mathbf{p}_0\right) G\left(\mathbf{p}_2 \leftrightarrow \mathbf{p}_1\right) \, dA\left(\mathbf{p}_3\right) dA\left(\mathbf{p}_2\right) \\ &\quad + \cdots \end{aligned}$$

Let

$$T\left(\overline{\mathbf{p}}_n\right) := \prod_{i=1}^{n-1} f\left(\mathbf{p}_{i+1} \rightarrow \mathbf{p}_i \rightarrow \mathbf{p}_{i-1}\right) G\left(\mathbf{p}_{i+1} \leftrightarrow \mathbf{p}_i\right)$$

be the throughput of path $\overline{\mathbf{p}}_n = (\mathbf{p}_n, \mathbf{p}_{n-1}, \cdots, \mathbf{p}_0)$. Also let

we can write $L\left(\mathbf{p}_1 \rightarrow \mathbf{p}_0\right)$ as a series:

$$L\left(\mathbf{p}_1 \rightarrow \mathbf{p}_0\right) = \sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_n\right).$$

This is known as the path integral form of the light transport equation. Notice that the right-hand side of the equation no longer contains $L$, so we have explicitly written the solution to the light transport equation in series form. $P\left(\overline{\mathbf{p}}_1\right)$ represents the light directly emitted from the source to the camera, $P\left(\overline{\mathbf{p}}_2\right)$ represents directly scattered light, and $P\left(\overline{\mathbf{p}}_3\right), P\left(\overline{\mathbf{p}}_4\right), \cdots$ represent multiple scattered indirect light.

Path Tracing Algorithm

The goal of the path tracing algorithm is to estimate

$$L\left(\mathbf{p}_1 \rightarrow \mathbf{p}_0\right) = \sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_n\right).$$

There are two problems to solve: the first is how to estimate the infinite sum, and the second is how to estimate the integral term $P\left(\overline{\mathbf{p}}_n\right)$.

Russian Roulette

The difficulty of the infinite sum can be solved by the so-called Russian roulette method.

Suppose we have an unbiased estimator $\widehat{P}\left(\overline{\mathbf{p}}_n\right)$ for any $P\left(\overline{\mathbf{p}}_n\right)$, then we can immediately write an unbiased estimator for

$$\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_n\right).$$

However, in reality, we can only obtain a finite number of terms for $\widehat{P}\left(\overline{\mathbf{p}}_{n}\right)$. The idea of Russian roulette is to discard the sum of all terms after the $N$th term with a certain probability $q$, i.e., $\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)$. If not discarded, we amplify $\sum\limits_{n=N}^{\infty} \widehat{P} \left(\overline{\mathbf{p}}_{n}\right)$ as compensation.

For the first two terms $\widehat{P}\left(\overline{\mathbf{p}}_{1}\right)$, $\widehat{P}\left(\overline{\mathbf{p}}_{2}\right)$, we do not use Russian roulette but calculate them directly. Starting from the third term, we perform Russian roulette. Let $\eta_1$ follow a 0-1 distribution,

$$\mathrm{P}(\eta_1=k) = q^{k}(1-q)^{1-k}, \quad k=0,1,$$

When $\eta_1$ is 1, we discard the subsequent terms, taking only the sum of the first two terms:

$$\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right),$$

When $\eta_1$ is 0, we scale the subsequent terms by a factor of $\frac{1}{1-q}$, thus obtaining:

$$\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1}{1-q} \sum_{n=3}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right).$$

Combining these two cases, our estimator should be written as:

$$\widehat{Q}_3 = \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1_{\eta_1=0}}{1-q} \sum_{n=3}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right).$$

Noting that

$$\mathrm{E}\left[\frac{1_{\eta_1=0}}{1-q}\right] = \frac{\mathrm{P}(\eta_1=0)}{1-q} = \frac{1-q}{1-q} = 1,$$

we have

$$\mathrm{E}\left[ \sum_{n=1}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) \right] = \mathrm{E}\left[\sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1_{\eta_1=0}}{1-q} \sum_{n=3}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right],$$

i.e., $\widehat{Q}_3$ is also an unbiased estimator of $\sum\limits_{n=1}^{\infty} {P}\left(\overline{\mathbf{p}}_{n}\right)$.

Let $\eta_2$ be independent and identically distributed with $\eta_1$. We can similarly construct

$$\widehat{Q}_4 = \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right) + \frac{1_{\eta_2=0}}{1-q} \sum_{n=4}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right).$$

Clearly, $\mathrm{E}\left[\widehat{Q}_{4}\right] = \mathrm{E}\left[\widehat{Q}_{3}\right]$, so $\widehat{Q}_4$ is also an unbiased estimator.

Generally, if $\eta_1, \eta_2, \cdots$ are independent and all follow a 0-1 distribution with success probability $q$, then

$$\begin{aligned} \widehat{Q}_{N+2} &= \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1_{\eta_1=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{3}\right) + \frac{1_{\eta_2=0}}{1-q}\left(\widehat{P}\left(\overline{\mathbf{p}}_{4}\right) + \cdots + \frac{1_{\eta_N=0}}{1-q}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right)\right)\right)\\ &= \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \frac{1_{\eta_1=0}}{1-q}\widehat{P}\left(\overline{\mathbf{p}}_{3}\right) + \frac{1_{\eta_1,\eta_2=0}}{(1-q)^2}\widehat{P}\left(\overline{\mathbf{p}}_{4}\right) + \cdots + \frac{1_{\eta_1,\eta_2,\cdots,\eta_N=0}}{(1-q)^N}\sum_{n=N+2}^{\infty} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) \end{aligned}$$

is an unbiased estimator of $\sum\limits_{n=1}^{\infty} {P}\left(\overline{\mathbf{p}}_{n}\right)$. Let $\widehat{Q} = \lim\limits_{N \to \infty} \widehat{Q}_N$, then $\widehat{Q}$ is also an unbiased estimator.

Let

$$\tau = \inf \{ i \mid \eta_i = 1 \},$$

i.e., all terms from the $\tau + 2$th term onward are discarded. We can write $\widehat{Q}$ as

$$\widehat{Q} = \sum_{n=1}^{2} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) + \sum_{n=3}^{\infty} \frac{1}{(1-q)^{n-2}} \widehat{P}\left(\overline{\mathbf{p}}_{n}\right) 1_{\tau > n-2}.$$

It is not difficult to find that $\tau$ follows a geometric distribution, with its probability mass function given by

$$\mathrm{P}(\tau = k) = (1-q)^{k-1} q, \quad k = 1, 2, \cdots.$$

Thus, we can use check the unbiasedness of $\widehat{Q}$ as follows:

$$\begin{aligned} \mathrm{E}\left[\widehat{Q}\right]&=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\mathrm{P}(\tau>n-2)\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)\sum_{k=n-1}^{\infty}(1-q)^{k-1}q\\ &=\sum_{n=1}^{2} P\left(\overline{\mathbf{p}}_{n}\right)+\sum_{n=3}^{\infty}\frac{1}{(1-q)^{n-2}} P\left(\overline{\mathbf{p}}_{n}\right)(1-q)^{n-2}\\ &=\sum_{n=1}^{\infty} P\left(\overline{\mathbf{p}}_{n}\right). \end{aligned}$$

Note that $\mathrm{P}(\tau<\infty)=1$, i.e., the probability that $\widehat{Q}$ is a finite sum is $1$. We can totally estimate the value of $\widehat{Q}$ in practice.

Pseudo Code

The pseudo code of the path tracing algorithm is shown below:

shade(p, wo):
    # Contribution from the light source.
    uniformly sample the light at x’ (pdf_light = 1 / A);
    shoot a ray from p to x’;
    L_dir = 0.0;
    if the ray is not blocked in the middle:
        L_dir = L_i * f_r * cos θ * cos θ’ / |x’ - p|^2 / pdf_light;

    # Contribution from other reflectors.
    L_indir = 0.0;
    test Russian Roulette with probability P_RR;
    if passing Russian Roulette test:
        uniformly sample the hemisphere toward wi (pdf_hemi = 1 / 2pi);
        trace a ray r(p, wi);
        if ray r hit a non-emitting object at q:
            L_indir = shade(q, -wi) * f_r * cos θ / pdf_hemi / P_RR;
    
    return L_dir + L_indir;