Computational and Mathematical Methods in Medicine

Volume 2017, Article ID 2953560, 15 pages

https://doi.org/10.1155/2017/2953560

## Mixed Total Variation and Regularization Method for Optical Tomography Based on Radiative Transfer Equation

^{1}Department of Mathematics, Harbin Institute of Technology, Harbin, Heilongjiang Province, China^{2}Department of Mathematics, University of Iowa, Iowa City, IA, USA^{3}School of Mathematics and Statistics, Northeast Petroleum University, Daqing, Heilongjiang Province, China

Correspondence should be addressed to Bo Han; nc.ude.tih@nahob

Received 22 March 2016; Accepted 27 September 2016; Published 9 February 2017

Academic Editor: Chuangyin Dang

Copyright © 2017 Jinping Tang 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

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 [1–4]. 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 [6–9]. 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 [12–14] and the references therein. To date, TV regularization has been the preferred regularization strategy of considerable recent works, like photoacoustic tomography (PAT) [15–17], 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, [26–30]. 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.