## Abstract

We present a computational study of diffuse optical tomography using the one-way radiative transfer equation. The one-way radiative transfer is a simplification of the radiative transfer equation to approximate the transmission of light through tissues. The major simplification of this approximation is that the intensity satisfies an initial value problem rather than a boundary value problem. Consequently, the inverse problem to reconstruct the absorption and scattering coefficients from transmission measurements of scattered light is simplified. Using the initial value problem for the one-way radiative transfer equation to compute the forward model, we are able to quantitatively reconstruct the absorption and scattering coefficients efficiently and effectively for simple problems and obtain reasonable results for complicated problems.

© 2015 Optical Society of America

## 1. Introduction

In diffuse optical tomography (DOT), one seeks to reconstruct tissue optical properties, such as the absorption and scattering coefficients, using non-invasive scattered light measurements. DOT requires the solution of two mathematical problems: the *forward model* and the *inverse problem*. The forward model predicts the scattered light measurements for a particular set of tissue optical properties. Effective forward models take into account specific details about the light sources, the complex interactions of light with tissues, and the characteristics of the detectors used to measure the scattered light. The inverse problem uses the forward model to reconstruct tissue optical properties from scattered light measurements. It is mathematically ill-posed due to several factors including insufficient data diversity and inherent instabilities to noise. Typical methods to solve the inverse problem find the optical properties which, when substituted into the forward model, yield the best fit (in a particular sense) of the predicted measurements to the actual measurements.

Light propagation in tissues is governed by the radiative transfer equation (RTE) [1, 2]. The RTE accurately takes into account absorption and multiple scattering of light in tissues. Using the RTE as a forward model for DOT involves prescribing a boundary value problem for it that incorporates the specific characteristics of the imaging setup. Posing this boundary value problem as a forward model is physically intuitive and poses no difficulty. However, computing solutions for this boundary value problem is challenging. Even computing solutions numerically requires sophisticated methods and substantial computing resources. While an individual computation of the forward model may not be too restrictive, solving the inverse problem usually requires computing the forward model several times for a single reconstruction. In fact, the forward model is typically the bottleneck in algorithms to solve the inverse problem. It is for this reason that forward models based on solving boundary value problems for the RTE are typically considered too computationally intensive to use for practical DOT problems.

For many DOT problems, the diffusion approximation to the RTE is used to compute the forward model. The diffusion approximation is a simplification of the RTE for optically thick media in which multiple scattering is dominant [1, 2]. Upon using the diffusion approximation for the forward model, one obtains a boundary value problem for the diffusion equation, a partial differential equation which is much simpler than the RTE. Even though the diffusion approximation suffers from large errors near sources and boundaries, it has been used in forward models to solve several important DOT problems (see the reviews [3, 4]).

An important problem in biomedical optics requires imaging tissue systems that are in the so-called mesoscopic scattering regime [5]. Mesoscopic scattering refers to tissue system sizes that are large enough that light undergoes multiple scattering, but not so large that this scattered light is diffusive. For mesoscopic scattering, the diffusion approximation is not sufficiently accurate to make its use in a forward model useful. Instead, one must either methods based on the RTE [6–9] or forward models that combine radiative transfer with diffusion [10,11]. A new forward model that is physically intuitive, accurate, and easy to solve would be useful to solve DOT problems in the mesoscopic scattering regime.

Here, we introduce a new forward model for DOT that uses the one-way RTE. The one-way RTE is a simplification of the RTE for transmission of a source through an anisotropic, forward-peaked scattering medium. It is derived simply by restricting limits of integration in the scattering operator of the RTE to the half-range of directions aligned with the source incident on the boundary. Since the derivation of this approximation is based only forward-peaked scattering, this approximation is valid across a broad range of absorption and scattering coefficients. Rather than having to solve a boundary value problem, this forward model requires the solution of a simpler initial value problem. This reduction from a boundary value problem to an initial value problem makes solving the inverse problem simpler and more efficient.

The purpose of this work is to show that using the initial value problem for the one-way RTE to compute the forward model for a DOT problem is effective and efficient. Using this forward model, we reconstruct quantitatively accurate images of the absorption and scattering coefficients from simulated measurements computed from the solution of the boundary value problem for the RTE. Our results show that this forward model is quite effective.

The remainder of this paper is as follows. In Section 2 we describe the imaging system that we will consider for this study. This imaging system is idealized to some extent, so we list our assumptions explicitly here. In Section 3 we describe the boundary value problem for the RTE for this imaging system. In Section 4, we derive the initial value problem for the one-way RTE that is an approximation of the forward model presented in Section 2. In Section 5 we give a discrete ordinate/finite difference method to compute the solution of the initial value problem for the one-way RTE. In Section 6 we describe the inverse problem and the method we use to solve it to obtain our reconstructions of the absorption and scattering coefficients. Section 7 gives the results of our simulations along with a discussion of those results. We give our conclusions in Section 8.

## 2. Imaging system

Consider the two-dimensional imaging system depicted in Fig. 1. An intensity-modulated plane wave source is incident normally on the source plane, *y* = 0. At the detector plane, *y* = *y*_{0}, we measure the transmittance by a slab containing a disk, each of which is composed of a uniform turbid medium except that the disk contains two obstacles. One of those obstacles (labeled as the “scatterer” in Fig. 1) scatters more than the background turbid medium, and the other (labeled as the “absorber” in Fig. 1) absorbs more than the background turbid medium. We obtain multiple transmittance measurements of the disk by rotating it. We focus our attention here on this two dimensional imaging system to keep the computations manageable and presentation of results simple.

For this imaging system, we make the following assumptions.

- The radius of the disk and its location in the slab are known.
- The optical properties of the slab exterior to the disk and the background optical properties within the disk are known.
- The refractive index is the same outside and inside the slab, including the disk, the absorber, and the scatterer.
- The anisotropy factor is fixed throughout the slab and is known.

With these assumptions, we seek to reconstruct quantitatively accurate images of the scatterer and absorber from the multiple transmittance measurements.

## 3. The radiative transfer equation

Light propagation in tissues is governed by the theory of radiative transfer [1,2]. The RTE takes into account absorption and scattering due to inhomogeneities. Consider the two dimensional slab, {(*x,y*) : *−*∞ < *x* < ∞,0 < *y* < *y*_{0}}. For the imaging system described above, we solve the following boundary value problem for the RTE in this slab

Here, *ω* is the intensity modulation frequency, *c* is the speed of light, and *μ _{t}* =

*μ*+

_{a}*μ*where

_{s}*μ*is the absorption coefficient and

_{a}*μ*is the scattering coefficient. The scattering phase function,

_{s}*p*, gives the fraction of light scattered in direction

**ŝ**= (cos

*θ*, sin

*θ*) due to light incident in direction ŝ′ = (cos

*θ*′, sin

*θ*′). For this study, we use the Henyey-Greenstein scattering phase function,

Here, *g* is the anisotropy factor. The scattering phase function is normalized according to

Typically *g ∼* 1 in tissues, so scattering is forward-peaked. Boundary condition (3.1b) prescribes a plane wave of intensity *I*_{0} incident on the boundary *y* = 0 in direction **ŝ** = ŷ. Boundary condition (3.1c) prescribes that no light enters into the domain from the boundary *y* = *y*_{0}.

Upon solving boundary value problem (3.1), we compute the transmittance, *T*(*x*), through evaluation of

*NA*is the set of angles corresponding to the numerical aperture of the detector. Because of the intensity modulation,

*I*is complex and so

*T*is complex. We will make use of both the real and imaginary parts of

*T*in the image reconstruction algorithm we describe below.

To compute simulated transmittance data for this study, we solve boundary value problem (3.1) numerically using a discrete ordinate method along with finite differences which we describe below in Section 5.4. Because this problem is restricted to two spatial dimensions, this solution can be computed using reasonable computing resources. However, solving the associated inverse problem to reconstruct quantitative images of *μ _{a}* and

*μ*requires the computation of several solutions of boundary value problem (3.1). The computations required to solve boundary value problem (3.1) in that case become prohibitively costly. This issue becomes even more pronounced for three dimensional computations. Rather than solving boundary value problem (3.1) for the the image reconstruction algorithm, we describe below the one-way radiative transfer equation to approximate the transmittance.

_{s}## 4. The one-way radiative transfer equation

The one-way RTE is an approximation to the RTE that leads to a simpler initial value problem [12]. To derive the one-way radiative transfer equation for this problem, we introduce the following half-range intensities,

where*I*corresponds to the

^{±}*±*ŷ directions, respectively. Substituting Eq. (4.1) into boundary value problem (3.1), we obtain

Here, we have introduced the operators ${\mathcal{P}}_{f}$, defined as

We call ${\mathcal{P}}_{f}$ the forward scattering operator because it gives the light scattered in the half-range of directions aligned with the incident direction. We call ${\mathcal{P}}_{b}$ the backward scattering operator because it gives the light scattered in the half-range of directions opposite to the incident direction. Boundary value problem (4.2) is just a reformulation of boundary value problem (3.1), so it is equivalent. In fact, we introduce in Section 5, a numerical method to solve boundary value problem (4.2) rather than one for boundary value problem (3.1). Once we have solved boundary value problem (4.2), we compute the transmittance through evaluation of

When scattering is forward peaked, the scattering phase function, *p*, is peaked about *θ* = *θ*′. Let *M _{f}* (

*θ*) denote the “mass” contained in

*θ −π/*2 ≤

*θ*′ ≤

*θ*+

*π*/2,

For isotropic scattering with *g* = 0, *M _{f}* = 1/2 since the mass is split evenly between the forward and backward directions. For forward peaked scattering,

*M*1 meaning that nearly all of the mass is contained in the forward scattering directions. In Fig. 2 we plot

_{f}*M*(

_{f}*θ*) for the Henyey-Greenstein scattering phase function given in (3.2) for

*θ*=

*π/*2.

In Fig. 2 we find that for *g ≥* 0.73 that over 90% of the mass is contained in the forward directions. For that case, we have

In light of Eq. (4.7), we neglect the term involving ${\mathcal{P}}_{b}{I}^{-}$ in Eq. (4.2a) and obtain

Equation (4.8) is decoupled from *I ^{−}*. Therefore,

*y*is a time-like variable that ranges from

*y*= 0 to

*y*=

*y*

_{0}. Hence, we need to supplement only Eq. (4.8) with Eq. (4.2c), which now serves an initial condition for Eq. (4.8). We call Eq. (4.7) with Eq. (4.2c) the initial value problem for the one-way RTE.

The one-way RTE will be useful for the image reconstruction algorithm we give in Section 6. In preparation for solving the inverse problem, we introduce here the adjoint problem for the one-way RTE:

This adjoint problem is a final value problem since Eq. (4.9b) gives a final condition at *y* = *y*_{0}. Hence, Eq. (4.9a) is to be solved “backwards” from *y* = *y*_{0} down to *y* = 0. Notice that this adjoint problem is the same as the the one-way RTE for *I ^{−}*.

## 5. Numerical method

In what follows, we give numerical methods to solve the initial value problem for the non-homogeneous one-way RTE,

In doing so, we will be able to describe a numerical method to solve boundary value problem (4.2) for the radiative transfer equation. The numerical methods involve the following components.

- Discrete ordinate method for
*θ*using an odd-ordered Gauss-Legendre quadrature rule mapped to 0 <*θ*<*π*. - Centered, second-order finite differences for
*x*with periodic boundary conditions. - Crank-Nicolson to advance the solution in
*y*.

In what follows, we describe each of the items above to solve initial value problem (5.1) and final value problem (5.2). Then we describe an iterative method to solve boundary value problem (4.2).

#### 5.1. Discrete ordinate method

The 2*M* + 1 point Gauss-Legendre quadrature rule takes the form

*x*are the quadrature abscissae and

_{m}*w*are the quadrature weights. Let

_{m}*θ*=

_{m}*π*(

*x*+ 1)/2. We use an odd number of points for this quadrature rule to ensure that

_{m}*θ*

_{M/}_{2}=

*π/*2 is a quadrature abscissa identically. Having this quadrature abscissa is convenient for this particular study involving the imaging system depicted in Fig. 1 since the plane wave is incident normally on the

*y*= 0 boundary corresponding to

*θ*=

*π/*2. We now introduce

For the discrete ordinate method, we replace
${\mathcal{P}}_{f}$ in Eqs. (5.1a) and (5.2a) by
${P}_{f}^{{}^{2M+1}}$ given in Eq. (5.4), and evaluate those resulting equations at the quadrature abscissae, *θ _{m}* for

*m*= 1, ⋯, 2

*M*+1. Let

*I*(

_{m}*x,y*)

*≈ I*(

*θ*) and

_{m},x,y*Ĩ*(

_{m}*x,y*)

*≈ I*(

*θ*) denote the approximate values of

_{m},x,y*I*and

*Ĩ*, respectively. Applying the discrete ordinate method to initial value problem (5.1) yields

*m*= 1

*,·*⋯,2

*M*+ 1. For final value problem (5.2), we obtain

*m*= 1

*,·*⋯,2

*M*+ 1.

#### 5.2. Finite differences

For the interval −*L _{x}*/2

*≤ x ≤ L*/2, we introduce the equi-spaced grid

_{x}*x*= −

_{i}*L*/2 + (

_{x}*i*−1)

*h*for

*i*= 1, ⋯

*,N*with

*h*=

*L*denoting the mesh width. Let

_{x}/N*I*(

_{m},i≈ I_{m}*x*) and

_{i},y*Ĩ*(

_{m},i ≈ Ĩm*x*) denote the approximations of

_{i},y*I*and

_{m}*Ĩ*, respectively, evaluated at grid point

_{m}*x*. We apply periodic boundary conditions so that

_{i}*I*

_{m,}_{0}=

*I*and

_{m,N}*Ĩ*

_{m,}_{0}=

*Ĩ*. Consequently,

_{m,N}*I*

_{m,N}_{+1}=

*I*

_{m,}_{1}and

*Ĩ*

_{m,N}_{+1}=

*Ĩ*

_{m,}_{1.}Using these numerical boundary conditions, we apply the centered, second order finite difference approximation in

*x*for initial value problem (5.5) and obtain

*m*= 1

*,·*⋯,2

*M*+ 1 and

*i*= 1

*,·*⋯

*,N*. Because

*μ*and

_{a}*μ*are functions of

_{s}*x*and

*y*, we have added subscripted indices to them to indicate that we are to evaluate them at grid point

*x*in Eq. (5.7a). Similarly, for final value problem (5.6), we obtain

_{i}*m*= 1

*,·*⋯,2

*M*+ 1 and

*i*= 1, ⋯

*,N*.

We now introduce some notation*·*that will allow us to proceed with the numerical method more easily. Let
${D}_{x}^{N}$ be the *N × N* matrix defined as

*T*and

_{c}*T*denote the following (2

_{s}*M*+ 1)

*×*(2

*M*+ 1) diagonal matrices:

Furthermore, let *A*(*y*) and *S*(*y*) denote the *N × N* diagonal matrices

The operator
${P}_{f}^{{}^{2M+1}}$, defined in Eq. (5.4), is a (2*M* + 1) *×* (2*M* + 1) matrix with entries [
${P}_{f}^{{}^{2M+1}}$ ]_{m,m}_{′} =
$\frac{\pi}{2}$*p*(*θ _{m} −θ_{m}*

_{′})

*w*

_{m}_{′}. We now rewrite initial value problem (5.7) for the (2

*M*+1)

*N*-vector

**I**(

*y*) = [

*I*

_{1}

_{,}_{1}(

*y*)

*,I*

_{2}

_{,}_{1}(

*y*)

*,·*⋯

*, I*

_{2}

_{M}_{+1}

_{,}_{1}(

*y*)

*, I*

_{1}

_{,}_{2}(

*y*)

*,·*⋯

*, I*

_{2}

_{M}_{+1}

*(*

_{,N}*y*)] and obtain

The (2*M* + 1)*N ×* (2*M* + 1)*N* matrices, *B* and *C*, are defined as

Here,
${\mathrm{I}}_{2M+1}$ and
${\mathrm{I}}_{N}$ are the (2*M* +1)*×*(2*M* +1) and *N×N* identity matrices, respectively, and ⊗ denotes the Kronecker product. Similarly, we can rewrite final value problem (5.8) in terms of the (2*M* + 1)*N* vector Ĩ as

#### 5.3. Crank-Nicolson

The Gauss-Legendre quadrature rule used above is an open rule, so the quadrature abscissas do not include the end points. Nonetheless, because sin*θ →* 0 as *θ →* 0*,π*, we find that the matrix *B* has very small diagonal entries. The contrast between these very small entries and the others for which sin*θ ∼* 1 cause initial value problem (5.14) and final value problem (5.17) to be stiff. Consequently, we apply the Crank-Nicolson method, which is an implicit, second order, and unconditionally stable method, to solve these problems.

Let *y _{j}* =

*j*Δ

*y*for

*j*= 0,⋯

*,J*with Δ

*y*=

*y*

_{0}

*/J*. Furthermore, let

**I**

^{j}≈**I**(

*y*) and Ĩ

_{j}*Ĩ(*

^{j}≈*y*) for

_{j}*j*= 1, ⋯

*,J*. Applying the Crank-Nicolson method to initial value problem (5.14), we obtain

**I**

^{0}=

**f**. The (2

*M*+1)

*N ×*(2

*M*+1)

*N*linear system given in Eq. (5.19) is to be solved for

*j*increasing starting at

*j*= 1 up to

*j*=

*J*. Similarly, applying the Crank-Nicolson method to final value problem (5.17), we obtain

*= $\tilde{\mathbf{f}}$. Equation (5.20) is to be solved for*

^{J}*j*decreasing starting at

*j*=

*J*down to

*j*= 1.

The matrices in the linear systems above are sparse. Exploiting their sparsity reduces the demand for computational resources and allows for more computationally efficient solutions.

#### 5.4. Iterative method to solve the radiative transfer equation

We can use the numerical method described above to solve the boundary value problem (4.2) for the RTE. To do so requires iterating forward and backward one-way RTEs for *I*^{+} and *I ^{−}*, respectively. We use this iterative method to compute simulated data from which we reconstruct quantitative images.

We denote the *k*th iteration of *I ^{±}* by

*I*

^{(}

^{k}^{)}

*. Let ${P}_{b}^{{}^{2M+1}}$ be the (2*

^{±}*M*+ 1)

*×*(2

*M*+ 1) matrix with entries [ ${P}_{b}^{{}^{2M+1}}$]

_{m,m}_{′}= $\frac{\pi}{2}$

*p*(

*θ*

_{m}−θ_{m}_{′}

*−π/*2)

*w*

_{m}_{′}. We set

*k*= 1 and

*I*

^{(0)}

*= 0. Then, we start the following iterative procedure.*

^{−}- Solve the initial value problem$$\begin{array}{c}\mathrm{sin}{\theta}_{m}\frac{\partial {I}^{(k)+}}{\partial y}+\mathrm{cos}\theta \frac{\partial {I}^{(k)+}}{\partial x}+\left({\mu}_{t}+\mathrm{i}\frac{\omega}{c}\right){I}^{(k)+}-{\mu}_{s}{\mathcal{P}}_{f}{I}^{(k)+}={\mu}_{s}{\mathcal{P}}_{b}{I}^{(k-1)-},\\ {I}^{(k)+}(\theta ,x,0)={I}_{0}\delta (\theta -\pi /2),\end{array}$$using the numerical method given in Eq. (5.19).
- Solve the final value problem$$\begin{array}{c}-\mathrm{sin}{\theta}_{m}\frac{\partial {I}^{(k)-}}{\partial y}-\mathrm{cos}\theta \frac{\partial {I}^{(k)-}}{\partial x}+\left({\mu}_{t}+\mathrm{i}\frac{\omega}{c}\right){I}^{(k)-}-{\mu}_{s}{\mathcal{P}}_{f}{I}^{(k)-}={\mu}_{s}{\mathcal{P}}_{b}{I}^{(k)+},\\ {I}^{(k)-}(\theta ,x,{y}_{0})=0,\end{array}$$using the numerical method given in Eq. (5.20).
- Check convergence. If convergence has not been reached, set
*k ← k*+ 1, and repeat. We compute the transmittance*T*(*x*) defined in Eq. (4.5) using the numerical approximation

The sum in Eq. (5.21) is over all *θ _{m}* values within

*NA*. Let

**t**

^{(}

^{k}^{)}denote the vector whose entries are given by

*T*

^{(}

^{k}^{)}(

*x*) defined in Eq. (5.21). To check for convergence, we check to see if $\sqrt{h}{\Vert {\mathbf{t}}^{(k)}-{\mathbf{t}}^{(k-1)}\Vert}_{2}<\epsilon $ for convergence where

_{i}*ε*is some user-defined tolerance value. In the results we show below in Section 7, we have set

*ε*= 10

^{−}^{8}.

In contrast to the conventional source iteration method [15], this iteration scheme is the same as the improved source-iteration introduced by Gao and Zhao [16]. The one-way RTE is the leading order approximation associated with this iteration method.

## 6. Inverse problem

We now have all of the components in place to specify and solve the inverse problem for the imaging system described in Section 2.

To simulate measured data, we solve the full RTE using the iterative method described in Section 5.4. In doing so, we obtain a collection of transmittance measurements for each orientation of the sample as it is rotated. Let **d** denote the vector of measured data whose components are *d _{i j}* =

*T*

_{RTE}(

*x*;

_{i}*φ*), the transmittance computed on the detector plane using the RTE sampled at the points

_{j}*x*, for

_{i}*i*= 1, ⋯

*,N*of the finite differences grid introduced in Section 5.2., for each of the

*A*rotations of the sample defined as

_{s}*φ*= 2

_{j}*π*(

*j −*1)

*/A*for

_{s}*j*= 1, ⋯

*,A*. Therefore,

_{s}**d**contains

*A*transmittance measurements. For the numerical results shown in Section 7, we have set

_{s}×N*L*= 12.8 and

_{x}*y*

_{0}= 5 and used

*M*= 16 for the quadrature rule,

*N*= 128 for the

*x*-grid, and

*J*= 50 for the

*y*-grid so that Δ

*x*= Δ

*y*= 0.1. In addition, we have used

*A*= 25 rotations for measurements. Also, we have set the source’s modulation frequency to

_{s}*ω*= 600MHz.

To model this measured data, we solve the one-way RTE, Eq. (4.8) using the method described in Sections 5.1 – 5.3 using the same grid used to compute the measurements. Even though we have used the same discretizations for the measurements and forward model, we are solving two completely different problems for the measurements and the forward model. The corresponding modeling error dominates over truncation error of the numerical methods used. Consequently, we do not anticipate artificially optimistic results in what follows. Let **m** denote the vector whose components are *m _{i j}* =

*T*

_{one-way}(

*x*;

_{i}*φ*), the transmittance computed on the detector plane using the one-way RTE sampled at the same points as for the measured data. To make explicit the dependence of

_{j}**m**on the absorption and scattering coefficients, we write the model as

**m**(

*μ*).

_{a},μ_{s}The inverse problem is as follows.

*Reconstruct the absorption and scattering coefficients of the sample given the measured data*, **d***, and the model*, **m**(*μ _{a},μ_{s}*).

Let

be the cost functional for the residual*R*, defined as

We compute reconstructions of the absorption and scattering coefficients by seeking the distributions of *μ _{a}* and

*μ*that minimize the cost functional given in Eq. (6.1). To minimize the cost functional, we start with an initial guess and update that guess along a descent direction of

_{s}*K*until we approach the minimum value within some tolerance. This descent direction is found by writing the parameters,

*μ*and

_{a}*μ*, as perturbations to their background values, ${\overline{\mu}}_{a}$ and ${\overline{\mu}}_{s}$:

_{s}Doing so leads to the following linearized residual operator,

*R*′ is the Fréchet derivative of

*R*. By setting the right hand side of Eq. (6.5) to zero, we obtain a set of equations for

*δ μ*and

_{a}*δ μ*. The minimum norm solution of these resulting equations is given by

_{s}Here
$\tilde{{R}^{\prime}}({\overline{\mu}}_{a},{\overline{\mu}}_{s})$ denotes the adjoint of
${R}^{\prime}({\overline{\mu}}_{a},{\overline{\mu}}_{s})$. It can be proved that Eq. (6.6) is a descent direction of the cost functional *K* (see Appendix A in [13]). However, the functional
${\left[{R}^{\prime}({\overline{\mu}}_{a},{\overline{\mu}}_{s})\tilde{{R}^{\prime}}({\overline{\mu}}_{a},{\overline{\mu}}_{s})\right]}^{-1}$ is very time consuming to calculate. Instead, we typically replace it by a multiple of the identity operator [13] to obtain

*K*.

To compute Eq. (6.7) we first solve Eq. (4.8) with optical parameters
${\overline{\mu}}_{a}$ and
${\overline{\mu}}_{s}$, subject to initial condition (4.2c) to obtain *I*^{+}. Then, we use that result to compute the residual *ρ* =
$R({\overline{\mu}}_{a},{\overline{\mu}}_{s})$ through comparison with the measured data. To compute
$\tilde{{R}^{\prime}}({\overline{\mu}}_{a},{\overline{\mu}}_{s})$*ρ*, we solve Eq. (4.9) using *ρ* in the final condition in Eq. (4.9b). In particular, the final condition prescribes that *Ĩ*(*θ,x,y*_{0}) is equal to the complex conjugate of *ρ* distributed uniformly over the set of directions making up the numerical aperture, *NA*, of the detectors, and zero for all other directions.

Upon solution of this final value problem, we obtain *Ĩ*^{+}. Having calculated *I*^{+} and *Ĩ*^{+}, we then compute the updates

With the updates given by Eqs. (6.8) and (6.9) defined, we now give the procedure for reconstructing images of *μ _{a}* and

*μ*.

_{s}For each *φ _{j}* with

*j*= 1,⋯

*,A*, we perform the following reconstruction procedure with ${\overline{\mu}}_{a}$ and ${\overline{\mu}}_{s}$ denoting the current absorption and scattering coefficients, respectively.

_{s}- Solve the initial value problem$$\begin{array}{c}\mathrm{sin}\theta \frac{\partial {I}^{+}}{\partial y}+\mathrm{cos}\theta \frac{\partial {I}^{+}}{\partial x}+\left({\overline{\mu}}_{a}+{\overline{\mu}}_{s}+\mathrm{i}\frac{\omega}{c}\right){I}^{+}-{\overline{\mu}}_{s}{\mathcal{P}}_{f}{I}^{+}=0,\phantom{\rule{1em}{0ex}}\text{in}\phantom{\rule{0.2em}{0ex}}0<y\le {y}_{0}\\ {I}^{+}=(\theta ,x,0)={I}_{0}\delta (\theta -\pi /2)\phantom{\rule{1em}{0ex}}0<\theta <\pi ,\end{array}$$to obtain
*I*^{+}and use it to model the measured transmittance by computing*T*one-way(*x*;_{i}*φ*) for_{j}*i*= 1, ⋯*,N*._{s} - Solve the adjoint problem$$\begin{array}{c}-\mathrm{sin}\theta \frac{\partial {\tilde{I}}^{+}}{\partial y}-\mathrm{cos}\theta \frac{\partial {\tilde{I}}^{+}}{\partial y}+\left({\overline{\mu}}_{a}+{\overline{\mu}}_{s}+\mathrm{i}\frac{\omega}{c}\right){\tilde{I}}^{+}-{\overline{\mu}}_{s}{\mathcal{P}}_{f}{\tilde{I}}^{+}=0,\phantom{\rule{1em}{0ex}}\text{in}0\le y<{y}_{0}\\ {I}^{+}(\theta ,x,{y}_{0})=\{\begin{array}{ll}{\rho}^{*}/\mathrm{\Delta}\theta & \text{on}-{x}_{0}/2\le x\le {x}_{0}/2,and(\pi -\mathrm{\Delta}\theta )/2\le \theta \le (\pi +\mathrm{\Delta}\theta )/2,\\ 0& \text{otherwise},\end{array}\end{array}$$to obtain
*Ĩ*. Here, Δ^{+}*θ*sets the angular width of the numerical aperture,*NA*, about*θ*=*π/*2, and*ρ*denotes the complex conjugate of^{∗}*ρ*.

We repeat this procedure until the value of the cost functional given in Eq. (6.1) reduces below a user-defined tolerance value. We initialize this reconstruction procedure by setting ${\overline{\mu}}_{a}$ and ${\overline{\mu}}_{s}$ to the known optical properties of the background medium in the disk.

## 7. Numerical results

We show reconstructions computed using the procedure described above in Section 6. For all of the results shown here, the absorption and scattering coefficients of the slab background are *μ _{a}* = 0.01cm

^{−}^{1}and

*μ*= 1.0cm

_{s}

^{−}^{1}, respectively. The disk containing the sample to be reconstructed has radius

*r*= 2cm. Within that disk the background absorption and scattering coefficients are

*μ*= 0.015cm

_{a}

^{−}^{1}and

*μ*= 1.1cm

_{s}

^{−}^{1}, respectively. The anisotropy factor is set to

*g*= 0.9.

Figure 3 shows reconstructions of a sample containing one absorber and one scatterer. Figure 3(a) shows the true image of *μ _{a}*, and Fig. 3(b) shows the true image of

*μ*for a fixed orientation of the sample. The absorption coefficient inside the absorber is

_{s}*μ*= 0.1cm

_{a}

^{−}^{1}and the scattering coefficient inside the scatterer is

*μ*= 1.5cm

_{s}

^{−}^{1}. Figure 3(c) shows the reconstructed image of

*μ*and Fig. 3(d) shows the reconstructed image of

_{a}*μ*. Although there are artifacts in both reconstructions, we find that they clearly identify the absorber and scatterer. The artifacts are more apparent in the reconstruction of the absorption coefficient. The quantitative accuracy of these reconstructions are good. In particular, the root mean square (RMS) error for the scattering coefficient reconstruction is 2%. The RMS error for the absorption coefficient reconstruction is 27% which is much higher. The reconstruction of the absorption coefficient suffers from cross-talk with the reconstruction of the scattering coefficient. Cross-talk in reconstructed images is typical in DOT [13]. It is not necessarily a feature associated with using the one-way RTE as the forward model. This cross-talk contributes to the higher RMS error for the reconstruction of the absorption coefficient. Nonetheless, we find that the one-way RTE provides an effective forward model to compute reconstructions of the absorption and scattering coefficients.

_{s}In Fig. 4, we show reconstructions for a sample containing two absorbers and one scatterer, each situated in a distinct region in the sample. The smaller of the two absorbers which appears near the top of Fig. 4(a) has absorption coefficient *μ _{a}* = 0.115cm

^{−}^{1}. The larger of the two absorbers which appears near the bottom of Fig. 4(b) has absorption coefficient

*μ*= 0.135cm

_{a}

^{−}^{1}. The scatterer shown in Fig. 4(b) has scattering coefficient

*μ*= 1.4cm

_{s}

^{−}^{1}. Figure 4(c) shows the reconstruction of

*μ*and Fig. 4(d) shows the reconstruction of

_{a}*μ*. For this more challenging case, we find that the artifacts in the reconstruction of

_{s}*μ*interfere more with the reconstruction. There still remains some cross-talk. Nonetheless, we can identify two absorbers and observe some quantitative differences between them with the lower absorber having larger values. The RMS error for the reconstruction of the scattering coefficient is 2% and for the reconstruction of the absorption coefficient is 21%.

_{a}In Fig. 5, we show reconstructions for two absorbers and one scatterer situated with overlapping regions in the sample. The absorption and scattering coefficients in each of the absorbers and scatterer are the same as those shown in Fig. 4. We have just changed their locations so that there is a non-trivial overlapping region. For this very challenging case, we obtain reasonable images. The scatterer can be identified and its quantitative accuracy is reasonable given the severe image artifacts. However, we have lost both qualitative and quantitative accuracy for the reconstruction of the absorbers. The RMS error for the reconstruction of the scattering coefficient is 6% and for the reconstruction of the absorption coefficient is 30%.

To attempt to reduce the severe image artifacts appearing Fig. 5, we use *L*_{2}-regularization using the heat kernel, as explained in [14]. Figure 6 shows reconstructions for the same sample as in Fig. 5, but using the regularization mentioned above. This regularization has improved the qualitative and quantitative results, especially for the reconstruction of the absorption coefficient. Some of the image artifacts seen in Fig. 5 have been suppressed. However, cross-talk with the reconstruction of the scattering coefficient still remains. The RMS error for the reconstruction of the scattering coefficient is 2% and for the reconstruction of the absorption coefficient is 32%. It is possible that using other regularization techniques may be more useful here. For example, using total variation (TV) regularization may be useful for segmenting distinct regions of the reconstructed images.

## 8. Conclusions

For transmission problems, the one-way RTE approximates the RTE when scattering is anisotropically forward-peaked. It is derived by neglecting backscattering thereby reducing the boundary value problem for the RTE to an initial value problem for the one-way RTE. This initial value problem is physically intuitive, accurate, and much more efficient to solve than the full boundary value problem.

The imaging system studied here is slightly idealized in that boundaries are index-matched. More practical systems include index-mismatched boundaries. These boundaries reduce the accuracy of the one-way RTE, but early results have shown that the one-way RTE is still effective for index-mismatched boundaries [12]. The reason for this is as follows. At the boundary where light is transmitted, partial reflection due to the index-mismatched boundary will redistribute a fraction of the forward flowing power backwards. For this light to eventually contribute to transmission, it must backscatter. But the argument in deriving the one-way RTE is that back-scattering is negligible due to anisotropic forward-peaked scattering. Thus, the contribution of light partially reflected at boundaries due to index-mismatch is a small, higher order correction to the one-way RTE.

Extending this work to three dimensions is conceptually straight-forward. For this case, the RTE is different because the scattering phase function is different. Nonetheless, the one-way RTE approximation is derived in exactly the same way shown here. One simply neglects the hemisphere of directions corresponding to backscattering and obtains an initial value problem. In addition, the numerical method used to solve the RTE and the one-way RTE can be readily extended for three dimensional problems. The main challenge associated with extending this work to three dimensions is that the size of the numerical computations increases significantly and may require modifications to the computational methods used. For example, the linear system given in Eq. (5.19) is small enough that it is solved directly. For the three dimensional problem, directly solving the corresponding linear system will likely be prohibitively expensive. Instead, one will need to use an iterative method along with an appropriate pre-conditioner that exploits the inherent sparsity of the system to obtain greater efficiency.

Using this initial value problem for the one-way RTE to compute a forward model, we have reconstructed quantitatively accurate images. Since the one-way RTE is a simplification of the RTE, it does not rectify issues in the inverse problem inherent also to the RTE such as cross-talk and the need for regularization. Moreover, its resolution is the same or worse than that inherent to the RTE. Nonetheless, we conclude from these results that the one-way RTE is an effective and efficient approximation for use in DOT problems.

## References and links

**1. **A. Ishimaru, *Wave Propagation and Scattering in Random Media* (IEEE Press, 1996).

**2. **L. V. Wang and H.-I. Wu, *Biomedical Optics* (Wiley, 2007).

**3. **S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. **15**, R41–R93 (1999). [CrossRef]

**4. **S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Probl. **25**, 123010 (2009). [CrossRef]

**5. **P. González-Rodríguez, A. Kim, and M. Moscoso, “Robust depth selectivity in mesoscopic scattering regimes using angle-resolved measurements,” Opt. Lett. **38**, 787–789 (2013). [CrossRef] [PubMed]

**6. **O. Dorn, “A transport-backtransport method for optical tomography,” Inverse Probl. **14**, 1107–1130 (1998). [CrossRef]

**7. **A. D. Klose and A. H. Hielscher, “Iterative reconstruction scheme for optical tomography based on the equation of radiative transfer,” Med. Phys. **26**, 1698–1707 (1999). [CrossRef] [PubMed]

**8. **K. Ren, G. Bal, and A. H. Hielscher, “Frequency Domain Optical Tomography Based on the Equation of Radiative Transfer,” SIAM J. Sci. Comput. **28**, 1463–1489 (2006). [CrossRef]

**9. **G. Bal, “Inverse transport theory and applications,” Inverse Probl. **25**, 053001 (2009). [CrossRef]

**10. **H. K. Kim and A. H. Hielscher, “A diffusion-transport hybrid method for accelerating optical tomography,” J. Innov. Opt. Health Sci. **3**, 1–13 (2010). [CrossRef]

**11. **T. Tarvainen, V. Kolehmainen, S. R. Arrdige, and J. P. Kaipio, “Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model,” J. Quant. Spectrosc. Radiat. Transf. **112**, 2600–2608 (2011). [CrossRef]

**12. **P. González-Rodríguez, B. Ilan, and A. D. Kim, “The one-way radiative transfer equation,” submitted for publication.

**13. **P. González-Rodríguez and A. Kim, “Comparison of light scattering models for diffuse optical tomography,” Opt. Express **17**, 8756–8774 (2009). [CrossRef] [PubMed]

**14. **P. González-Rodríguez, O. Dorn, M. Kindelan, and M. Moscoso, “History matching problem in reservoir engineering using the propagation-backpropagation method,” Inverse Probl. **21**, 565–590 (2005). [CrossRef]

**15. **E. E. Lewis and W. F. Miller Jr., *Computational Methods of Neutron Transport* (American Nuclear Society, La Grange Park, IL, 1993)

**16. **H. Gao and H. Zhao, “A fast forward solver of radiative transfer equation,” Transport Theory Statist. Phys. **38**, 149–192 (2009). [CrossRef]