Computational and Mathematical Methods in Medicine

Volume 2015 (2015), Article ID 286161, 23 pages

http://dx.doi.org/10.1155/2015/286161

## Image Reconstruction for Diffuse Optical Tomography Based on Radiative Transfer Equation

^{1}Department of Mathematics, Harbin Institute of Technology, Harbin, Heilongjiang 150006, China^{2}School of Mathematics and Statistics, Northeast Petroleum University, Daqing, Heilongjiang 163318, China^{3}Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA^{4}Imaging Diagnosis and Interventional Center, State Key Laboratory of Oncology in South China, Sun Yat-sen University Cancer Center, Guangzhou, Guangdong 510060, China

Received 5 October 2014; Revised 8 December 2014; Accepted 17 December 2014

Academic Editor: Reinoud Maex

Copyright © 2015 Bo Bi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Diffuse optical tomography is a novel molecular imaging technology for small animal studies. Most known reconstruction methods use the diffusion equation (DA) as forward model, although the validation of DA breaks down in certain situations. In this work, we use the radiative transfer equation as forward model which provides an accurate description of the light propagation within biological media and investigate the potential of sparsity constraints in solving the diffuse optical tomography inverse problem. The feasibility of the sparsity reconstruction approach is evaluated by boundary angular-averaged measurement data and internal angular-averaged measurement data. Simulation results demonstrate that in most of the test cases the reconstructions with sparsity regularization are both qualitatively and quantitatively more reliable than those with standard regularization. Results also show the competitive performance of the split Bregman algorithm for the DOT image reconstruction with sparsity regularization compared with other existing algorithms.

#### 1. Introduction

Diffuse optical tomography (DOT) is an emerging imaging modality that has attracted much attention in clinical diagnosis, for example, in breast cancer detection, monitoring of infant brain tissue oxygenation level, and functional brain activation studies, cerebral hemodynamic, and so forth; compare [1–3]. With the use of near-infrared (NIR) light, DOT probes the optical properties, mainly the absorption coefficient and the scattering coefficient of human tissues. In experimental systems, a set of optical fibres and optodes are attached to the boundary of the object as measurement detectors and sources. NIR light as the inflow current is emitted by laser and guided by some fibre optics into the object, one position at a time. The light is transmitted, and then the outflow current is measured from all the measurement positions using light sensitive detectors.

By employing a setup comprising a set of external light sources and detectors, the optical properties of a tissue can be recovered by applying the principles of tomography. The difference in absorption or scattering between the normal and abnormal tissues can provide the imaging contrast for tissue diagnostics. In this context, the aim of diffuse optical tomography is to provide the spatial distribution of the absorption and scattering coefficients of a tissue.

The forward problem in DOT describes the photon propagation in tissues and the inverse problem involves estimating the absorption and scattering coefficients of tissues from light measurements on the surface. There are many models describing the light propagation, for example, the Fokker-Planck equation, differential approximations of the radiative transfer equation (RTE) and so forth (cf. [4, 5]), and the diffusion equation (DA), which is accepted as a popular forward model describing light propagation in tissues. DA is a low-order approximation of RTE. However, the validation of DA breaks down in the following situations. First, since DA is only valid in the tissues with high scattering and low absorption, while in many tissues of human body with low scattering, such as skeleton, joint, DA is largely limited. Second, DA is not suitable to describe the light propagation in regions near to the light source, which will result in the large error when using DA as the forward model to image small animal or small tissue. For more on DA, we refer to [1, 6, 7]. Due to these limitations of DA, RTE is a more natural choice to model the light propagation. In this paper, we will present numerical evidence to show that DOT based on RTE can achieve satisfactory resolution.

Since DOT suffers from severe ill-posedness caused by noise and incomplete measurement data, its efficient, stable, and accurate numerical treatment is very challenging. Due to its immense range of prospective applications, there has been a vast amount of research work on mathematical as well as practical aspects of the inverse problem. In particular, the design of efficient and stable numerical algorithms has received considerable attention.

In the linearized Jacobian-based methods, a popular approach for solving the minimization problem is to treat it as a nonlinear least-square problem that can be solved with many standard optimization techniques [8], among which Levenberg-Marquardt method (LM) [9] is perhaps the most commonly used one, which is a Gauss-Newton method with regularization. Based on RTE, the authors in [10] discussed in detail the LM method for parameter estimation problems. We will describe the LM method briefly in Section 3. In the nonlinear gradient-based methods, we directly minimize the nonlinear functional; in the optimization process only the gradient of the nonlinear functional needs to be computed, for example, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) method and the L-BFGS method [8, 11].

The appropriate choice of regularization depends on* a priori* knowledge of the domain and the inclusions. Although regularization is usually a natural choice for its simplicity, it is not the optimal strategy. For example, when the coefficient distribution is sparse or discontinuous, it is well known that the regularization method [12, 13] or the bounded total variation (TV) regularization method [14, 15] is more efficient than regularization.

Sparse reconstruction has attracted much researchers’ attention, especially since Donoho et al. [16, 17] established the theory of compressive sensing (CS). CS describes the sparse reconstruction problem as an quasinorm optimization [18]. However, Donoho in [19] proved that the regularization can also obtain the sparsest solution and proposed the equivalent condition between the regularization and the regularization. In [20], Candes et al. studied the stable signal recovery from incomplete and inaccurate measurements and reduced a computationally difficult problem to the basis pursuit problem. Because the sparsity regularization is nonsmooth, it is still a challenge to find efficient methods to solve this convex basis pursuit optimization problem, and the choice of techniques for solving it becomes crucial. Classical gradient-based methods usually bring high computational burden [21]. Motivated by the inverse problem in imaging [22–25], the authors used and developed the Bregman iteration technique for the regularization problems and proved that the Bregman iteration method is an effective way to solve the norm minimization problems.

In this paper, we adopt the split Bregman method for solving sparsity regularization problems. The split Bregman method is a simple and efficient algorithm which can split the minimization into and functionals and significantly reduce the computational burden [24]. We will introduce the split Bregman algorithm in detail in Section 3.

The organization of this paper is as follows. Section 2 introduces the DOT forward problem, briefly presents the radiative transfer equation (RTE), and describes the imaging modality of the forward problem. Section 3 presents in detail the reconstruction algorithm, including the Levenberg-Marquardt algorithm and the split Bregman algorithm. Section 4 provides several simulations to show the validity of the sparsity regularization reconstruction. Section 5 gives some conclusion statements.

#### 2. Forward Problem

In this section, we formulate our problem of RTE based diffuse optical tomography.

##### 2.1. Radiative Transfer Equation

Photon propagation in tissues can be described by the radiative transfer equation. Let , or , denote the physical domain of the medium with boundary , the unit sphere, the unit outer normal vector, and the outgoing and incoming boundaries defined by

The variables and denote the spatial position and the angular direction. Then we consider the following boundary value problem (BVP) of RTE: where is the total attenuation coefficient, describes the probability that a photon is absorbed in unit length, its reciprocal being the absorption mean free path, and is the scattering coefficient, describing the probability that a photon is scattered in unit length, its reciprocal being the scattering mean free path. Further, is the radiance and is the internal light source. In this paper, we consider the case with no light source inside ; . is the inflow current on , and the boundary condition (3) implies that once a photon escapes the domain , it does not reenter it; compare [26]. is the scattering operator. Denote by the infinitesimal area element on the unit sphere . Then is defined as with a nonnegative normalized phase function: In many applications, the function is independent of . However, in our general study, we allow to depend on . Indeed, we can even allow to be a general function of , , and , that is, in the form . The scattering phase function describes the probability that a photon with an initial direction will have a direction after a scattering event. In DOT, one typical example is the Henyey-Greenstein phase function; compare [27]. Consider where the parameter is the anisotropy factor of the scattering medium. Note that for isotropic scattering, for forward scattering, and for backward scattering. In biomedical imaging problems, the scattering is strongly forward peaked and is close to 1.

##### 2.2. Forward Problem

Similar to the X-ray CT, DOT experiments acquire the current distribution of detectors on the boundary under the multi-incidents [28]. The experimental procedure of acquiring potential measurements is as follows. First, set a set of laser devices and detectors on the boundary of the object. Then launch an incident impulse from one laser device and record the resulting measurements from all the detectors. In order to gather enough data information, repeat this procedure on other laser devices.

We can model this procedure mathematically. Before doing this, let us introduce some assumptions on the domain and the coefficients, assumed to be valid throughout the rest of this paper.

We assume that the absorption and scattering coefficients are approximated piecewise-constant. It is known that the BVP (2) and (3) has a unique solution [29], and the solution depends continuously on the input current on the incoming boundary .

The forward problem of DOT is to determine the outgoing current on the detectors when the incident impulse and the absorption and scattering coefficients are known. Excite the domain with a sequence of incident impulses , , and get a sequence of measurements corresponding to each incident impulse, with its component , , being the measurement value on detectors. Then a mathematical description of such an experiment is provided by a sequence of forward operators: which maps prescribed optical parameters to the corresponding measurements. Here, denotes the th forward operator with respect to the th incident impulse and the resulting detected measurement data on detectors. Then, for , the domain of the operator is defined as is a column vector representing the measurement data on detectors.

We mention that, in DOT, the measured quantity is the excitance on the boundary of the domain. Due to the limit of measurement techniques of the optical devices, the angularly resolved measurement data is not practical. In this study, for , we use the boundary angularly averaged data as our measurements [30]:
where is the solution of the BVP (2) and (3) corresponding to the optical coefficients and the th incident impulse and is the unit outer normal vector, whereas the set
In [30], a detailed analysis is given on properties of the forward operator, including the Lipschitz continuity and the* Fréchet* differentiability. The BVP (2) and (3) needs to be solved by numerical methods, such as the finite difference method or the finite element method in spatial discretization, and the finite element method or approximation in angular discretization. For the simulation results in this paper, we use the RTE2DMATLAB codes [6], in which the finite element method in both spatial and angular spaces is used to discretize the forward mapping. We refer the reader to [4, 6, 31–34] for details about the finite element implementation in both spatial and angular space.

#### 3. Inverse Problems

In practice, due to the limitations of the experimental environment and the laboratory equipment, the measurement we get usually contains noise. Here, we assume that, for , the actual measurements have the noise level ; that is, , where represents the actual measurement data and represents the true data corresponding to the true optical coefficients. Then our inverse problem of DOT is to determine such that the following nonlinear equations hold: for .

In this paper, we only reconstruct the scattering coefficient, while assuming that the distribution of the total attenuation coefficient is known.

As is typical for many inverse problems, the DOT inverse problem is ill-posed. In order to reconstruct the optical parameter stably, regularization is required. We minimize the following Tikhonov functional:
over the admissible set
for the coefficient . Here, is a regularization penalty functional that enforces* a priori* knowledge on the optical parameter to be reconstructed and is the regularization parameter used to trade off the discrepancy term (the first item of ) that incorporates the information contained in the data and [35]. Then we analyse the minimization problem:

We consider two popular regularization formulations: the norm penalty and the sparsity constraint regularization. We will describe these two approaches below.

##### 3.1. Standard Reconstruction

First, we use the traditional norm squared penalty, which consists of minimizing the following functional:
with a specified based on* a priori* information on the solution. In the statistical inversion framework, the corresponding prior constructions are known as smoothness priors [36]. To compute a minimizer of problem (14), many iterative regularization methods are available. We use the Levenberg-Marquardt method based on linearization. Specifically, for every , the forward operator is linearized around some initial guess ; that is,
where is the* Fréchet* derivative of with respect to the coefficient at and denotes the Taylor remainder for the linearization around . Then substituting the above linearized expression into the functional and ignoring the higher-order remainder , we get a linearized problem
The Euler equation of the discrete problem is
that is,
with the identity matrix. The system can be solved directly to get a new estimate for based on the initial guess . Then we iteratively update the reconstruction by taking the solution as an initial guess. In practice, the iterative procedure achieves required accuracy within a few iterations. The complete algorithm is given in Algorithm 1. The stopping criterion can be defined based on monitoring the relative change of consecutive iterations.