Table of Contents Author Guidelines Submit a Manuscript
Computational and Mathematical Methods in Medicine
Volume 2015 (2015), Article ID 286161, 23 pages
http://dx.doi.org/10.1155/2015/286161
Research Article

Image Reconstruction for Diffuse Optical Tomography Based on Radiative Transfer Equation

1Department of Mathematics, Harbin Institute of Technology, Harbin, Heilongjiang 150006, China
2School of Mathematics and Statistics, Northeast Petroleum University, Daqing, Heilongjiang 163318, China
3Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA
4Imaging 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 [13]. 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 [2225], 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, 3134] 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.

Algorithm 1: Reconstruction algorithm based on the regularizing Levenberg-Marquardt method.

3.2. Sparsity Reconstruction

In the sparsity reconstruction, the functional to be minimized is of the form We will apply the Bregman framework to solve (20). The key to our method is to “decouple” the and portions in (20). Rather than (20), we will consider the constrained problem [23]: To solve the above minimization problem, the corresponding unconstrained optimization problem is where is the split parameter. Then we can iteratively solve the following subproblems [37, 38]: The minimization of the above subproblems can be iteratively solved by splitting it into the minimizations of and separately. This suggests the following steps.

Step 1. Consider

Step 2. Consider

Step 3. Consider

For the solving of Step 1, we can use an iterative method, for example, the Landweber iteration method, the Levenberg-Marquardt method, and so forth. Since the Levenberg-Marquardt method has a higher convergence rate than the conventional Landweber iteration method, we use the Levenberg-Marquardt method on Step 1. Thus, we will solve a minimization problem as follows.

Step  1. Consider To solve Step  1, we solve the explicitly given variational equation as follows: Step 2 is an norm regularization problem and it can be solved efficiently through the shrinkage operator; that is, where the shrinkage operator

From the three steps, we can see that the speed of the split Bregman method is largely dependent on the speed of solving Step 1, where the computation of the Jacobian matrix is very time-consuming. We list the split Bregman method in Algorithm 2.

Algorithm 2: Reconstruction algorithm based on the split Bregman method.

3.3. Analysis of Standard Reconstruction and Sparsity Reconstruction

The appropriate choice of regularization depends on a priori knowledge of the solutions. Although -norm regularization is a natural choice for its simplicity, it is not always the optimal strategy.

For the DOT reconstruction, the sought-for optical coefficients distribution usually consists of an essentially uninteresting background plus some small inclusions. Thus we will require the solution of the reconstruction in a sparse coefficient vector form; that is, the coefficient vector of the optical parameter contains only a finite number of nonzero elements.

The standard reconstruction and sparsity reconstruction add -norm penalization and -norm penalization into the problem, respectively, which allow additional constraints or prior information towards the approximate solution. Figure 1(a) illustrates how -norm penalization leads to sparse solutions. Take a linear problem as an example; suppose we are looking for a solution of the linear equation in a model space with two degrees of freedom and the line consists of that satisfy . To find the solution with smallest -norm, we can imagine taking a small circle around the origin and increase its radius until it first touches the solution line : the tangent point is the minimum; see Figure 1(b). Obviously, the point of solution has both and components. In a similar way, to find the solution with smallest -norm, we take a small -ball around the origin and increase its radius until it first touches : the touching point is the minimum -norm solution, which is sparser than the solution achieved with smallest -norm, because it is on one axis only; only one component is nonzero.

Figure 1: regularization and regularization.

There is also an extreme penalty term, say, -norm penalization, which is simply the number of nonzero coefficients and leads to the sparsest solution. Nevertheless, the optimization problem with this penalty is not computationally tractable. Hence -norm penalization is preferred in practical problems. In addition, it has been proven that, for some large matrices , if there exist sufficiently sparse solutions, the sparsest solution can be achieved by the -norm minimization [19, 20]. Under certain conditions, the penalty term can provide an accurate result even with limited observations [17]. These phenomena are further investigated in different fields such as electrical impedance tomography and tomographic inversion [39, 40].

4. Numerical Implementations

In this section, we report simulation results on numerical examples in two dimensions (2D). We perform the simulations on a 3.0 GHz PC with 8 GB RAM in MATLAB 2013b environment under Windows 7. In our simulations, the scattering coefficient is reconstructed both with the standard regularization reconstruction and with the sparsity regularization reconstruction. We use two kinds of measurements for reconstruction: boundary angular-averaged measurements and internal angular-averaged measurements. The boundary measurements are the excitance received by detectors attached to the boundary of tissue. Due to the limit of measurement techniques of the optical devices, we cannot accurately receive the excitance from all angles; instead we receive a boundary angular-averaged data which can be regarded as the integration of radiance on the boundary of the domain in all outgoing directions. The internal angular-averaged measurements are the same as boundary angular-averaged measurements except for the location of detectors which we assumed to be inside the domain.

We design different simulations to demonstrate the following points.

First, for sparse coefficient distribution, the sparsity regularization reconstruction localizes the location of the inclusion better than standard regularization reconstruction, especially when there is noise and the data are incomplete. See Figures 4, 9, and 10.

Second, the proposed sparsity regularization reconstruction works better when than . See Figures 6 and 7.

Third, when the internal angular-averaged measurements are available, the proposed sparsity regularization reconstruction works better with internal angular-averaged measurements than with boundary angular-averaged measurements. See Figures 13 and 14.

Last, the split Bregman algorithm can efficiently solve the DOT image reconstruction problem with sparsity regularization. The results reconstructed by the Bregman algorithm are more accurate than those achieved by three other algorithms. See Figures 16 and 17.

We mention here that, in the DOT reconstruction problems, the measurement data are usually synthesized from the numerical solution of the forward problem. Under this situation, the phenomenon of inverse crime will happen especially when the same discretization is used for the forward and inverse process, because it will make the ill-posedness of the inverse problem not evident [41, 42]. Hence, in order to avoid the inverse crime, we will use different discretization meshes in the forward and inverse problems.

Now let us state our numerical experiments. Before stating the details of every experiment, we first state some common experiment settings that will be used in all of our experiments. Our purpose is to reconstruct the scattering coefficient based on the BVP (2) and (3), where we assume that the following 2D simulations are all performed on a unit circular domain centered at (0 mm, 0 mm), corresponding to the the internal light source , the inflow current for every incident impulse that is settled as The boundary value is a piecewise linear function whose spatial support is and achieves the value at the center node of ; here denotes the finite element through which the th incident impulse passes. The direction of the th incident impulse is node of the angular mesh that points approximately from the center of to the center of .

If we denote the domains that contain the inclusions as , then the background domain is . Then the true value of absorption coefficient in BVP (2) and (3) is defined as Since the absorption coefficient is assumed to be known in our scattering reconstruction problems, we only give the mathematical formulation of the absorption coefficient. While in order to compare explicitly the true value and the reconstructed value of scattering coefficient, we demonstrate the true value and the reconstructed results of scattering coefficient in figures.

In order to make discretization in angle, we can divide the angular space uniformly into directions with equal interval length as shown in Figure 2. And we set in all examples; that is, there are 32 angular directions.

Figure 2: Angular discretization in 2D.

In the first example, Figure 3, a small inclusion with 0.1 mm diameter, is centered at . 12 sources and 12 detectors are equally spaced attached to the boundary of the circle domain. The simulated boundary angular-averaged measurements performing on the boundary are generated with a spatial mesh (a) of 1097 nodes and 2104 elements and reconstruction mesh (b) having 286 nodes and 526 elements.

Figure 3: Standard regularization reconstruction for a small spot with 132 boundary angular-averaged measurements. (a) Mesh for the forward problem, (b) mesh for the inverse problem, (c) the true scattering distribution, and (d) the image reconstructed from standard regularization.
Figure 4: Sparsity regularization reconstruction for a small spot with 132 boundary angular-averaged measurements, images with (a) 0%, (b) 0.1%, (c) 1%, and (d) 10% Gaussian noise, respectively.

The image from standard regularization reconstruction without noise is displayed in (d) in Figure 3. The images from sparsity regularization reconstruction without or with noise are displayed in Figure 4 with the noise defined as , where is the signal-to-noise ratio and is a Gaussian random variable with zero mean and unity variation. In Figure 4, we vary from 0 percent to 10 percent and plot the sparsity regularization reconstruction image. Comparing images from standard regularization and sparsity regularization, we conclude that sparsity regularization can efficiently localize the sparse inclusion with a few measurements, even in the presence of different degrees of noise, while standard regularization fails to reconstruct the sparse inclusion even without noise.

In the second example, we compare the sparsity regularization reconstruction and the standard regularization reconstruction on squares of three different sizes.

In Figure 5, (a), (c), and (e) show the forward mesh of three different sizes squares. (b), (d), and (f) show the true distribution of the scattering coefficient in three different experiments.

Figure 5: Meshes for square cases. (a), (c), and (e) are forward meshes for three square cases, respectively. (b), (d), and (f) are the true scattering coefficient distribution for big square case, middle square case, and small square case, respectively.
Figure 6: Comparison of sparsity regularization reconstruction and standard regularization reconstruction with 132 boundary angular-averaged measurements when . (a), (c), and (e) are reconstructed images with sparsity regularization for big square, middle square, and small square, respectively. (b), (d), and (f) are reconstructed images with standard regularization for big square, middle square, and small square, respectively.
Figure 7: The same reconstructions as in Figure 6 but with . (a), (c), and (e) are reconstructed images with sparsity regularization for big square, middle square, and small square, respectively. (b), (d), and (f) are reconstructed images with standard regularization algorithm for big square, middle square, and small square, respectively.

Comparing (a), (c), and (e) in Figures 6 and 7 and (b), (d), and (f) in Figures 6 and 7, clearly, the sparsity regularization can localize the position of the inclusion better than the standard regularization in all the three kinds of squares and has a clear contrast with the backgrounds. But the standard regularization also has its advantages. Comparing (a) and (b), (c) and (d), and (e) and (f) in Figure 6, when , there are many blurred dots in the standard regularization reconstructed images; we can not identify the distribution domain of the reconstructed scattering coefficient, but the value of the scattering coefficient in standard regularization reconstruction is much closer to the true value of the scattering coefficient than that in the sparsity regularization reconstruction. However, this advantage is broken when ; see Figure 7.

On the other hand, the proposed sparsity regularization reconstruction performs better when than when in localizing the position of the inclusions as well as identifying the value of the scattering coefficient.

In conclusion, when the anisotropic factor is small, the standard regularization can identify the value of the scattering coefficient better than the sparsity regularization. Although there are some blurred dots in standard regularization reconstructed images, one can remove them by using multilevel approach. On the other hand, the sparsity regularization can localize the position of the inclusion better than standard regularization and has a clear contrast, especially in the forward-peaking regime with big anisotropic factor .

In the third example (Figure 8), three-letter inclusions (a) are put in the domain to evaluate our reconstruction algorithms; the forward mesh (a) has 1077 nodes and 2072 elements; the reconstruction mesh (b) has 280 nodes and 518 elements.

Figure 8: Meshes for three-letter case. (a) is the forward mesh. (b) is the inverse mesh. (c) is the true scattering coefficient distribution.
Figure 9: Reconstruction with boundary data when . (a) and (c) are reconstructed images by sparsity regularization with 132 boundary angular-averaged measurements and 380 boundary angular-averaged measurements. (b) and (d) are reconstructed images by standard regularization with 132 boundary angular-averaged measurements and 380 boundary angular-averaged measurements.
Figure 10: The same reconstructions as in Figure 9 but with . (a) and (c) are reconstructed images by sparsity regularization with 132 boundary angular-averaged measurements and 380 boundary angular-averaged measurements. (b) and (d) are reconstructed images by standard regularization with 132 boundary angular-averaged measurements and 380 boundary angular-averaged measurements.

In Figures 9 and 10, the reconstructions with and 0.9, respectively, are carried out to demonstrate that sparsity regularization ((a), (c) in Figures 9 and 10) in general reconstructs better images than standard regularization ((b), (d) in Figures 9 and 10).

Comparing (a) and (c) in Figures 9 and 10, respectively, we find that the proposed sparsity regularization reconstruction performed better when the inverse problem is more severely ill-posed due to fewer measurements.

In conclusion, efficient algorithm and proper regularization, for example, sparsity regularization reconstruction for sparse coefficient distribution, are essential to recover high-resolution images. On the other hand, we can notice that there are blurred dots in the middle of the reconstructed images in (a) and (c) in Figures 9 and 10, due to the ill-posedness of inverse problem. These blurred dots have small contrast and thus can be removed through sparsity regularization reconstruction by choosing proper and or, beyond our paper, through multilevel approach [43].

Next, using sparsity regularization reconstruction, we will show the superiority of the internal angular-averaged measurements over the boundary angular-averaged measurements with two relatively complicated cases. The detectors for internal angular-average measurements in two cases lie on a 0.95-radius circle equidistant from the center at (0 mm, 0 mm). Case 1 is shown in the first column in Figures 13 and 14; Case 2 is shown in the second column in Figures 13 and 14.

Figure 11 shows the meshes for Case 1. (a) shows the forward mesh; it has 1131 nodes and 2160 elements. (b) shows the inverse mesh; it has 296 nodes and 540 elements. (c) shows the true distribution of the scattering coefficient. Figure 12 shows the meshes for Case 2. (a) shows the forward mesh; it has 1347 nodes and 2592 elements. (b) shows the inverse mesh; it has 350 nodes and 648 elements.

Figure 11: Meshes for Case 1. (a) is the forward mesh. (b) is the inverse mesh. (c) is the true scattering coefficient distribution.
Figure 12: Meshes for Case 2. (a) is the forward mesh. (b) is the inverse mesh. (c) is the true scattering coefficient distribution.
Figure 13: Sparsity regularization reconstruction when . (a) and (c) are reconstructed images for Case 1 with 132 boundary angular-averaged measurements and 132 internal angular-averaged measurements, respectively. (b) and (d) are reconstructed images for Case 2 with 240 boundary angular-averaged measurements and 240 internal angular-averaged measurements, respectively.
Figure 14: The same reconstructions as in Figure 13 but with . (a) and (c) are reconstructed images for Case 1 with 132 boundary angular-averaged measurements and 132 internal angular-averaged measurements, respectively. (b) and (d) are reconstructed images for Case 2 with 240 boundary angular-averaged measurements and 240 internal angular-averaged measurements, respectively.

Comparing (a) and (c) in Figures 13 and 14, we find that, for Case 1, the internal angular-averaged measurements ((c) in Figures 13 and 14) can localize the position of the outside inclusion, although the internal inclusion is localized blurred while the boundary angular-averaged measurements ((a) in Figures 13 and 14) can not completely localize any inclusion.

For the complicated Case 2 ((b) and (d) in Figures 13 and 14), we increase the measurement data by increasing the number of detectors to 16, respectively, in order to get more information from the boundary measurements. Both of the two kinds of measurements can not accurately identify the inclusion, but the internal angular-averaged measurements ((d) in Figures 13 and 14) performed relatively better than the boundary angular-averaged measurements ((b) in Figures 13 and 14). A justification of this phenomenon is that the energy of the incident current is decayed much due to the large and relatively frequently change of the scattering coefficient of the inclusion. One can alleviate this phenomenon by increasing the number of detectors or measurement data. In conclusion, the internal data can better-pose the inverse problem.

As the last simulation, we compute Case 3 in Figure 15 to investigate the performance of split Bregman algorithm for DOT image reconstruction with sparsity regularization. The results are compared with Gauss-Newton algorithm [44] and the state-of-the-art algorithms including GPSR [45] and YALL1 [46]. To deal with the nondifferentiability of the absolute value at in Gauss-Newton algorithm, we replace by ,  . We choose in Gauss-Newton algorithm.

Figure 15: Meshes for Case 3. (a) is the forward mesh. (b) is the inverse mesh. (c) is the true scattering coefficient distribution.
Figure 16: Sparsity regularization reconstruction when with the measurement noise level of 0.1% and 132 boundary angular-averaged measurements, using four different algorithms. (a)–(d) are reconstruction results by the split Bregman, YALL1, GPSR, and Gauss-Newton algorithm, respectively.
Figure 17: Sparsity regularization reconstruction when with the measurement noise level of 0.1% and 132 boundary angular-averaged measurements, using four different algorithms. (a)–(d) are reconstruction results by the split Bregman, YALL1, GPSR, and Gauss-Newton algorithm, respectively.

In Case 3, two circle inclusions are embedded at the top and bottom of the circle domain with 1137 nodes and 2168 elements for the forward mesh and 298 nodes and 542 elements for the inverse mesh. We use 132 boundary angular-averaged measurements with 0.1% Gaussian noise for reconstruction. The reconstructed images are shown in Figures 16 and 17, and we summarize the results with in Table 1, containing the parameters for the four methods, computational time, data misfit , relative solution error norm (RE), and signal-to-noise ratio (SNR). The relatively optimal parameters are chosen empirically.

Table 1: Comparison of sparsity reconstructions for Case 3 when with the regularization parameter , is the split parameter of split Bregman method, represents the tolerance for stopping criterion of YALL1 and GPSR method, and is the control parameter of Gauss-Newton method.

The RE is calculated as and SNR is calculated as

The reconstruction results in Figures 16 and 17 show that all of the four algorithms can locate the position of two inclusions. The split Bregman and GPSR algorithm can achieve more accurate shape of the inclusions, and the reconstructed scattering coefficients achieved by the two algorithms are closer to the true value. The Gauss-Newton algorithm can find the approximate location of the inclusions, but the reconstructed circle inclusion is smaller in size and the reconstructed scattering coefficient departs from actual ones.

The results displayed in Table 1 also demonstrate that the split Bregman algorithm performs better. It leads to the lowest RE and data misfit with less calculation. The GPSR algorithm achieves higher SNR at the expense of about 25% extra computational time. The Gauss-Newton algorithm spends least computational time, but the reconstruction error is not very satisfactory.

These results are justified by the fact that the split Bregman algorithm decouple the sparsity reconstruction problem into and portions leading to a better compromise in the efficiency and quality of the reconstructed optical parameters.

5. Conclusion

In this paper, by using the image modalities in DOT, we employ the sparsity regularization method on the RTE-based coefficient identification problems, which is proven to perform in general better than the standard regularization reconstruction, especially for sparse distribution coefficient and large noise, and in the forward-peaking regime with big anisotropic factor . On the other hand, we construct cases and compare the reconstruction results with boundary and internal measurement to test the validity of the proposed method; results show that the proposed method is practicable and feasible; it performs steadily with various measurement data; meanwhile, the internal measurement can better-pose the inverse problem and achieve more accurate results. However, we usually cannot obtain the internal measurement data in practice. Hence, our method cannot reconstruct inclusions with complicated internal structure accurately with small amount of boundary angular-averaged measurements. One can alleviate this phenomenon by increasing the number of detectors or measurement data. In the further work we will consider seeking for multi-imaging modality which can further improve the inversion quality with boundary angular-averaged measurements. Results show the competitive performance of the split Bregman algorithm for the DOT image reconstruction compared with other existing algorithms.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to express their gratitude to the referees for their careful reading and detailed suggestions, which helped to improve the revised version of this paper. This work is supported by the National Nature Science Foundation of China under Grant no. 91230119, no. 81071207, and no. 81271622, and by “the Fundamental Research Funds for the Central Universities” (HIT.IBRSEM.B.201416).

References

  1. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems, vol. 15, no. 2, pp. R41–R93, 1999. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  2. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Problems, vol. 25, no. 12, Article ID 123010, 2009. View at Publisher · View at Google Scholar · View at Scopus
  3. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Physics in Medicine and Biology, vol. 50, no. 4, article R1, 2005. View at Publisher · View at Google Scholar · View at Scopus
  4. W. Han, J. A. Eichholz, X. L. Cheng, and G. Wang, “A theoretical framework of x-ray dark-field tomography,” SIAM Journal on Applied Mathematics, vol. 71, no. 5, pp. 1557–1577, 2011. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  5. W. Han, J. A. Eichholz, and G. Wang, “On a family of differential approximations of the radiative transfer equation,” Journal of Mathematical Chemistry, vol. 50, no. 4, pp. 689–702, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  6. H. Gao and H. K. Zhao, “A fast-forward solver of radiative transfer equation,” Transport Theory and Statistical Physics, vol. 38, no. 3, pp. 149–192, 2009. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  7. K. Ren, G. Bal, and A. H. Hielscher, “Transport- and diffusion-based optical tomography in small domains: a comparative study,” Applied Optics, vol. 46, no. 27, pp. 6669–6679, 2007. View at Publisher · View at Google Scholar · View at Scopus
  8. J. Nocedal and S. J. Wright, Numerical Optimization, Springer Series in Operations Research, Springer, New York, NY, USA, 1999. View at Publisher · View at Google Scholar · View at MathSciNet
  9. M. Hanke, “The regularizing Levenberg-Marquardt scheme is of optimal order,” Journal of Integral Equations and Applications, vol. 22, no. 2, pp. 259–283, 2010. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  10. T. Feng, P. Edstrom, and M. Gulliksson, “Levenberg-Marquardt methods for parameter estimation problems in the radiative transfer equation,” Inverse Problems, vol. 23, no. 3, pp. 879–891, 2007. View at Publisher · View at Google Scholar
  11. H. Gao, S. Osher, and H. K. Zhao, “Quantitative photoacoustic tomography,” in Mathematical Modeling in Biomedical Imaging II, Lecture Notes in Mathematics, pp. 131–158, Springer, Berlin, Germany, 2012. View at Publisher · View at Google Scholar
  12. I. Daubechies, M. Defrise, and C. de Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, no. 11, pp. 1413–1457, 2004. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  13. R. Ramlau, “Regularization properties of Tikhonov regularization with sparsity constraints,” Electronic Transactions on Numerical Analysis, vol. 30, pp. 54–74, 2008. View at Google Scholar · View at MathSciNet · View at Scopus
  14. S. Osher, M. Burger, D. Goldfarb, J. Xu, and W. Yin, “An iterative regularization method for total variation-based image restoration,” Multiscale Modeling & Simulation, vol. 4, no. 2, pp. 460–489, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  15. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena, vol. 60, no. 1–4, pp. 259–268, 1992. View at Publisher · View at Google Scholar · View at Scopus
  16. Y. Tsaig and D. L. Donoho, “Extensions of compressed sensing,” Signal Processing, vol. 86, no. 3, pp. 549–571, 2006. View at Google Scholar
  17. D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  18. B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM Journal on Computing, vol. 24, no. 2, pp. 227–234, 1995. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  19. D. L. Donoho, “For most large underdetermined systems of linear equations the minimal l1-norm solution is also the sparsest solution,” Communications on Pure and Applied Mathematics, vol. 59, no. 6, pp. 797–829, 2006. View at Google Scholar
  20. E. J. Candes, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Communications on Pure and Applied Mathematics, vol. 59, no. 8, pp. 1207–1223, 2006. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  21. J. Chamorro-Servent, J. F. P. J. Abascal, J. Aguirre et al., “Use of split bregman denoising for iterative reconstruction in fluorescence diffuse optical tomography,” Journal of Biomedical Optics, vol. 18, no. 7, Article ID 076016, 2013. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus
  22. J.-F. Cai, S. Osher, and Z. W. Shen, “Linearized Bregman iterations for compressed sensing,” Mathematics of Computation, vol. 78, no. 267, pp. 1515–1536, 2009. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  23. T. Goldstein and S. Osher, “The split Bregman method for L1-regularized problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323–343, 2009. View at Publisher · View at Google Scholar · View at MathSciNet
  24. W. T. Yin, S. Osher, J. Darbon, and D. Goldfarb, “Bregman iterative algorithms for l1-minimization with applications to compressed sensing,” SIAM Journal on Imaging Science, vol. 1, no. 1, pp. 143–168, 2008. View at Google Scholar
  25. W. T. Yin, “Analysis and generalizations of the linearized Bregman model,” SIAM Journal on Imaging Sciences, vol. 3, no. 4, pp. 856–877, 2010. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  26. T. Tarvainen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio, “Image reconstruction in diffuse optical tomography using the coupled radiative transport-diffusion model,” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 112, no. 16, pp. 2600–2608, 2011. View at Publisher · View at Google Scholar · View at Scopus
  27. L. C. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” The Astrophysical Journal, vol. 93, pp. 70–83, 1941. View at Publisher · View at Google Scholar
  28. H. B. Jiang, Diffuse Optical Tomography: Principles and Applications, CRC Press, Boca Raton, Fla, USA, 1st edition, 2010.
  29. V. Agoshkov, Boundary Value Problems for Transport Equations, Modeling and Simulation in Science, Engineering and Technology, Birkhäuser, Boston, Mass, USA, 1998. View at Publisher · View at Google Scholar · View at MathSciNet
  30. J. Tang, W. Han, and B. Han, “A theoretical study for RTE-based parameter identification problems,” Inverse Problems, vol. 29, no. 9, Article ID 095002, 2013. View at Publisher · View at Google Scholar
  31. H. Gao and H. K. Zhao, “Analysis of a numerical solver for radiative transport equation,” Mathematics of Computation, vol. 82, no. 281, pp. 153–172, 2013. View at Publisher · View at Google Scholar · View at PubMed · View at MathSciNet · View at Scopus
  32. W. Han, J. A. Eichholz, J. Huang, and J. Lu, “RTE-based bioluminescence tomography: a theoretical study,” Inverse Problems in Science and Engineering, vol. 19, no. 4, pp. 435–459, 2011. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  33. W. Han, J. Huang, and J. A. Eichholz, “Discrete-ordinate discontinuous Galerkin methods for solving the radiative transfer equation,” SIAM Journal on Scientific Computing, vol. 32, no. 2, pp. 477–497, 2010. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  34. E. E. Lewis and W. F. Miller, Computational Methods for Neutron Transport, Wiley-Interscience, 1984.
  35. M. Gehre, T. Kluth, A. Lipponen et al., “Sparsity reconstruction in electrical impedance tomography: an experimental evaluation,” Journal of Computational and Applied Mathematics, vol. 236, no. 8, pp. 2126–2136, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  36. J. Kaipio and E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, NY, USA, 2005. View at MathSciNet
  37. J. Bush, Bregman algorithms [M.S. thesis], University of California, Santa Barbara, Calif, USA, 2011.
  38. J. Wang, J. Ma, B. Han, and Q. Li, “Split Bregman iterative algorithm for sparse reconstruction of electrical impedance tomography,” Signal Processing, vol. 92, no. 12, pp. 2952–2961, 2012. View at Publisher · View at Google Scholar · View at Scopus
  39. B. Jin, T. Khan, and P. Maass, “A reconstruction algorithm for electrical impedance tomography based on sparsity regularization,” International Journal for Numerical Methods in Engineering, vol. 89, no. 3, pp. 337–353, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  40. I. Loris, G. Nolet, I. Daubechies, and F. A. Dahlen, “Tomographic inversion using 1-norm regularization of wavelet coefficients,” Geophysical Journal International, vol. 170, no. 1, pp. 359–370, 2007. View at Google Scholar
  41. D. L. Colton and R. Kress, Inverse Acoustic and Elctromagnetic Scattering Theory, Springer, Berlin, Germany, 2nd edition, 1998.
  42. M. Jiang, T. Zhou, J. Cheng, W. X. Cong, and G. Wang, “Image reconstruction for bioluminescence tomography from partial measurement,” Optics Express, vol. 15, no. 18, pp. 11095–11116, 2007. View at Publisher · View at Google Scholar · View at Scopus
  43. H. Gao and H. K. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation. Part 2: total variation and l1 data fidelity,” Optics Express, vol. 18, no. 3, pp. 2894–2912, 2010. View at Publisher · View at Google Scholar · View at PubMed
  44. M. Schweiger, S. R. Arridge, and I. Nissilä, “Gauss-Newton method for image reconstruction in diffuse optical tomography,” Physics in Medicine and Biology, vol. 50, no. 10, pp. 2365–2386, 2005. View at Publisher · View at Google Scholar · View at PubMed · View at Scopus
  45. M. A. T. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: application to compressed sensing and other inverse problems,” IEEE Journal on Selected Topics in Signal Processing, vol. 1, no. 4, pp. 586–597, 2007. View at Publisher · View at Google Scholar · View at Scopus
  46. J. Yang and Y. Zhang, “Alternating direction algorithms for l1 problems in compressive sensing,” SIAM Journal on Scientific Computing, vol. 33, no. 1, pp. 250–278, 2011. View at Google Scholar