Volume Rendering
Introduction
Volume rendering is a technique used to display a 2D projection of a 3D scalar field. It’s widely used for 3D data visualization, such as medical imaging.
In this blog post, we will discuss mathematically the basic optical model for volume rendering and how to estimate the observed color of a ray casting on the 2D image.
Volume Rendering
The process of volume rendering maps a 3D scalar field of optical properties to a 2D image. The visual appearance on the 2D image is computed by integrating the optical properties along the ray.
In the basic optical model, each point in the volume is an infinitesimal particle with a certain opacity $\tau$ that emits the light represented by a color $c$.
The probability of a ray $\mathbf{r}$ hitting a particle along an infinitesimal interval $ds$ is $\tau(\mathbf{r}(s)) ds$, where $\mathbf{r}(s) = \mathbf{o} + s \mathbf{d}$ is the point on the ray $\mathbf{r}$ at distance $s$ from the origin $\mathbf{o}$ in the direction $\mathbf{d}$.
The transmittance of the ray $\mathbf{r}$ from the origin $\mathbf{o}$ to the point $\mathbf{r}(s)$, $T_{\mathbf{r}}(s)$, is defined as the probability of the ray traveling from the origin to the point $\mathbf{r}(s)$ without hitting any particle.
$$
\begin{align}
T_{\mathbf{r}}(s + ds) &= T_{\mathbf{r}}(s) (1 - \tau(\mathbf{r}(s)) ds) \\
\end{align}
$$
The infinitesimal change in the transmittance of the ray $\mathbf{r}$ at distance $s$ from the origin $\mathbf{o}$ is
$$
\begin{align}
d T_{\mathbf{r}}(s)
&= T_{\mathbf{r}}(s + ds) - T_{\mathbf{r}}(s) \\
&= - T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) ds \\
\end{align}
$$
$$
\begin{align}
\frac{d T_{\mathbf{r}}(s)}{T_{\mathbf{r}}(s)}
&= - \tau(\mathbf{r}(s)) ds \\
\end{align}
$$
This is an ordinary differential equation. To solve this, the definite integral could be applied to both sides of the equation to get the transmittance of the ray $\mathbf{r}$ from the origin $\mathbf{o}$ to the point $\mathbf{r}(s)$.
$$
\begin{align}
\int_{0}^{s} \frac{d T_{\mathbf{r}}(s)}{T_{\mathbf{r}}(s)}
&= - \int_{0}^{s} \tau(\mathbf{r}(u)) du \\
\ln T_{\mathbf{r}}(s) - \ln T_{\mathbf{r}}(0)
&= - \int_{0}^{s} \tau(\mathbf{r}(u)) du \\
\end{align}
$$
Note that because $T_{\mathbf{r}}(0) = 1$,
$$
\begin{align}
\ln T_{\mathbf{r}}(s)
&= - \int_{0}^{s} \tau(\mathbf{r}(u)) du \\
T_{\mathbf{r}}(s)
&= \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) \\
\end{align}
$$
The probability density of the ray $\mathbf{r}$ terminating exactly at the point $\mathbf{r}(s)$ is
$$
\begin{align}
P_{\mathbf{r}}(s)
&= T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) \\
\end{align}
$$
The observed color of the ray $\mathbf{r}$ on the 2D image is the expected value of the color $c(\mathbf{r}(s))$ of the particle along the ray $\mathbf{r}$ at distance $s$ from the origin $\mathbf{o}$.
$$
\begin{align}
\mathbb{E}_{s \sim P_{\mathbf{r}}(s)}[c(\mathbf{r}(s))]
&= \int_{0}^{\infty} P_{\mathbf{r}}(s) c(\mathbf{r}(s)) ds \\
&= \int_{0}^{\infty} T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
&= \int_{0}^{\infty} \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
\end{align}
$$
The integral in $T_{\mathbf{r}}(s)$ is usually not analytically solvable. However, it could be approximated by Riemann sum.
$$
\begin{align}
\int_{0}^{s} \tau(\mathbf{r}(u)) du &\approx \sum_{i = 0}^{n - 1} \tau(\mathbf{r}(s_i)) \Delta s_i \\
\end{align}
$$
where $\Delta s_i = s_{i + 1} - s_i$ and $s_n = s$. Thus,
$$
\begin{align}
T_{\mathbf{r}}(s)
&= \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) \\
&\approx \exp\left(- \sum_{i = 0}^{n - 1} \tau(\mathbf{r}(s_i)) \Delta s_i\right) \\
&= \exp\left( - \sum_{i = 0}^{n - 1} \tau(\mathbf{r}(s_i)) \left(s_{i + 1} - s_i\right)\right) \\
\end{align}
$$
$$
\begin{align}
T_{\mathbf{r}}(s_k)
&\approx \exp\left( - \sum_{i = 0}^{k - 1} \tau(\mathbf{r}(s_i)) \left(s_{i + 1} - s_i\right)\right) \\
\end{align}
$$
The observed color of the ray $\mathbf{r}$ on the 2D image is usually integrated over a finite interval $[0, s]$ rather than an infinite interval $[0, \infty)$.
Assuming for $\forall s \in [s_i, s_{i + 1}]$, $\tau(\mathbf{r}(s))$ and $c(\mathbf{r}(s))$ are constant, the observed color of the ray $\mathbf{r}$ on the 2D image approximately becomes
$$
\begin{align}
\mathbb{E}_{s \sim P_{\mathbf{r}}(s)}[c(\mathbf{r}(s))]
&= \int_{0}^{\infty} T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
&= \int_{0}^{s} T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
&= \sum_{i = 0}^{n-1} \int_{s_i}^{s_{i + 1}} T_{\mathbf{r}}(s) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
&= \sum_{i = 0}^{n-1} \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) \tau(\mathbf{r}(s)) c(\mathbf{r}(s)) ds \\
&\approx \sum_{i = 0}^{n-1} \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{0}^{s} \tau(\mathbf{r}(u)) du\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{0}^{s_i} \tau(\mathbf{r}(u)) du\right) \exp\left(- \int_{s_i}^{s} \tau(\mathbf{r}(u)) du\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) \int_{s_i}^{s_{i + 1}} T_{\mathbf{r}}(s_i) \exp\left(- \int_{s_i}^{s} \tau(\mathbf{r}(u)) du\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{s_i}^{s} \tau(\mathbf{r}(u)) du\right) ds \\
&\approx \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \int_{s_i}^{s_{i + 1}} \exp\left(- \int_{s_i}^{s} \tau(\mathbf{r}(s_i)) du\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \int_{s_i}^{s_{i + 1}} \exp\left(- \tau(\mathbf{r}(s_i)) \int_{s_i}^{s} du\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \int_{s_i}^{s_{i + 1}} \exp\left(- \tau(\mathbf{r}(s_i)) (s - s_i)\right) ds \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \left( \frac{1}{-\tau(\mathbf{r}(s_i))} \exp\left(- \tau(\mathbf{r}(s_i)) (s - s_i)\right) \right) \Bigg|_{s_i}^{s_{i + 1}} \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \left( \frac{1}{\tau(\mathbf{r}(s_i))} \exp\left(- \tau(\mathbf{r}(s_i)) (s_{i + 1} - s_i)\right) - \frac{1}{-\tau(\mathbf{r}(s_i))} \exp\left(- \tau(\mathbf{r}(s_i)) (s_i - s_i)\right) \right) \\
&= \sum_{i = 0}^{n-1} \tau(\mathbf{r}(s_i)) c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \left( \frac{1}{-\tau(\mathbf{r}(s_i))} \exp\left(- \tau(\mathbf{r}(s_i)) (s_{i + 1} - s_i)\right) - \frac{1}{-\tau(\mathbf{r}(s_i))} \right) \\
&= \sum_{i = 0}^{n-1} c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \left( 1 - \exp\left(- \tau(\mathbf{r}(s_i)) (s_{i + 1} - s_i)\right) \right) \\
&= \sum_{i = 0}^{n-1} c(\mathbf{r}(s_i)) T_{\mathbf{r}}(s_i) \left( 1 - \exp\left(- \tau(\mathbf{r}(s_i)) \Delta s_i\right) \right) \\
\end{align}
$$
where $T_{\mathbf{r}}(s_i)$ for $\forall i \in [0, n - 1]$ could be computed by the recursive formula relative efficiently.
$$
\begin{align}
T_{\mathbf{r}}(s_i)
&\approx \exp\left(- \sum_{j = 0}^{i - 1} \tau(\mathbf{r}(s_j)) \left(s_{j + 1} - s_j\right)\right) \\
\end{align}
$$
This means that the observed color of the ray $\mathbf{r}$ on the 2D image could be estimated by quadrature and the computation could be parallelized.
Conclusions
Based on the basic optical model and some approximations, we could easily estimate the observed color of a ray casting on the 2D image using the opacity $\tau$ and the color $c$ of the infinitesimal particles in the volume and this rendering process is fully differentiable with respect to the opacity $\tau$ and the color $c$ of the infinitesimal particles. In this way, by reverse rendering, we could estimate the opacity $\tau$ and the color $c$ of the infinitesimal particles in the volume from the observed color of the ray casting on multiple 2D images of the volume. This is also called the 3D reconstruction from 2D images.
References
Volume Rendering