Improved Reconstruction Quality of Bioluminescent Images by Combining SP3 Equations and Bregman Iteration Method

Bioluminescence tomography (BLT) has a great potential to provide a powerful tool for tumor detection, monitoring tumor therapy progress, and drug development; developing new reconstruction algorithms will advance the technique to practical applications. In the paper, we propose a BLT reconstruction algorithm by combining SP3 equations and Bregman iteration method to improve the quality of reconstructed sources. The numerical results for homogeneous and heterogeneous phantoms are very encouraging and give significant improvement over the algorithms without the use of SP3 equations and Bregman iteration method.


Introduction
As an emerging molecular imaging technique, bioluminescence imaging (BLI) is potentially well suited for early detection, clinical drug development and monitoring, and regeneration research [1][2][3][4][5]. erefore, this imaging modality has received increasingly intense research interest worldwide over the recent years.
To date, planar BLI is commonly used because of its ease of implementation and operational simplicity, but it also suffers from signi�cant limitations, including the low resolution, the lack of quanti�cation, and the incapacity of accurately providing depth information [6]. In contrast, bioluminescence tomography (BLT) could overcome these limitations by using accurate reconstruction algorithms coupled with theoretical models of photon propagation in biological tissues, providing higher resolution, quanti�cation accuracy, and depth information [7]. In comparing BLT to planar BLI, planar BLI is a qualitative analysis and BLT is a quantitative analysis [8]. erefore, scientists are now paying more attention to the advancement of BLT research.
e objective of BLT is to recover the unknown bioluminescent source distribution based on the noisy surface measurements Φ meas [6,7]. Indeed, the problem is also called the inverse problem. However, a major difficulty in recovering the bioluminescent source distribution is imposed by multiple scattering which occurs when light propagates through biological tissues. is makes the inverse problem severely ill-posed [7]. Furthermore, the number of recovered unknown source distributions is usually far more than the number of detected boundary measurements, that is, (in many cases, ). Hence, BLT is also a typically underdetermined problem. To obtain a meaningful solution, regularization techniques are usually adopted, which consist of solving the following constrained optimization problem [9]: where ( ) is a properly chosen regularization term, represents regularization parameter, and is a linear operator, typically formed by discretizing diffusion equation with �nite element methods [10].
When ( ) 2 2 , the above regularized problem becomes the popular Tikhonov regularization, which inherently provides smoothed solutions and therefore offers compromised accuracy in localizing bioluminescent sources [11]. Recently, 1 -regularized problems, that is, ( ) 1 , have received an increasing amount of attention in optical imaging, which allow high-quality images to be reconstructed from a small amount of boundary measurements [11][12][13][14]. However, 1regularized problems can sparsify the bioluminescent source distribution, which affects the quality of reconstructed images [13,15].
Furthermore, in order to obtain the matrix in (1), the diffusion approximation (DA) to radiative transfer equation (RTE) is widely used as the forward model for BLT reconstructions. Although the DA is one of the most important approximation methods in BLT [6][7][8][9][10][11], it suffers from some limitations [12][13][14]. Firstly, the scattering is dominated over absorption and secondly, the DA fails in modeling light propagation in the vicinity of those highly vascularized tissue parts [12][13][14]. erefore, the DA will introduce signi�cant error in some BLT cases [14]. In contrast, the RTE is widely accepted as an accurate model for light propagation in biological tissues. However, the use of the RTE as the forward model for BLT is oen not feasible due to the facts that analytical solutions cannot exist for biological tissues with spatially nonuniform scattering and absorption properties and the computation of numerical approximations for the solution is extremely time consuming [16,17]. A generalized delta-Eddington phase function was recently presented to simplify the RTE, and the more accurate solution was obtained relative to the DA [18,19]. However, the parameter used in the model is difficult to compute [18,19]. In addition, the system matrix for the model is also difficult to construct for complex heterogeneous geometries. ese factors seriously limit the utilization of the model in BLT. e use of simpli�ed spherical harmonics (SP N ) equations to approximate the RTE has been demonstrated to signi�cantly improve the diffusion solution in domains with high absorption and small geometries [5, 12-14, 16, 20]. Meanwhile, the SP N methods are computationally less expensive than the RTE ones.
Large efforts in combining multiple types of a priori information to develop BLT reconstruction algorithms to improve the quality of reconstructed images, particularly the permissible source region and multispectral information, have formed the grounds of BLT reconstructions [9][10][11][20][21][22][23][24][25][26]. Despite the recent advances in BLT reconstruction algorithms and light propagation models, it is necessary to develop and re�ne reconstruction methods to improve image quality.
Bregman iteration method has been studied recently and is widely used in compressed sensing [27,28]. e idea is to add the residual, that is, the error produced at the current iteration, back to the data for the next iteration to be corrected [27]. e method is particularly attractive for sparse reconstruction, but so far it has not been fully investigated and analyzed in BLT, and this is the goal of this paper.
In this paper, we propose a BLT algorithm to improve the quality of reconstructed images. In the algorithm, SP 3 equations are adapted to model light propagation, and Bregman iteration method is used to solve the inverse problem for BLT. Numerical results demonstrate that the quality of reconstructed images is improved greatly. e rest of the paper is organized as follows. In the following section, we described SP 3 equations as light propagation model and Bregman iteration method. Last, numerical experiments were performed to evaluate the proposed algorithm, and corresponding conclusions were made. 3 Equations. e propagation of light in biological tissues can be well modeled by SP 3 equations. SP 3 equations are two coupled diffusion equations for the moments 1 and 2 [16,17]:
e boundary conditions are given by e coefficients 1 , … , 1 , … , 2 , … , 2 can be found in [16]. Furthermore, the partial current can be obtained from solutions 1 and 2 : e coefficients 0 , 1 , … , 3 can also be found in [16]. Solving the above equations by �nite element methods, a linear operator can be established [29].

Bregman Iteration Method.
Bregman iteration method is based on the de�nition of Bregman distance. e Bregman distance associated with a convex function at the point is given as [27] ( , ) = ( ) − ( ) − , − , where is in the subgradient of at . Clearly, this is not a distance in the usual sense because it is not in general symmetric. However, it does measure closeness in the sense that ( , ) 0 and ( , ) ( , ) for on the line segment between and [27].
where ∈ ) and is the adjoint operator of . Since the operator is linear in BLT reconstructions, the above complicated iteration can be transformed to the following two-stage iteration procedure with [27]: is is done by iteratively solving the optimization problem (7) and then modifying the measured value of meas used in the next iteration. And (8) is usually referred as "adding back the noise" [30]. In the paper, ) is �xed as the regularizer. e implementation of (7) was performed by a gradient projected (GP) algorithm [31]. e proposed algorithm was depicted in Algorithm 1.

Results
To fully evaluate the performance of the proposed algorithm, homogeneous and heterogeneous experiments were performed. In the experiments, the parameters and max were set to × −3 and 10, respectively. e parameters in GP algorithm set default values, except the maximum iteration number is �xed at 50000 to ensure the convergence of the algorithm unless otherwise is speci�ed.

Homogeneous Phantom Experiments.
In this section, 2D numerical simulations were used to investigate the performance of the proposed algorithm since less computational time was required for 2D data. Here, two individual cases were considered. In the �rst case, numerical simulations were performed on a homogenous circle with 10 mm radius. Within this circle, two sources (source 1 and source 2) were placed in (−5, 0) mm and (0, 5) mm, respectively and each T 1: Optical properties for different bands [22]. source had a radius of 1.0 mm. e corresponding optical parameters were listed in Table 1. e boundary data were generated for two wavelengths (600 and 620 nm) with �nite element methods, and different levels of Gaussian noise (0%, 10%, and 30%) were added to the datasets. BLT reconstructions were performed without and with Bregman iteration method. Corresponding results were shown in Figure 1. In this case, the ratios of ′/ are larger than 10; therefore, the circular phantom has high-scattering characteristics. Hence, the DA is suitable for the simulation. For comparison, we carried out BLT reconstructions with the DA as the forward model; reconstructed images were also illustrated in Figure 1. From Figure 1, we can see that the results with SP 3 equations are better than those obtained with the DA and Bregman iteration method can improve the quality of reconstructed images. e best results are obtained by combing SP 3 equations and Bregman iteration method. In addition, quantitative results were summarized in Table 2. Data in Table 2 show that reconstructed position errors can be signi�cantly reduced when SP 3 equations are used together with Bregman iteration method Furthermore, we tested the proposed algorithm by using experiments with multiple bioluminescent sources. e optical properties of a real mouse muscle for different wavelengths (580 and 620 nm) were assigned as listed in Table 3 [29]. Four identical sources with 1 mm radii were placed different positions. First, the sources were placed near the surfaces, and the distance to the center of the circle was 7.07 mm. e boundary measurements were also produced by �nite element methods, and 20% Gaussian noise was added into the simulated data. Note that in the test, ′ / for two wavelengths are less than 10; therefore, the condition ′ ≫ does not hold and the DA is less valid. Hence, BLT reconstructions with the DA were not implemented. e results with SP 3 equations are shown in Figures 2(a) and 2(b). Next, the sources were placed at 5 mm positions off the center. en BLT reconstructions were performed, as shown in Figures 2(c) and 2(d). Furthermore, quantitative results were shown in Table 4. It is worthy of mentioning  that BLT reconstructions without and with Bregman iteration method use the same regularization parameter (i.e., 3×10 −6 ), but the reconstructed results are different. From Figure 2 and Table 4, it is easily concluded that better images can be obtained by combining SP 3 equations and Bregman iteration method.

Heterogeneous Phantom.
In the subsection, a micro-MRI-based heterogeneous mouse model (MOBY) was used to validate the proposed algorithm [32]. About 2/3 of the entire phantom was used for mesh generation, and a volumetric mesh with 17661 nodes and 93312 tetrahedron elements was obtained by iso2mesh [33], as shown in Figure 3.  e optical properties of different tissues were assigned according to Table 5, reproduced from Alexandrakis et al.
[21�. e forward simulation data was produced by �nite element methods, and 10% Gaussian noise was added. en BLT reconstructions were performed without and with Bregman iteration method. e regularization parameters used in the two methods were the same, and the value was 0.1. e maximum iteration number in the GP algorithm was set to 5000, and other parameters remained unchanged. e reconstructed results without and with Bregman iteration method were shown in Figure 4. From the images, we can see that the quality of reconstructed images can be improved with the

Conclusion
We have presented a BLT reconstruction algorithm by combing SP 3 equations and Bregman iteration method as a competitive method for reconstructing bioluminescent sources and validated the proposed algorithm using homogeneous and heterogeneous experiments. It has been demonstrated that the proposed algorithm can enhance the recovery of bioluminescent sources in terms of the quality of reconstructed images and localization error. e use of SP 3 equations is a helpful technique to improve BLT reconstructions. Our experiments have illustrated that the appearance of artifacts can be reduced when SP 3 equations are used as the forward model. However, the computation of the system matrix by solving SP 3 equations is very expensive, especially when the imaged objects are very complex, irregular, and heterogeneous. Fortunately, with the fast development of graphics processing unit (GPU), the computation of can be signi�cantly accelerated.
One merit of the proposed algorithm is that the improved results are obtained by making use of the available boundary measurements and thus do not require increased number of boundary measurements and do not bring more hardware requirements. Meanwhile, the proposed algorithm is relatively easy to implement. erefore, the algorithm is suitable for in vivo applications. As a sacri�ce, the computational burden for the proposed algorithm is greatly increased, especially for the heterogeneous mouse experiment, since solving (1) brings extra cost through Bregman iteration method, and each iteration of which is equivalent of solving a standard "L1" problem. To increase computational efficiency for mouse experiments, developing fast large-scale optimization algorithms is essential.
In conclusion, we have developed a BLT reconstruction algorithm by combing SP 3 equations and Bregman iteration method and indicated its feasibility and merits. In the near future, we expect to accelerate the proposed algorithm based on GPU and extend it to in vivo mouse experiments.