3.1.

Forward Model

In the forward problem, as we mentioned before, we use a path integral to build a mathematical model for the light transport. Here, we follow the notation developed in the computer graphics literature^17,23,53,60 to introduce the path integral. Sections 3.2 to 3.6 will show the simplified model we propose.

Given a space ${\mathfrak{R}}^3$, a light source is located at $x_0\in {\mathfrak{R}}^3$ and a detector at $x_{M+1}\in {\mathfrak{R}}^3$, and in between them is the participating media $\nu \subset {\mathfrak{R}}^3$ with boundary $\partial \nu$ and interior volume ${\nu }_0≔\nu \setminus \partial \nu$. A light path $\tilde{x}$ connecting $x_0$ and $x_{M+1}$ of length $M+2$ consists of $M+2$ vertices $x_m\in {\mathfrak{R}}^3$ for $m=0,1,\dots ,M+1$, denoted by $\tilde{x}=x_0,x_1,\cdots ,x_M,x_{M+1}$. Thus, absorption, scattering, or reflection events happen at $x_1,\dots ,x_M$. The set of all paths of length $M$ is denoted by ${\mathrm{\Omega }}_M$. The path space $\mathrm{\Omega }$ is the countable set of all paths ${\mathrm{\Omega }}_M$ of finite length.

A direction is denoted by $\omega \in S^2$, where $S^2$ is a unit sphere in ${\mathfrak{R}}^3$. A unit vector ${\omega }_{x_m,x_{m+1}}$ is the direction from vertex $x_m$ to vertex $x_{m+1}$ in a path $\tilde{x}$.

Veach²⁰ introduced a framework representing the rendering equation in the form of a path integral for scenes without participating media (i.e., no scattering), and later, Pauly et al.¹⁷ extended it to the volume rendering equation with scattering. The amount of light $I$ observed by the detector is given by the path integral

Here, the bidirectional scattering distribution function $f_s(x_{m-1},x_m,x_{m+1})$ is used for locations on the surface of objects, and the scattering coefficient ${\sigma }_s(x_m)$ at $x_m$ and the phase function $f_p(x_{m-1},x_m,x_{m+1})$ are used for those inside the medium. $G(x_m,x_{m+1})$ is a generalized geometric term.

Putting all together, we have a path integral of the following infinite sum of all possible path contributions.

Note that all vertices $\{x_m\}$ depend on a path $k$; different paths have different sets of vertices. In the equation above, however, we omit the path index $k$ for simplicity. Later, we will again use $k$ as the path index.

3.2.

Assumptions on the Path Integral Formulation

As our target is optical tomography, we restrict the model to deal with inside the participating media. To do so, we assume that the light source $x_0$ and detector $x_{M+1}$ are located on the surface, and the other vertices $x_1,x_2,\dots ,x_M,x_{M+1}$ are inside the medium, that is, $x_0,x_{M+1}\in \partial \nu$ and $x_1,\dots ,x_M\in {\nu }_0$. Then the transmittance is simplified as

Furthermore, we assume that the observations are ideal and the camera response function is the identity, $W_e(x_M,x_{M+1})=1$.

Apart from the assumptions above, we rewrite the geometric term and the differential measure. The definitions above use area measures $dA(x_m)$ and volume measures $dV(x_m)$ along with the squared distance geometric term;^17,23,53 however, steradian measures $d\omega (x_m)$ and the identity geometric term is equivalent and also widely used.^10,12,60

Therefore, we employ the steradian measures and rewrite it as follows:

Now, Eq. (12) is written as

3.3.

Discretization of the Forward Model

For numerical computation, we first discretize the medium into voxels of a regular grid, where each voxel has its own extinction coefficient ${\sigma }_t[b]$ ($b$ is the index of the voxel) as shown in Fig. 1.

Fig. 1

Illustration of a discretization example. (a) Voxelization of the medium into a regular grid of size $5\times 5$. Voxels are indexed in raster scan order in this example, from left to right, and top to bottom. Each voxel $b$ has extinction coefficient ${\sigma }_t[b]$. (b) A path segment between vertices $x_1$ and $x_2$. Voxels involved in the segment are shaded. (c) Lengths $d_{12}[b]$ of the involved voxels $b=2,3,8,9$. Here we denote $d_{12}[b]$ instead of $d_{x_1,x_2}[b]$ for simplicity.

With this voxelization, the paths of light are also divided into segments, as explained below. First, we explain the integral [Eq. (11)] along a single segment $x_m,x_{m+1}$ of a path $\tilde{x}$. It describes the attenuation of light along the segment due to the extinction coefficients of the voxels involved. Because of the discretization of the medium, Eq. (11) can be written as a sum of voxel-wise multiplications.

For the second equality, $b$ is the index of a set ${\mathcal{B}}_{x_m,x_{m+1}}$ of all voxels involved by segment $x_m,x_{m+1}$, and $d_{x_m,x_{m+1}}[b]$ is the length of the part of the segment $x_m{x}_{m+1}$ passing through voxel $b$. This is illustrated in Fig. 1(c). The extinction coefficient ${\sigma }_t$ is now a piece-wise constant function because of the voxelization; then the integral turns into a sum (the idea that this integral can be turned into a sum has been discussed before,⁶¹ however, not in the context of tomography).

This simplifies the computation; however, the sum over a set ${\mathcal{B}}_{x_m,x_{m+1}}$ is not preferable in terms of implementation and optimization. We propose here to use a vector representation of both extinction coefficients and segment lengths, which is the third equality of the above equation. The first vector ${\mathit{\sigma }}_t$ stores the values of the extinction coefficients ${\sigma }_t[b]$ of all voxels. This vector can be generated by serializing the voxels on the grid in a certain order. The second vector ${\mathit{d}}_{x_m,x_{m+1}}$ contains the values of the lengths $d_{x_m,x_{m+1}}[b]$ for all voxels. We should note that this vector is very sparse; most of the voxels have no intersection with the segment $x_m,x_{m+1}$. Hence, only a few elements in ${\mathit{d}}_{x_m,x_{m+1}}$ have nonzero values, and the other elements are zero because those voxels $b$ have no intersection and $d_{x_m,x_{m+1}}[b]=0$.

This sparsity of the vector facilitates the construction of a whole path $\tilde{x}$ because path segments can be added as follows:

Using this notation to rewrite Eq. (17), we have

3.4.

Two-Dimensional Layered Model of Forward Scattering

As a first attempt, we design a two-dimensional (2-D) layered grid, instead of the three-dimensional (3-D) one. Since we voxelize the medium into a regular grid, the 2-D medium consists of parallel layers. Hereafter, a 3-D direction $\omega$ between vertices is written as a 2-D direction $\theta$ and a steradian measure $d\omega$ as an angular measure $d\theta$.

As shown in Fig. 2, we assume a particular layer scattering having the following properties. First, vertices $x_1,\cdots ,x_M$ of path $\tilde{x}$ are located at the centers of each voxel. The light source $x_0$ is located on the boundary of the top surface of the voxels in the top layer. Similarly, the detector $x_{M+1}$ is located on the boundary of the bottom surface of the voxels in the bottom layer. Second, directions ${\theta }_{x_0,x_1}$ and ${\theta }_{x_M,x_{M+1}}$ at the beginning and end of a path are perpendicular to the boundary. This means that scattering begins at $x_1$ and ends at $x_M$. Third, forward scattering happens layer by layer. More specifically, light is scattered at the center of a voxel in a layer and then goes to the center of a voxel in the next (below) layer. Scattering is assumed to happen every time the light traverses voxel centers. Even if the next voxel is just below the current voxel and the path segment is straight, it is regarded as scattering. Fourth, the scattering coefficient is uniform, ${\sigma }_s(x)={\sigma }_s$.

Fig. 2

Proposed two-dimensional layered model of scattering. This example shows path $\tilde{x}$ consisting of vertices $x_1,\cdots ,x_M$ located at the centers of voxels in a grid with $M$ parallel layers. $x_0$ is a light source located on the top surface, and $x_{M+1}$ is a detector at the bottom. At each vertex, the light scatters to voxels in the next layer, and possible scattering directions are indicated by arrows.

By ignoring paths exiting from the sides of the grid, the number of all possible paths is $N^M$, where $M$ is the number of layers and $N$ is the number of voxels in one layer.

3.5.

Approximating the Phase Function with a Gaussian

We use a Gaussian model $f_p(\theta ,{\sigma }^2)$ as an approximation of the phase function

First, existing phase function models^{10,12,35–38} are those for 3-D scattering, not for 2-D. This means that those functions are normalized for integrals over the unit sphere $S^2$: ${\int }_{S^2}{f}_p(\omega )\mathrm{d}\omega =1$. Most of the phase functions assume isotropy (rotational symmetry), and hence, the function has a form taking angle $\theta$ as an argument; however, ${\int }_{-\pi }^{\pi }{f}_p(\theta )\mathrm{d}\theta \ne 1$. These functions, therefore, are not adequate for our case.

Second, our assumption of layer-wise forward scattering does not allow scattering to happen backwards or sideways, and the Gaussian model is suitable for it. As shown in Fig. 3, the Gaussian model has the form of forward-only scattering (no backwards or sideways) in a reasonable range of ${\sigma }^2$, and it is almost normalized; ${\int }_{-\pi /2}^{\pi /2}{f}_p(\theta ,{\sigma }^2)\mathrm{d}\theta \approx 1$. Other 2-D phase functions exist which are not forward-only. For example, Heino et al.⁶² introduced a 2-D analog of Henyey-Greenstein’s phase function,³⁵ shown in Fig. 3. Although the parameters are different, the two functions in Fig. 3 have similar shapes. The most important difference is that Heino’s function has backward scattering, but our Gaussian model does not. More realistic scattering rather than the layer-wise forward scattering introduced here needs Heino’s or Henyey-Greenstein’s phase function.

Fig. 3

Comparison of two-dimensional phase functions. The upward vertical direction is $\theta =0$, and horizontal directions are $\theta =\pm (\pi /2)$. (a) Gaussian approximated phase functions with ${\sigma }^2=0.1,0.2,\dots ,1.0$. The tallest and narrowest shapes correspond to ${\sigma }^2=0.1$, and the shape becomes shorter and rounder for larger values of ${\sigma }^2$. (b) Heino’s two-dimensional analogs⁶² of Henyey-Greenstein’s phase function with parameter $g=0.1,0.2,\dots ,1.0$. The tallest and narrowest shapes correspond to $g=1.0$, and the shape becomes shorter and close to a hemisphere for smaller values of $g$.

We should note one further simplification in our layer-wise forward scattering model. The angle ${\theta }_m$ in the phase function is usually defined between ${\theta }_{x_{m-1},x_m}$ and ${\theta }_{x_m,x_{m+1}}$, that is, the difference of directions changed by the scattering event. Instead of dealing with such an exact difference of directions, we use the angle between ${\theta }_{x_m,x_{m+1}}$ and the vertical (downward) direction for efficiency of computation. This assumption enables us to discretize the Gaussian phase function much more easily. Since $f_p(\theta )$ integrates to (approximately) one, such a normalization can be discretized with a sum as follows:

Fig. 4

An illustration of angle measure $\mathrm{\Delta }{\theta }_b$ for voxel $b$ in the next layer. For the center voxel of the upper layer, voxel $b$ (shaded) in the next layer subtends an angle of $\mathrm{\Delta }{\theta }_b$, which is used for the angle measure in Eq. (24).

The above equation can be considered as the energy distribution from a voxel in one layer to the voxels in the next layer. For a voxel $b$ at direction ${\theta }_b$, the value of $f_p({\theta }_b,{\sigma }^2)\mathrm{\Delta }{\theta }_b$ describes what percentage of the energy will be scattered to this voxel. Figure 5 shows plots of the values corresponding to two phase functions with different parameters. We can see that, due to forward scattering, most of the energy is concentrated in the voxel just below, and a small part goes to the adjacent voxels.

Fig. 5

(a) The phase functions with parameter ${\sigma }^2=0.2$ (dashed line) and ${\sigma }^2=0.4$ (solid line). Plot of the value $f_p({\theta }_b,{\sigma }^2)\mathrm{\Delta }{\theta }_b$ for each voxel $b$ for (b) ${\sigma }^2=0.2$ and (c) ${\sigma }^2=0.4$. Note that index $b$ is relative to the voxel in the next layer just below the voxel in consideration. The voxel just below is $b=0$, the voxel on its right side is $b=1$, and that on the left side is $b=-1$.

The contribution $H_k$ in Eq. (22) now needs to be rewritten so that it deals with the Gaussian phase function and the discretized energy distribution discussed above. First, we reorder the measure

Note that the factor $dA(x_0)\mathrm{\Delta }{\theta }_{x_0,x_1}{\sigma }_s^M$ is common for all paths because we assumed that the grid is uniform so that $dA(x_0)$ is constant, and the direction ${\theta }_{x_0,x_1}$ (or ${\omega }_{x_0,x_1}$) is perpendicular to the top surface, and ${\sigma }_s$ is constant.

3.6.

Observation Model

Suppose the 2-D layered medium is an $M\times N$ grid; it has $M$ layers, each of which is made of $N$ voxels. We now construct an observation model of the light transport between a light source and a detector: emitting light to each of the voxels at the top layer, and capturing light from each voxel from the bottom layer. More specifically, let $i\in {\mathcal{B}}_1$ and $j\in {\mathcal{B}}_M$ be voxel indices of the light source and detector locations, respectively. By restricting the light paths to only those connecting $i$ and $j$, the observed light $I_{ij}$ is written as follows:

In the above equation, $k$ indexes the light paths, which share the same $i$ and $j$. Due to the layered scattering model in the $N\times M$ grid, the number of different paths between $i$ and $j$ is $N_{ij}=N^{M-2}$. This is, however, too large even for small $N$ and $M$, e.g., $N=M=10$. Therefore, we exclude paths having small contributions from the computation. This is done by a simple thresholding while computing $H_{ijk}$ as shown in Algorithm 1. This results in generating fewer paths; $N_{ij}\le N^{M-2}$. For example, there are $N_{ij}=742$ paths for $N=M=20$ with ${\sigma }^2=0.4$ when $\mathrm{th}=0.001$, which enables us to reduce the computation cost.

Algorithm 1

Computing contribution Hijk and omitting low contribution path by thresholding.

4.1.

Cost Function

In the $M\times N$ 2-D layered medium described in Sec. 3.6, we had assumed a configuration of a light source and detector similar to the left-most one shown in Fig. 6; the light source is located above the medium and the detector is below, and the observed light is $I_{ij}$, where $i,j$ are the voxel indices of the light source and detector locations. By sliding the light source and the detector, we can obtain $N^2$ observations, resulting in the following least-squares equation:

Fig. 6

Four configurations of light sources and detectors. From left to right, we call configurations T2B (top-to-bottom), L2R (left-to-right), B2T (bottom-to-top), and R2L (right-to-left), which represent locations of light sources and detectors.

Furthermore, as shown in Fig. 6, we have four configurations of light sources and detectors by changing their positions. This gives us four different sets of observations $I_{ij}$ and paths $ijk$. These four different sets lead to four objective functions ($f_{T2B}$, $f_{L2R}$, $f_{B2T}$, $f_{R2L}$) as shown in Fig. 6. Since the four objective functions share the same variables ${\mathit{\sigma }}_t$, we can use all of them at the same time by adding them to form a new single function $f_0$ at the expense of additional (factor of four) computation cost.

4.2.

Optimization Problem with Inequality Constraints

Since the inverse problem [Eq. (31)] is nonlinear, we employ an interior point method,²⁶ an iterative optimization algorithm for problems with constraints. Here, we first review several key points in optimization; then we will develop an algorithm to solve Eq. (31) along with the required first- and second-order derivatives of the cost function.

4.3.

Algorithm for Solving the Inverse Problem

Algorithm 3 shows our algorithm, which uses a barrier method with Quasi-Newton for solving the inverse problem. We should mention the following parts where we have modified the original algorithm.²⁶

Algorithm 3

Barrier method of interior point with Quasi-Newton solver.

Warm start: For each inner loop, the Quasi-Newton method needs an initial guess of the inverse Hessian $B_0$. Instead of fixing $B_0$ for every inner loop, we reuse the $B_k$ of the last inner loop to accelerate the convergence (shown in lines 4 and 19 in Algorithm 3).

Checking feasibility: Since the Quasi-Newton method and line search estimate without constraints, the next estimate ${\mathit{\sigma }}_t^{k+1}$ may go beyond the constraints; in our case, each element ${\sigma }_t^{k+1}[b]$ in ${\mathit{\sigma }}_t^{k+1}$ must be inside $[0,u]$ after the step size has been determined. Therefore, in line 8, we check the feasibility of the estimate ${\mathit{\sigma }}_t^{k+1}$ for the current step size ${\alpha }_k$. If it exceeds the boundary of the feasible region, we pull the estimate back into the feasible region by halving the step size. If it is still outside the feasible region, then the step size is halved again. Why do we not just set the step size so that ${\mathit{\sigma }}_t^{k+1}$ is exactly on the boundary? The reason is the log-barrier: if ${\mathit{\sigma }}_t^{k+1}$ is on the boundary, in other words, ${\sigma }_t^{k+1}[b]$ is either 0 or $u$, then $\log ({\sigma }_t[b])$ or $\log (u-{\sigma }_t[b])$ becomes infinite, which results in numerical instability. Therefore, the procedure described above is needed.

Checking for positive definiteness: The BFGS update rules guarantee $B_k$ to be positive definite if ${\mathit{y}}^T\mathit{s}>0$ and $B\succ 0$ are satisfied. While the latter is satisfied by giving an appropriate initial guess, the former depends on the updates at each iteration. If it is not satisfied, then the BFGS updates are no longer valid, and we reset the inverse Hessian $B_k$ to a scaled identity⁶³ at line 16.