1. Introduction

In diffuse optical tomography, one seeks to reconstructs the optical properties of tissues using measurements of scattered near-infrared light. This imaging modality has applications in neonatal brain imaging and breast imaging, for example (see [1, 2] and references contained therein). The major challenge in this field is taking into account the strong multiple scattering of light in tissues.

Light propagation in tissues is governed by the theory of radiative transport [3]. This theory takes into account absorption and scattering due to inhomogeneities. The radiative transport equation is an integral-differential equation for the specific intensity. This equation is valid when distances between scatterers are much larger than the wavelength. Consequently, this equation does not take into account interference effects between scatterers. There are very few specific cases for which we know the analytical solution to the radiative transport equation. In addition, it is difficult to solve the radiative transport equation numerically, especially when the medium is optically thick or scatters with a sharp peak in the forward direction. Both of these situations are relevant to light propagation in tissues. Forward-peaked scattering requires restrictively large computational resources to approximate the integral operator in the transport equation accurately.

There are several approximations to the radiative transport equation that take into account forward-peaked scattering [4]. Comparisons of the reflectance and transmittance of light by tissues computed using these approximations show good agreement with those computed by the radiative transport equation. These approximations require less computational resources to solve than the radiative transport equation. Hence, these approximations will be useful for diffuse optical tomography if they produce good quality image reconstructions.

To investigate the quality of image reconstructions that these approximations produce, we solve inverse problems using these approximations as models for the measured data and examine the quality of the images they produce. In particular, we compute image reconstructions produced by five light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation. Our simulated data is computed using the radiative transport equation with forward-peaked scattering governed by the Henyey-Greenstein scattering phase function.

For the solution of the inverse problem, we use the transport-backtransport method [5, 6]. The transport-backtransport method is one example of the more general propagation-backpropagation method. The propagation back-propagation method has been used widely in medical imaging [1, 7, 8], geophysical applications [9, 10, 11, 12] and other general applications as shown in [13, 14]. This inversion method is comprised of two parts: (1) solving the direct problem and (2) calculating updates for the parameters that we seek to reconstruct. We use the five light scattering models mentioned above for part (1) of the inversion method to reconstruct five different sets of images. Through comparison of these five sets of reconstructed images, we determine which of the approximations are suitable for diffuse optical tomography.

Because we are using the transport-backtransport method to solve the inverse problem, we solve only one direct and one adjoint problem per iteration regardless of the number of unknowns in the problem. In contrast, direct gradient calculation methods require the solution of one direct problem for each parameter in the problem. For the problems we study here, we are dealing with approximately 2500 parameters. Hence, direct gradient calculation methods would require restrictively large computations.

The inverse problem involves the reconstruction of the spatial distribution of the absorption and scattering coefficients from scattered light measurements at the boundary. This problem is severely ill-posed due to multiple scattering. As a result, we do not anticipate a perfect, detailed reconstruction. In fact, one artifact that we will see in our reconstructions is the so-called “crosstalk” in which images of the absorption coefficient show artifacts from the scattering coefficient, and vice versa. As a result, one method to compare the performance of the approximations in image reconstructions is to evaluate the amount of cross-talk in the images.

The remainder of this paper is as follows. In Section 2 we discuss the direct and inverse problems for diffuse optical tomography. Included in that discussion is a description of the radiative transport equation, which governs light propagation in tissues and the transport-backtransport method which we use to solve the inverse problem. In Section 3 We discuss the different approximations of the radiative transport equation for forward-peaked scattering. We give the numerical method we use to solve the direct problem in Section 4. In Section 5 we show the results of the numerical simulations. Section 6 gives our conclusions. Appendix A gives a proof of a result used in the transport-backtransport method in Section 2. Appendix B discusses the self-adjoint properties of the light scattering models we consider here.

2. Diffuse optical tomography

For diffuse optical tomography, we study two mathematical problems: the direct problem and the inverse problem. For the direct problem, we seek to model measurements of scattered light for a given set of optical properties. For the inverse problem, we compare our model for the measurements with the actual measurements to reconstruct the optical properties. In what follows, we discuss both the direct and inverse problems.

2.1. The direct problem

To model light propagation in tissues, we solve the radiative transport equation:

Here, S ^n-1 denotes the set of unit direction vectors, Ω ⊂ Rⁿ denotes the spatial domain with boundary ∂Ω and [0,T] is the interval of time we consider. The absorption coefficient is denoted by μ_a, the scattering coefficient is denoted by μ_s and the source is denoted by Q. The scattering operator 𝓛 is defined as

The scattering phase function f gives the fraction of light scattered in direction ŝ due to light incident in direction ŝ′. We discuss specific choices for 𝓛 in Section 3.

In addition to Eq. (1), we must prescribe boundary and initial conditions. Boundary conditions give the light for all directions pointing into the spatial domain for all points at the boundary over all time. For our boundary conditions, we prescribe that

with ν(r _b) denoting the unit outward normal at r _b ∈ ∂D. The initial condition gives the intensity at t = 0 for all directions and all points in the spatial domain. For the initial condition, we prescribe that

Hence, the only source of light in the direct problem is given by the source Q. This source Q may correspond to a source at the boundary of the spatial domain.

2.2. The inverse problem

For measurements, we consider angle-averaged quantities at the boundary that take the form

Here, S _out ^n-1 = {ŝ ∈ S ^n-1,ŝ · ν(r _d) < 0} with ν(r _d) denoting the unit outward normal at r _d ∈ ∂Ω, is the set of all angles pointing out of the spatial domain. Suppose we have collected a set of measurements of the form given by Eq. (5) for each source within a set of sources. For the inverse problem, we seek to reconstruct μ_a and μ_s using this data and the model for measurements Eq. (5) with I denoting the solution of the direct problem given by Eq. (1) subject to Eq. (3) and Eq. (4).

This inverse problem is nonlinear since the solution of the direct problem depends nonlinearly on the parameters μ_a and μ_s. Moreover, this inverse problem is severely ill-posed due to multiple scattering. Finally, we will study inverse problems that are underdetermined since it is typical in practical situations to have only a few detectors placed around the spatial domain.

To address these challenges, we seek to solve the inverse problem by solving an associated optimization problem. We denote the space of μ_a or μ_s by P and the space of data M _μa,μs by D. The residual is given by the difference between the real data F generated by experimental measurements and the modeled data (i.e. the whole set of sources and receivers) corresponding to the parameter distribution (μ_a,μ_s):

We introduce the cost functional:

Now, we turn our attention toward determining the parameter distributions for μ_a and μ_s that minimize K. That minimum norm solution yields the image reconstructions for μ_a and μ_s that we are seeking.

2.3. The transport-backtransport method

We use the transport-backtransport method developed by O. Dorn [5, 6] to determine the parameter distributions for μ_a and μ_s that minimize K given by Eq. (7). We do not describe all of the details of this method here. Rather, we give only the main results.

By perturbing the parameters about μ̅_a and μ̅_s as μ_a = μ̅_a+δμ_a and μ_s = μ̅_s+δμ_s, we find that the linearized residual operator about (μ̅_a,μ̅_s) is given by

with R′ denoting the Fréchet derivative of R [15]. By setting the right-hand side of Eq. (8) equal to zero, we obtain a linear system for the perturbations (δμ_a, δμ_s). The minimum norm solution of this linear system is given by

Here, R̃′(μ̅_a,μ̅_s) denotes the adjoint of R′(μ̅_a, μ̅_s). We prove in Appendix A that this minimum norm solution is a descent direction of the functional K(μ_a, μ_s).

Included in the minimum norm solution Eq. (9) is the operator [R′(μ̅_a,μ̅_s)R̃′(μ̅_a,μ̅_s)]^-1. This operator is time consuming to calculate in practice because it requires solving as many direct and adjoint problems as there are data points. Instead, we approximate this operator by the scalar factor β multiplied to the identity operator, so that

This approximation is substantially faster to compute. Moreover, we show in Appendix A that this approximate solution is a descent direction of the functional K(μ_a, μ_s) as well. To compute the update Eq. (10), we compute first r = R(μ̅_a,μ̅_s). This residual defined by Eq. (6) requires the solution of the direct problem with parameters μ̅_a and μ̅_s. Next, we must compute R̃′(μ̅_a,μ̅_s)r. This operation is given in terms of Ĩ, the solution of the adjoint problem [5, 6] given by

This form for the adjoint problem relies on 𝓛 being self-adjoint. We consider this matter in detail in Appendix B.

Solving Eq. (11) is called backtransport. Hence, the transport-backtransport method is comprised of computing the direct problem followed by the adjoint problem. Upon solution of the direct and adjoint problems, we compute the components of Eq. (10) and thus obtain updated values given by

2.4. Procedure to solve the inverse problem

We now describe the procedure we take to solve the inverse problem. Suppose we have chosen an initial distribution for the absorption and scattering coefficients. Typically, we choose constant values for the initial distributions of μ_a and μ_s. Next, we choose an initial step size λ. Included within λ is the constant β given in Eq. (10). In the procedure below, this parameter λ changes to accelerate the convergence.

Once we have set these initial values, we proceed to solve the inverse problem in the following way.

It is possible that the line search method described above may affect the overall performance of a particular light scattering model in a way that is different for another light scattering model. Hence, using this algorithm in making absolute comparisons between the different light scattering models may be subject to objection. For that reason, we have taken care in analyzing the results that follow to avoid over-statement of the results.

3. Models for scattering operators

For the radiative transport equation Eq. (1), we must specify the scattering operator 𝓛. Tissues typically scatter light with a sharp forward-peak. For that case, the scattering phase function f given in Eq. (2) is is peaked about ŝ · ŝ ~ 1. Forward-peaked scattering makes the radiative transport equation difficult to solve. For that reason, we discuss several approximations that model forward-peaked scattering that we will consider in this study. These approximations are derived in detail in Ref. [4]. Here, we just present the approximations.

3.1. Henyey-Greenstein scattering

The Henyey-Greenstein scattering phase function,

is used typically for light propagation in tissues. Here, g denotes the asymmetry parameter or anisotropy factor and P_n denotes the nth order Legendre polynomial. For forward-peaked scattering, the anisotropy factor is g ~ 1. Using Eq. (13) in Eq. (2), we obtain 𝓛_HG which we call the Henyey-Greenstein scattering operator.

3.2. The delta-Eddington approximation

The delta-Eddington approximation was introduced by Joseph et al. [16] as an extension to the two-term Eddington approximation to study forward-peaked and large-angle scattering. It approximates the scattering phase function as the sum of a delta function and a smooth background. When applying the delta-Eddington approximation to the Henyey-Greenstein scattering phase function, we obtain the following delta-Eddington scattering operator:

3.3. The Fokker-Planck approximation

The Fokker-Planck equation was first proposed by A. D. Fokker [17] and M. Planck [18] to describe the Brownian motion of a particle. In the context of light propagation in tissues, the Fokker-Planck approximation is given by [19]

with ∆_ŝ denoting the spherical laplacian.

3.4. The Fokker-Planck-Eddington approximation

Recently, the authors have combined the delta-Eddington approximation and the Fokker-Planck approximation to obtain the Fokker-Planck-Eddington approximation [4]:

For the Henyey-Greenstein scattering phase function defined in Eq. (13), the coefficients in Eq. (16) are given by

3.5. The generalized Fokker-Planck-Eddington approximation

In addition to the Fokker-Planck-Eddington approximation, the authors have introduced the generalized Fokker-Planck-Eddington approximation [4]. This approximation uses a higher-order Fokker-Planck approximation along with the two-term Eddington approximation. This approximation leads to the generalized Fokker-Planck-Eddington scattering operator 𝓛_gFPE defined as

with ℐ denoting the identity operator. For the Henyey-Greenstein scattering phase function defined in Eq. (13), the coefficients in Eq. (18) are given by

3.6. Two dimensional scattering

The approximations explained above are for the three dimensional scattering. For our computations, we consider only two spatial dimensions x and y and scattering is restricted to the xy-plane. For this case, the scattering phase function is slightly different and so are the corresponding approximations.

According to J. Heino et al [20], the Henyey-Greenstein phase function in two dimensions is:

The corresponding Henyey-Greenstein scattering operator is given by

Consequently, the approximations for the scattering operator are given by the following.

4. Numerical solution of the direct scattering problem

We describe the numerical solution of

in [0,2π]×[0,L_x]×[0,L_y]×[0,T]. Here, we describe the numerical solution for a general scattering operator 𝓛. Later, we give numerical approximations that we use for computing numerical approximations of 𝓛_HG defined by Eq. (21), 𝓛_δE defined by Eq. (22), 𝓛_FP defined by Eq. (23), 𝓛_FPE defined by Eq. (24) and 𝓛_gFPE defined by Eq. (26).

4.1. Source iteration method

To solve Eq. (28), we use the source iteration method [21]. For the source iteration method, we solve numerically

for l = 1,2,⋯ until we reach convergence:

with ε denoting a user-defined termination parameter. Here, the solution we seek I is given by I^l as l → ∞. The source J^l is defined as

For each iteration, the numerical method is comprised of two parts: computing the source update Eq. (31) and solving the partial differential equation Eq. (29). In what follows, we describe each of these two parts in detail.

4.2. Computing the source update

Depending on the scattering operator we are considering, the source update involves either an integral operator, a differential operator or a combination of both. For integral operators, we use a quadrature rule and for the differential operators, we use a finite difference scheme.

For both types of operators, we consider the angular grid of N points defined as θ_n = (n - 1)∆θ with ∆θ = 2π/N for n = 1, ⋯,N. We now introduce the approximant I_n(x,y,t) ≈ I(θ_n,x,y,t) that we will use in the numerical approximations below.

For integral operations, we use a product quadrature rule [22] of the form:

Here, the weights f_nk are defined as

For the derivative operations, we use a centered finite difference scheme of the form

with periodic numerical boundary conditions: I _-1 = I_N and I _N+1 = I ₁.

By using either Eq. (31), Eq. (33) or a combination of the two, we form the discrete scattering operator 𝓛^N. Using 𝓛^N, we compute the source update at θ_n through evaluation of

4.3. Diamond differencing scheme

Using the angle grid defined above, we seek to solve numerically the partial differential equation

with J_n known. To solve this problem, we use the diamond differencing scheme [21]. An equispaced spatial grid defines spatial cells each of size ∆x by ∆y. We use integer indices i and j to denote each cell. Furthermore, we assume that μ_a and μ_s are constant within each cell. Then, we seek the cell-centered average of I_n(x,y,t) over each cell defined as

Here, we have denoted the edges of spatial cells by the fractional indices i ± 1/2 and j ± 1/2.

To derive an equation for I_nij(t), we integrate Eq. (36) over [x _i-1/2,x _i+1/2] × [y _j-1/2,y _j+1/2], divide by ∆x∆y and obtain the semi-discrete system

Here, we have introduced the cell-edge quantities:

Next, we seek a method to evolve the solution from time t_k to t _k+1 = t_k+∆t. Integrating Eq. (38) with respect to t over [t _k-1/2,t _k+1/2] and dividing by ∆t, we obtain

Here, we have introduced the notation I_nij ^k+1/2 ≈ I_nij(t _k+1/2), I_nij ^k-1/2 ≈ I_nij(t _k-1/2) and

The diamond differencing scheme uses the following approximations to go back-and-forth from cell-averaged quantities (fractional indices) to cell-edge quantities (integer indices):

By combining Eq. (40) and Eq. (42), we obtain a discrete system of equations whose solution yields the numerical solution we seek. For the diamond differencing scheme in time, we substitute 2I_nij^k - I_nij ^k+1/2 for I_nij ^k-1/2 in Eq. (40) and obtain an update formula for I_nij ^k+1/2. Once we determine that update, we compute the full update through evaluation of I_nij ^k+1 = 2I_nij ^k+1/2 -I_nij ^k. The actual form of the update formula for I_nij ^k+1/2 depends on the signs of cos θ_n and sin θ_n. In particular, we use the following.

So for example, we solve the following system of equations to compute the I_nij ^k+1/2 update for cos θ_n ≥ 0 and sin θ_n ≥ 0:

To ensure stable time-stepping, we must set $c\mathrm{\Delta }t\le \genfrac{}{}{0.1ex}{}{1}{2}min\left\{∆x,∆y\right\}.$. In addition, we must set c∆tmax_ij{μ_a,ij} < 1 [5, 21]. To avoid this second restriction, we have implemented a semi-implicit time-step. Thus, instead of solving Eq. (43), we solve

5. Numerical results

We investigate the quality of image reconstructions produced by the five light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation. For all numerical experiments, we obtain simulated boundary measurements by solving the radiative transport equation with the Henyey-Greenstein scattering phase function with g = 0.9. Using those simulated data, we use the transport-backtransport method to reconstruct images of the absorption and scattering coefficients. We perform the transport-backtransport method using each of the five light scattering models to obtain five sets of image reconstructions.

For our numerical calculations, the spatial domain is a 2 cm × 2 cm square with constant background optical properties: μ_a = 0.1mm^-1 and μ_s = 1.0 mm^-1. Within this medium, there is an absorbing obstacle whose absorption coefficient is 1.0 mm^-1 and a scattering obstacle whose scattering coefficient is 5.0mm^-1. These parameter values are reasonable for studies of light propagation in tissues. In particular, we have chosen these parameters to focus on the regime in which the scattering mean free path l_s = μ_s ^-1 is significantly smaller than the transport mean free path I ^* = [μ_s(1 - g)]^-1. For that case, the medium scatters light strongly, but with a strong forward peak. This regime is where the approximations are designed to work well. Figures 1(a) and 2(a) show the placement and size of the absorbing and scattering obstacles, respectively.

For all of our computations of the direct and adjoint problems, we have set the number of discrete angles to be N = 64 so that ∆θ = π/32. The spatial grid has 51 points in x and y so that ∆x = ∆y = 0.04cm. We have taken a time step of ∆t = 0.1 sampling an 8 second interval of time (we have set c = 1 here for convenience).

Along the boundary of the spatial domain, we have placed forty, equi-distant sources and detectors surrounding the entire spatial domain. In fact, these sources and detectors are collocated spatially. Thus, for a particular source, the set of detectors correspond to all source/detector locations except the ones located in the same side of the domain.

The results of these numerical simulations are shown in Figs. 1 and 2. In all of the reconstructed images, we observe the presence of “cross-talk” in which the scattering obstacle appears as an artifact in the absorption coefficient image and the absorbing obstacle appears as an artifact in the scattering coefficient image. Even the so-called inverse crime solution in which the radiative transport equation is used as the light scattering model exhibits cross-talk. This cross-talk is due to the ill-posedness of the inverse problem. In particular, the data we have are limited to just a few spatial locations (albeit surrounding the entire domain) of angle-averaged quantities. Using a priori information or incorporating more data into the inverse problem may alleviate this problem.

Even though all of the light scattering models produce images with substantial cross-talk, all of the reconstructed images are reasonably good. For example, the absorbing and scattering obstacles are localized correctly. We observe that the reconstructions of the scattering coefficient are better overall than the respective images for the absorption coefficient. Differences occurring between the images reconstructed using the different light scattering models include better localization and less image artifacts. For example, the reconstructions for delta-Eddington, Fokker-Planck and Fokker-Planck-Eddington models show image artifacts near the source positions. Additional regularization near the sources may remove these artifacts provided we the reasonable assumption that the obstacles to be reconstructed are not close to the sources.

For the absorption coefficient image reconstructions shown in Fig. 1, the cross-talk is more severe than for the scattering coefficient. Not surprisingly, the inverse crime solution in which the image is produced by the radiative transport equation is the best one. The worst qualitative and quantitative image is produced by the Fokker-Planck approximation. The delta-Eddington approximation produces a slightly better image and the Fokker-Planck-Eddington approximation produces an even better image. Besides that inverse crime solution, the generalized Fokker-Planck-Eddington scattering model produced the best image.

For the scattering coefficient, we find again that the inverse crime solution produces the best image. Besides that inverse crime solution, we find that the generalized-Fokker-Planck-Eddington approximation produces the best image. The Fokker-Planck approximation produces a slightly better qualitative result than the delta-Eddington approximation since there are less image artifacts. However, the delta-Eddington approximation yields a better quantitative result. The Fokker-Planck-Eddington approximation is better than delta-Eddington and Fokker-Planck, but worse than the generalized Fokker-Planck-Eddington approximation.

In Fig. 5 we show the evolution of the cost functional Eq. (7) as a function of the successful iterations using the different light scattering models. This plot gives another perspective on the performance of the different light scattering models. For example, we see clearly that the radiative transport equation performs best by a substantial amount. The wide variety of cost functional values are due largely to modeling error. When the light scattering model agrees better with the radiative transport equation, we find that the cost values are smaller than those from light scattering models that do not agree as well. That is the reason why we find that the generalized Fokker-Planck-Eddington cost values are smaller than the delta-Eddington or Fokker-Planck ones. In fact, we find that the generalized Fokker-Planck-Eddington approximation performs best among the approximations. In these computations, we have set a limit on the number of consecutive iterations in which the residual does not decrease rather than setting a tolerance value (Step 3 (b) in the procedure given in Section 2.4). It is for that reason that we see that the delta-Eddington approximation gives the poorest performance in this context – it reaches the maximum number of rejected iterations faster than the others without reducing the cost functional substantially.

 

Fig. 1. The actual absorption coefficient (a) and reconstructions of the absorption coefficient produced by the five light scattering models (b)–(f).

Download Full Size | PDF

 

Fig. 2. The actual scattering coefficient (a) and reconstructions of the scattering coefficient produced by the five light scattering models (b)–(f).

Download Full Size | PDF

 

Fig. 3. Comparison of the values of the cost functional as a function of successful iterations for the radiative transport equation (RTE), the generalized Fokker-Planck-Eddington approximation (gFPE), the Fokker-Planck-Eddington approximation (FPE), the Fokker-Planck approximation (FP) and delta-Eddington approximation (DE).

Download Full Size | PDF

Table 1. Time required to compute the numerical solution of a single direct problem in seconds.

View Table

An important consideration in using one of the various approximations for image reconstruction is the computational savings they afford over using the radiative transport equation. For the computational results shown here, we have used a workstation with two dual core Intel Xeon (3.0 GHz) processors. In Table 1 we give the time required to solve a single direct problem for each of the light scattering models. Using any of the approximations instead of the radiative transport equation is at least 25% faster. This computational savings is compounded even further since the transport-backtransport method requires the computation of several direct and adjoint problems. The implementation of the numerical solution may be implemented on parallel cluster thereby making these calculations faster overall. For example, one may run each solution corresponding to an individual source on different processors. Doing so would nearly linearly with the number of processors. However, the relative speeds among the different light scattering models are still an important factor in the overall computational time.

6. Conclusions

We have investigated the use of approximations of the radiative transport equation with forward-peaked scattering for producing image reconstructions. Using the transport-backtransport method for solving the inverse problem, we have used five different light scattering models: (1) the radiative transport equation, (2) the delta-Eddington approximation, (3) the Fokker-Planck approximation, (4) the Fokker-Planck-Eddington approximation and (5) the generalized Fokker-Planck-Eddington approximation to produce five sets of image reconstructions. By comparing those results, we have evaluated the potential use of approximations for reconstructing images in diffuse optical tomography. These approximations provide substantial computational savings compared with the computational resources required by the radiative transport equation.

All of our results exhibit cross-talk in which the absorbing obstacle appears in the scattering coefficient image and the scattering obstacle appears in the absorption coefficient image. Crosstalk occurs even in the inverse crime scenario in which the radiative transport equation computes the data and is the light scattering model for the inverse problem solution. It is an artifact of the ill-posedness inherent in the inverse problem. Hence, it has become clear through this investigation that eliminating or reducing cross-talk, regardless of the light scattering model used in the inverse solution, is an important topic. We will address cross-talk in our future work.

Overall, we have found that all of the light scattering models produce reasonable results. Besides the inverse crime solution that uses the radiative transport equation, we have found that the generalized Fokker-Planck-Eddington approximation produces the best qualitative and quantitative images. As expected, the quality of image reconstructions is directly proportional to the computational time required to solve a single direct problem. Nonetheless, since the approximations are at least 25% faster to solve than the radiative transport equation, they are certainly worthwhile to consider for practical applications.