Abstract

Optical tomography is an emerging and important molecular imaging modality. The aim of optical tomography is to reconstruct optical properties of human tissues. In this paper, we focus on reconstructing the absorption coefficient based on the radiative transfer equation (RTE). It is an ill-posed parameter identification problem. Regularization methods have been broadly applied to reconstruct the optical coefficients, such as the total variation (TV) regularization and the regularization. In order to better reconstruct the piecewise constant and sparse coefficient distributions, TV and norms are combined as the regularization. The forward problem is discretized with the discontinuous Galerkin method on the spatial space and the finite element method on the angular space. The minimization problem is solved by a Jacobian-based Levenberg-Marquardt type method which is equipped with a split Bregman algorithms for the regularization. We use the adjoint method to compute the Jacobian matrix which dramatically improves the computation efficiency. By comparing with the other imaging reconstruction methods based on TV and regularizations, the simulation results show the validity and efficiency of the proposed method.

1. Introduction

Optical tomography with near-infrared light is a promising technique for noninvasive studying of the functional characters of human tissues. One can find many applications of this technique, for example, the early detection of breast cancer, cervical cancer screening, monitoring of infant brain tissue oxygenation level and functional brain activation, and the study of dosimetry in photon dynamic therapy; see [14]. Unlike the X-ray tomography, optical tomography uses the near-infrared (NIR) light as the probing radiation. Different human tissues have different absorption and scattering properties which affect the transmission of the NIR light. Therefore, we can identify the location and quantity of abnormal tissues by the emerging light measured on the boundary.

Optical tomography is a high contrast imaging modality. However, currently only low resolution reconstructions are possible, especially when using an unmodulated continuous wave source [5]. Another barrier of optical tomography is that reconstructed images have a poor quality particularly when abnormal targets are located deep in tissues. On the other hand, due to the limited measurement data, as well as the high scattering and high absorption properties, the identification problem of optical tomography usually is ill-posed and underdetermined. The ill-posedness makes the coefficient reconstruction sensitive to small perturbation from measurements, such as noise and computational error. Therefore, various reconstruction algorithms based on regularization are developed to obtain reasonable and stable reconstructions [69]. Recently, the regularization method with norm has been used in optical tomography [10, 11]. In [10], Levenberg-Marquardt strategy is applied for solving the regularization step of split Bregman algorithm. The contrast tests show the superiority of regularization over regularization. The numerical experiments therein also show that the regularization has a better utility when the independent measurements are much more limited. In [11], linearized Bregman iteration based on the Bregman distance is exploited to minimize the sparse regularization. The experimental results in numerical simulation of an in vivo mouse demonstrate the effectiveness of the algorithm. However, the sparsity regularization may oversparsify the distribution of the coefficients.

Among many types of regularizations, TV regularization is often used due to its abilities of well reconstructing the discontinuous and piecewise constant distributions. TV regularization was first introduced in the field of image processing and image reconstruction; see [1214] and the references therein. To date, TV regularization has been the preferred regularization strategy of considerable recent works, like photoacoustic tomography (PAT) [1517], bioluminescence tomography (BLT) [18, 19], fluorescence tomography (FT) [20], optical coherence tomography (OCT) [21], and so on. In [18], the split Bregman algorithm for TV regularization (SBRTV) is applied to the reconstruction of source distribution in BLT. The algorithm is evaluated with 2D and 3D simulations and 3D in vivo experiments. In the reported results, the source distribution can be reconstructed with better accuracy for the location with SBRTV regularization than with or regularizations. In [19], the synergistic combination of Bregman method and TV regularization is utilized to quantitatively improve the reconstruction of absorption and scattering coefficients for both Jacobian-based and gradient-based methods in quantitative PAT. The feasibility of the algorithm is checked with 3D simulations.

For optical tomography, the coefficients are usually taken as piecewise constant. Since TV regularization is effective for piecewise constant reconstruction, it is a natural choice for reconstructing coefficients of piecewise constant distributions in optical tomography. However, due to the nonsmoothness and nondifferentiability, TV regularization is difficult to realize computationally. TV regularization for optical tomography based on the radiative transport equation has been studied in [22], where inexact Gauss-Newton method is used to solve the TV regularization problem. In [23], norm is chosen as the penalty in the regularization strategy for optical tomography based on the frequency radiative transport equation.

In this paper, we use regularization as a combination of TV and norm. To reduce the computational complexity, we apply the split Bregman method to solve the joint regularization problem. It is reasonable to combine regularization and TV regularization to improve the reconstruction quality in optical tomography. Moreover, based on the split Bregman method, the iterative algorithms can be computationally efficient. Various experiments in 2D are performed to evaluate the performance of the algorithm.

The rest of this paper is organized as follows. In Section 2, we briefly describe the forward and inverse problems for optical tomography. In Sections 3 and 4, we show the implementation details about computing the Fréchet derivative of the forward operator and the iterative procedure. In Section 5, numerical results are presented.

2. Optical Tomography

We consider two mathematical problems: the forward problem and the inverse problem. For the forward problem, based on the physical model of light propagation in tissues, for a given set of optical properties, we model the measurements on the boundary. For the inverse problem, the optical properties can be reconstructed by matching the predictions calculated from the forward problem and the measurements from the detectors.

2.1. Forward Problem

The light propagation in tissues is described by the radiative transport equationHere, , or , denotes a bounded convex domain with a boundary and denotes the unit sphere of . The variables and denote the spatial position and the angular direction. describes the density of photons. The expression denotes the directional derivative at position along the direction . The nonnegative normalized phase function is the probability that photons traveling in the direction are scattered into the direction . In optical tomography, the phase function usually is taken as the Henyey-Greenstein phase function (cf. [24]): in two dimensions, it is of the form where the parameter is the anisotropy factor of the scattering medium. The absorption coefficient is denoted by , and the scattering coefficient is denoted by . For the optical parameters and , the following conditions are assumed to hold throughout this presentation.

Assumptions(A1)The function is uniformly positive and bounded; that is, there exist two positive constants and such that a.e. in .(A2)The function is uniformly positive and bounded; that is, there exist two positive constants and such that a.e. in .

Equation (1) is supplemented by boundary conditions. Similar to the X-ray CT, optical tomography experiments acquire the current distribution of detectors on the boundary under multi-incidents. Let , , be disjoint, connected subsets of . Corresponding to incident sources on , define bywhere is the unit outward normal vector at .

Corresponding to each independent incident , the measurable quantity is the outgoing light on the boundary of domain, which can be written asIf we assume there are detectors located at different positions , then the optical tomography experiment consists of exciting the domain with a sequence of incident source and recording the corresponding measurements data . Here, and represent the row index and column index in matrix , respectively.

With above notations, a mathematical description of such an experiment is provided by the following nonlinear forward operator: which maps prescribed optical parameters to the corresponding measurements data. Here, denotes the th forward operator corresponding to the th incident source and the resulting measurement data on detectors. The forward operator is well defined for and in the set(cf. [25]).

There are many references on the discretization of RTE; see, for instance, [2630]. In this paper, we use continuous linear elements and discontinuous Galerkin method with piecewise linear functions to discretize the angular variables and the spatial variables, respectively. An upwind numerical flux is used to approximate the incoming flux through the surface of the control element and inflow boundary. After assembling the full discretization formulation and forming a system of tebra equation, we solve the linear system by Gauss-Seidel method. We use piecewise constants to approximate the absorption and scattering coefficients. We can express aswhere denotes the number of the elements, denotes the character function corresponding to the th element, and and denote the values of absorption and scattering coefficients on the th element.

2.2. Inverse Problem

The inverse problem of optical tomography is to determine the unknown coefficients and from the boundary detector readings. In this paper, we only reconstruct the absorption coefficient assuming that the scattering coefficient is known; then the forward operator in fact acts on the unknown only. Thus, the inverse problem of optical tomography is to determine from the following system of nonlinear equations:Then the optical tomography can be formulated as minimizing the difference between the measurement data and the model predictionsknown as data fidelity. The inverse problem (9) is ill-posed in the sense that the amount of the measurements is quite limited compared with the number of the unknowns and that the measurements contain noises. To overcome the ill-posedness, the data fidelity should be combined with appropriate regularization. In the regularization strategy, we minimize the following objective functional:Here, is the regularization penalty functional that enforces a prior information on and is the regularization parameter that trades off the weight between the discrepancy term and the penalty functional.

According to different ways of treating the derivative of data fidelity, there are two sorts of approaches for solving (10). One is to linearize the forward operator near the th iteration of , which is denoted as , as follows:where is the derivative of the forward operator with respect to . The linearized formulation provides a good approximation when is close to the true value. As a result, the minimized problem (10) is treated by an iterative procedure as follows:This Jacobian-based minimization problem can be solved with many iteratively optimization techniques, such as the Levenberg-Marquardt type method (LM). LM type method is the special case of Gauss-Newton type method. The standard LM method (or Gauss-Newton method) uses as the penalty. The LM type method defines the update in a region around , while the Gauss-Newton method always defines in a neighbourhood of the initial guess for each . Hence, from the optimization point of view, LM type method is more favourable in nature. Based on the Jacobian-based method, many optimization techniques such as split Bregman method can be used to solve (12). On the other hand, when is close to the true value, the Jacobian-based method shares the typical quadratic convergence from Newton method.

Alternatively, (10) can be solved directly by minimizing the nonlinear functional with some gradient-based methods, like Quasi-Newton method, nonlinear conjugate gradient method, limited-memory BFGS method, and so on. Thus, the computation of the linearized Jacobian matrix is avoided and only the gradient of the nonlinear minimizing functional needs to be computed. The gradient-based method in general has superlinear convergence and is economic in memory storage which is suitable for large scale problems, such as 3D problems. Since our numerical experiments are all done in two dimensions and the problem scale is not so large as that in three dimensions, we use the Jacobian-based LM type method to solve (10).

3. Adjoint Problem

Since we adopt the Jacobian-based minimization method to solve (10) in this paper, the derivative of the forward operator becomes a matrix which is called Jacobi matrix. For convenience of representing the element of Jacobi matrix hereinafter, we use to represent the Jacobi matrix of , with respect to , ; here denotes the th element of the discretized . For each , is defined as follows: can be computed by direct method which differentiates with respect to the perturbation of on each element. In order to compute , for each , we need to solve forward problems times. For a total of sources, we need to compute times, so the direct method is very time consuming. Therefore, we use the adjoint formulation to compute instead of direct method.

We first consider the analytic Fréchet derivative of the forward operator for the boundary value problem of RTE given by [25].where denotes the adjoint operator of and is the solution of the adjoint RTEwith boundary conditionFrom the adjoint RTE (15) and the boundary condition (16), it seems that the adjoint problem should be solved in the reverse direction of propagation with a completely new computation solver. However, by the standard reciprocity theorem for the Boltzmann equation given in [31],which states that the angular density at in direction due to a source at in direction is the same as the angular density at in direction due to a source at in direction . Here, presents the Green function of (1) for isotropic point source on the boundary. Then the adjoint problems (15) and (16) can be transformed into the same form as the forward problem, simply replacing direction with . Then we just need to solve the following equation for the radiance by the same forward solver:with adjoint boundary conditionThis means we need to solve (18) and (19) with the forward solver firstly and then reverse all directions on solution.

If we consider one source position and one detector position , then the th column of the Jacobi matrix can be computed as follow:where is the solution of the following boundary value problem:Thus for the adjoint method to compute , for each , we only need to solve forward problems times. For sources, we need to solve forward problems times, which dramatically reduces the computation burden and improves the efficiency of the algorithm.

4. Iterative Procedure

The unknown can be estimated through the regularization method. The quality of reconstructed image strongly depends on the choice of the penalty term . If we choose the norm as penalty, that is, , the reconstructed image usually blurs with a low resolution. If we choose norm as the penalty, that is, , the reconstructed image tends to find a sparse solution [10, 18]. To treat the discontinuity and the edges of different distribution regions, total variation is usually chosen as the penalty functional, that is,The symbol denotes the total variation seminorm [32] of as follows:To improve the reconstruction quality, we consider the total variation mixed with the norm as the penalty. The functional to be minimized is of the formwhere are the regularization parameters. We will apply the split Bregman method to solve (24) [33]. Instead of (24), we consider the following constrained optimization problem:Solution of the above minimization problem can be obtained by solving the unconstrained optimization problemwhere is the split parameter. Now let us iteratively solve the following subproblems [34]:with the following update for : The minimization of subproblems in (1) can be iteratively solved by splitting it into the minimizations of and separately. This suggests the following steps.

Step 1.

Step 2.

Step 3.

For the solution of Step 1, we use the Levenberg-Marquardt method, which has a high convergence rate. Thus we solve a minimization problem as follows.

Step where denotes the Jacobi matrix of the forward operator with respect to . For solving Step  , in order to avoid the nondifferentiability of the total variation term at a zero point, we approximate it with a smooth functional defined asThe parameter is a constant, and it cannot be too large or too small. The Euler-Lagrange equation for Step iswhere denotes the transpose, is the short notion of , the term is defined asThen we can get the iterative update by solving the above equation with Newton method as follows:Step 2 is an norm regularization problem and it can be solved efficiently through the shrinkage operator, that is,where the shrinkage operatorImplementation of the split Bregman formulation for our - regularization is described in Algorithm 1.

Input: Given the initial guess ; the regularization parameter , , ;
  , and the tolerance .
while do
   For , compute with (36);
   Compute ;
   Compute .
end while
Output: Output an approximation .

5. Reconstruction Results and Discussion

In this section, the - regularization with the split Bregman formulation is applied to 2D test problems. We assume the distribution of the scattering coefficient is known. Reconstructions of spatially dependent distributions of absorption coefficient inside the medium are performed and discussed. In our simulation, a circular domain which contains different inclusions is investigated. The radius of the circle is 10 mm. sources and detectors are located on the boundary of the domain with equal space. This yields totally source-detector pairs to be used in the inversion.

Noise-free synthetic data are generated by solving the forward problem on triangular meshes with the method we mentioned in Section 2.1. Note that the meshes for the inverse problem are coarser than the meshes for the forward problem in order to avoid the “inverse crime” [35] and for the purpose of testing the robustness of the proposed algorithm. We will describe the meshes in each case. In all the cases, the angular space is discretized into directions which equally divide the interval .

The background medium has scattering coefficient of 10 m and absorption coefficient of 0.01 m, respectively, and these values keep the same throughout this paper. The initial guess of reconstruction is set to be identical to the properties of the background medium. That is, the iterative procedure started with the background value of the absorption coefficient.

The incident impulse on the inflow boundary is settled aswhere is a piecewise linear function whose spatial support is and achieves the value 1 at the center node of , where is the element through which the ith incident impulse passes. The direction of points approximately from the center of to the center of .

Based on above settings, we use various simulations to validate the proposed split Bregman algorithm for our - regularization. Our purpose is to show the following results.

First, by comparing the reconstruction with different anisotropic factors , the proposed algorithm works better when is bigger.

Second, for small sparse inclusions, the - regularization can reconstruct the location and the quality of the coefficient more accurately than the regularization and the regularization.

Third, the proposed split algorithm for - needs less computation time than the Levenberg-Marquardt algorithm for regularization. Even though our proposed algorithm needs more computation time than the split Bregman algorithm for regularization [10], the reconstruction quality with the proposed algorithm is much better than that with the split Bregman algorithm for regularization.

5.1. Simulation 1: Reconstruction with Different Anisotropic Factors

In this simulation, we reconstruct one circle inclusion with the radius 0.5 mm centered at (7.0 mm, 0.0 mm). In this circle, the absorption coefficient and . Our goal is to reconstruct . For solving the inverse problem, we use a mesh of 772 nodes and 1484 triangular elements. To avoid the inverse crime, we use a mesh of 1296 nodes and 2440 triangular elements for the synthetic data. The simulated absorption distributions and the meshes can be seen in Figure 1. We solve the minimizer of by using the split Bregman method for - method with exact data, that is, , in this example. Then we compare the reconstruction results for different anisotropic factors . For the comparison purpose, we use the same parameters α, β, η, and ε for different . Noticing that the so-called exact data in fact contains noise, so the regularization parameter should not be too small. Since the inclusion is very small in this example, we enhance more weight of penalty than the weight of TV penalty. Here we take and . The other two parameters are taken as and , respectively.

The reconstruction results can be seen in Figure 2. The first row of Figure 2 is the reconstruction results for and from the left to the right. The second row of Figure 2 is the reconstruction results for and from the left to the right. As it can be seen, the reconstructions for the bigger are more clear than that for smaller . There are some apparent blurred parts in the left side of the inclusion in Figures 2(a) and 2(b). The blurred parts are much smaller in Figure 2(c) than in Figure 2(b) and almost disappeared in Figure 2(d).

From this example, we can find that the proposed algorithm works better for the bigger anisotropic factor . The following simulations in this paper are all performed with .

5.2. Simulation 2: Reconstruction with TV- Regularization and the Standard TV Regularization

There are two groups of experiments in this section. In the first group of experiment, we design three types of inclusions to compare the performance of the - regularization and the regularization with noise-free synthetic data. In the first case, a small inclusion with radius 0.5 mm centered at (7.0 mm, 0.0 mm) is designed for and , respectively. In the second case, a middle inclusion with radius 2 mm centered at (5.0 mm, 0.0 mm) is designed for and , respectively. In the third case, a middle inclusion with radius 4 mm centered at (3.0 mm, 0.0 mm) is designed for and , respectively. The simulated true absorption coefficient can be seen on the first column of Figure 3. For the small inclusion case, we use a mesh of 772 nodes and 1484 triangular elements for the inverse problem and a mesh of 1267 nodes and 2440 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV- regularization and TV regularization are shown in Figure 3. For the middle inclusion case, we use a mesh of 575 nodes and 1080 triangular elements for the inverse problem and a mesh of 813 nodes and 2642 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV- regularization and TV regularization are shown in Figure 4. For the big inclusion case, we use a mesh of 583 nodes and 1096 triangular elements for the inverse problem and a mesh of 871 nodes and 3124 triangular elements for the forward problem, the simulated true absorption coefficient and reconstruction results with TV- regularization and TV regularization are shown in Figure 5. The values of the parameters are shown in Table 1. Since we use the same as the last section, we will not present it repeatedly in Table 1.

Observing from the results, we can see that, for small and middle inclusions, the - regularization performs better than the regularization in both localizing the location and quantifying the values; see Figures 3(b) and 4(b) for - regularization and Figures 3(c) and 4(c) for regularization. In fact, these results are reasonable and can be interpreted. The - regularization imposes both penalty and penalty on . The penalty tends to find edges and the penalty tends to find the sparse details of the inclusion. As can be seen from Table 1, for the big inclusion, we enhance the weight of penalty. When the inclusion gets smaller, we enhance the weight of penalty.

For the big inclusion, - regularization and regularization perform no big differences, but we still can see that the blurred parts in Figure 5(b) with - regularization are smaller than that in Figure 5(c) with regularization. We can see from Figures 5(b) and 5(c) that there are large area of blurs in the big inclusion case. It is reasonable in the sense that large area of absorption inclusion means more photons are absorbed in propagation process. Hence, we can alleviate this phenomenon by increasing the source-detector pairs. In other words, by increasing the measurements, one can alleviate the effect of absorption to some extent.

In the second group of experiments of this section, three types of complicated and multi-inclusions are designed for comparing the reconstruction results by - regularization, regularization, and regularization.

In the first case, we reconstruct two small circle inclusions centered at (7 mm, 0 mm) and (0 mm, 0 mm) with the same radius 0.5 mm. In both of the two inclusions, and . We use a mesh of 1277 nodes and 2488 triangular elements for the inverse problem and a mesh of 1861 nodes and 3616 triangular elements for the forward problem. The true distributions of and reconstruction results with various regularizations can be seen in Figure 6.

In the second case, we reconstruct three small circles centered at (7 mm, 0 mm), (0 mm, −7 mm), and (0 mm, 0 mm) with the same radius, 0.5 mm. In all the three inclusions, and . We use a mesh of 1671 nodes and 3264 triangular elements for the inverse problem and a mesh of 2141 nodes and 4176 triangulations for the forward problem. The true distributions of and reconstruction results with various regularizations can be seen in Figure 7.

In the third case, we design three inclusions. One big circle centered at (−5 mm, 2 mm) with radius 2 mm. Two small circles centered at (0 mm, −7 mm) and (7 mm, 0 mm) with the same radius 0.5 mm. In all the three cases, we take the same coefficient values and . We use a mesh of 1323 nodes and 2568 triangulations for the inverse problem and a mesh of 1809 nodes and 3512 triangulations for the forward problem. The true distributions of and reconstruction results with various regularizations can be seen in Figure 8.

In Table 2, for all the three cases, we present the parameter settings for the iterative procedure (α, β, η), the steps of the iterations , the relative residual error , and the relative solution error in the sense of norm corresponding to each regularization method. and are defined as

We say that in the first and the second case, compared with the inclusions near the boundary, the inclusion deep in the domain is difficult to identify. In the third case, it is more difficult to reconstruct the small inclusion than the big inclusion. From the reconstruction results, we can see that the - regularization can preferably identify the location and the value of the inclusions. In other words, the advantages with - regularization are indeed apparent over regularization and regularization; see Figures 6(b), 7(b), and 8(b) for the reconstructions results with - regularization. The LM type algorithm for the regularization converges slowly especially when the inclusions are very small; see Figures 6(c) and 7(c). When big inclusion is included, the LM for regularization can well identify it, but it still cannot perfectly reconstruct the value of the small inclusions; see Figure 8(c). The split Bregman algorithm for regularization is an efficient algorithm with high convergent rate; we can see this from the number of iterations in Table 2. But, in our simulations, although the algorithm converges quickly, it loses validity when the inclusions are very small; see Figures 6(d) and 7(d). In Figure 8(d), the big inclusion is identified, but we only can see two blurred circles at where the two small inclusions are located. Moreover, we can find that the reconstruction may be oversparsified by using the regularization.

5.3. Simulation 3: High Absorbing and Low Scattering Inclusions with TV-

As the last simulation, our purpose is to show the validity of our algorithm in the situation that high absorbing inclusion and low scattering inclusion are contained in the domain . Moreover, we investigated the reconstruction under different noise level to show the robustness of the proposed algorithm. The noises added to the exact data are by the following rule: , where is the signal-to-noise ratio; is a Gaussian random variable with zero mean and unity variation. The background values of the absorption and the scattering are and , respectively. The absorbing inclusion for which we set the absorption coefficient is which is centered at (3.5 mm, 3.5 mm) with radius 2 mm. The scattering inclusion for which we set the scattering coefficient is which is centered at (0 mm, −5 mm) with radius 2 mm. For the discretization, we use a mesh of 463 nodes and 856 triangulations for the inverse problem and a mesh of 681 nodes and 1288 triangulations for the forward problem. See the first row of the Figure 9(a) for the true distribution of . In Figure 9(b), the reconstruction of with exact data is presented. The parameter value in this case is chosen as , , and . Figure 9(c) presents the reconstruction with noise level , in which the parameter is taken as , , and . In Figure 9(d), the reconstruction with noise level is shown, in which the parameter value is taken as , , and . From the four reconstructions under four different noise levels, we can find that as the noise level increased, the quality of the reconstruction gets worse. Here, we take the noise level under , since the iteration will not converge if the noise level is bigger than which reflects the severe ill-posedness of the problem. We also find that in order to obtain the reconstruction results with the relative error in norm no more than the noise level on the synthetic data should be less than .

6. Summary

In this paper, forward and inverse problems of the radiative transfer equation are considered. An image reconstruction method based on the - regularization is proposed. The forward problem is discretized with the discontinuous Galerkin method on the spatial space and the finite element method on the angular space which both are implemented on the piecewise linear basis. We discretize the absorption and scattering coefficients on the piecewise constant basis. The minimization problem is solved by a Levenberg-Marquardt type method which is equipped with a split Bregman algorithm for the penalty. The adjoint method is used to compute the Jacobian matrix which is the discretized Fréchet derivative of the forward operator. We numerically compare the proposed reconstruction method with the other imaging reconstruction methods based on and regularizations. The simulation results show the validity and robustness of the proposed method.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This work is supported by the National Nature Science Foundation of China under Grants nos. 41474102 and 61307023.