1. Introduction

In the past decade, bioluminescence imaging (BLI) has been one of the hot topics in optical imaging and has been successfully applied in tumorigenesis studies, cancer diagnosis, metastasis detection, drug development, gene therapy, etc. [1–5]. However, the intrinsic planar imaging mode of BLI is just a qualitative imaging tool, which provides no absolute quantitative information and detects limited depths of a few millimeters inside tissues [4,5]. The bioluminescence tomography (BLT), the three-dimensional (3D) tomographic counterpart of BLI, has received much attention because it can provide localization and quantitative information of the bioluminescent source inside a living small animal [6–9].

In BLT, the tumor cells of the living mouse model are usually transfected with a luc-gene. Then, the oxidation of the substrate and luciferin results in a red-shifted light emission with a wavelength in the range from 500 nm to about 760 nm [3]. Based on the surface bioluminescent images and the heterogeneous tissue structure information captured by, for example, the micro-CT imaging system, BLT can reconstruct the location and distribution of the internal embedded bioluminescent light source and acquire quantitative volumetric reconstruction information [8].

The inverse source reconstruction is crucial for recovering three dimensional information of the light source. Many efforts have been devoted to studies on efficient inverse reconstruction [10–28]. Generally, some prior knowledge, such as the permissible source region and multispectral information, are used to reduce the ill-posedness of inverse reconstruction [10–14,16,17]. Meanwhile, the development of a fast and robust algorithm is another challenging problem. The common schemes include the regularization method, iterative method, stochastic method and so on [15,18–26]. Since the compressive sensing (CS) theory was introduced in 2006, which is a signal processing technique for efficiently acquiring and reconstructing images with sparse representation from far fewer measurements or signals below the Nyquist sampling theorem demands, it has been utilized by various areas of applied mathematics, computer science, electrical engineering and medicine [29,30]. As a new imaging technique, the sparse bioluminescent source and insufficient measuring signal are the characteristics of BLT, which can benefit by using the CS technique. Along with the development of the compressive sensing (CS) theory, the l_p regularization method has become the mainstream strategy to obtain the optimal solution of the BLT inverse problem [17,22,24–26].

Generally speaking, regularization is used to construct the approximation of an ill-posed problem by a family of neighboring well-posed problems [31]. So far, the regularization methods have been used for solving various types of inverse problems in mathematics, statistics, machine learning, signal and image processing, heat conducting, and so on. The most commonly used form of regularization is the Tikhonov-type regularization. This method achieves a compromise between the minimization of the residual norm and the penalty term [32]. The l₂-norm is the integrant penalty term for the Tikhonov-type regularization. The l₂-norm penalty term imposes a small weight on small residuals, but a large weight on others. As a result, the l₂-norm approximation has many modest residuals and relatively few larger ones, which generates smooth solutions and brings multipseudo sources to surround the true source in BLT reconstruction [22]. Considering the source sparsity characteristics and the insufficiency of the measured data, the l₁-norm penalty term has been adopted in the regularization method and was successfully applied to BLT [24–26,33]. Unlike l₂-norm regularization, l₁-norm regularization tends to produce a sparse solution with a large number of components equal to zero. However, in many practical applications, the solutions of the l₁ regularizer are often less sparse than those of the l₀ regularizer. The l₀ regularizer aims to find the number of non-zero elements in the solving domain, but it is an impossible mission in practical applications [34]. To find more sparse solutions than the l₁ regularizer, l_p (0 < p <1) regularization methods have been developed and employed to solve the minimization problems [34–36] and the BLT reconstruction problem [28].

Although the l_p (0 < p ≤ 1, p = 2) regularization methods have been applied to the BLT reconstruction problem, so far, there is no generally accepted method which can be suitable for all of the reconstruction cases. This paper intends to fill the gap in the existing studies to systematically benchmark the performance of the l_p-regularization-based BLT reconstruction algorithms. Our group has been conducting research in the BLT reconstruction method for years. In this paper, we evaluated six reconstruction algorithms, including the Truncated Total Least Square method (TTLS) [27], the Incomplete Variables Truncated Conjugate Gradient method (IVTCG) [25], the Truncated Newton Interior-Point method (TNIPM) [24], the Primal-Dual Interior-Point method (PDIP) [26], and the Weighted Iterative Shrinkage/Thresholding algorithm (WISTA) [28] that were developed by our group, and the Tikhonov regularization that we used for the BLT reconstruction in its standard way [31]. In order to investigate the response of these algorithms to the permissible source region, measurement noise, optical properties, and tissue specificity, we conducted a series of single source numerical phantom experiments. Then, the double source numerical phantom experiment and the in vivo mouse experiment were carried out to further test their performance.

This paper is organized as follows. The next section describes the diffusion approximation model of BLT and the inverse reconstruction formula. The essence of the six regularization algorithms are then introduced briefly. In Section 3, the response of these algorithms to several influence factors are investigated and discussed based on the numerical phantom and in vivo mouse experiments to demonstrate the performance of the different algorithms. Finally, the discussion and conclusions are addressed in Section 4.

2. Methodology

2.1. Inverse reconstruction formula of BLT

By assuming light as photons, the propagation of light in biological tissues is an extremely complex process that involves a series of optical behaviors, including the absorption, scattering, internal reflection and transmission. The radiative transfer equation (RTE) is often used to analytically model the process of light propagation. Although the accuracy of the RTE model is indubitable, it is generally difficult to directly solve it in complex biological tissues. Furthermore, the complicated implementation in the numerical settings and the large numbers of the resulting equations makes its efficiency low for BLT [37]. In infrared optical imaging, most of the biological tissues exhibit the characteristics of high scattering and low absorption. As the first order approximation model of RTE, the diffusion approximation (DA) is widely used to depict the process [8]:

where $S\left(r\right)$ is the source distribution, r is the position vector, $\Phi \left(r\right)$ is the photon flux density, ${\mu }_a(\text{r})$ is the absorption coefficient, $D(\text{r})=1/(3({\mu }_a(\text{r})+(1-g){\mu }_s(\text{r})))$ is the diffusion coefficient, ${\mu }_s(\text{r})$ is the scattering coefficient, g is the anisotropy parameter, and $\Omega \subset R^3$ is the domain of interest. Equation (2) is the Robin boundary condition, and $v\left(r\right)$ is the unit outer normal to $\partial \Omega$ at r, $A_n(r)=(1+R(r))/(1-R(r))$ is the boundary mismatch factor and $R(r)$ depends on the refractive index n of the surrounding medium that can be approximated by $R\approx -1.4399n^{-2}+0.7099n^{-1}+0.6681+0.0636n$ [38]. In the optical imaging experiment, the outgoing flux density on $\partial \Omega$ is measured by the highly sensitive CCD camera and can be expressed as follows [8]:

To solve Eqs. (1) and (2), the finite element method (FEM) was adopted and the equations were discretized by the finite element basis functions. A linear relationship between the measurable boundary data b and the unknown source distribution S is then established and the BLT forward problem can be obtained as [39]:

where the weight matrix A establishes the relationship between the measurements and the unknowns.

Because the measurement is insufficient compared with the whole solving domain, the inverse reconstruction is an ill-posed problem. As a result, it is difficult and impossible to solve Eq. (4) directly. Usually, the regularization method is employed to convert Eq. (4) into the following optimization problem and S is replaced by x to make the notation more general:

where ${\Vert x\Vert }_{\text{p}}={\left({\displaystyle \sum _{i\text{=}1}^n{\left|x_i\right|}^{\text{p}}}\right)}^{1\text{/p}}$ and p denotes the parameter used in the l_p regularization and satisfies 0 < p ≤ 1 or p = 2. Selecting an appropriate optimization method, the localization and distribution of the bioluminescent probes can be obtained from Eq. (5).

2.2. Regularization methods

2.2.1. Tikhonov regularization

Tikhonov regularization method is one of the most commonly used methods to solve the ill-posed inverse problems. When p = 2, Eq. (5) is converted into the following minimization problem:

where $\text{λ>0}$ is a properly chosen regularization parameter and controls the weight given to the minimization of the regularization penalty term relative to the minimization of the residual norm term. In the BLT application, the right-hand side b is always contaminated by various types of errors, such as the measurement errors, approximation errors, and rounding errors. Assuming the right-hand side b satisfies the discrete Picard condition, an explicit solution of (6) is given by

By using the singular value decomposition (SVD) explicitly, the system matrix A can be expressed as $A=U\Sigma V^T$ and the Tikhonov regularized solution is given by

where $U=(u_1,\mathrm{....},u_n)$ and $V=(v_1,\mathrm{....},v_n)$ are matrices with orthonormal columns, $U^{T}U=V^{T}V=I_n$, and where $\Sigma =diag({\sigma }_1,\mathrm{...},{\sigma }_n)$ has non-negative diagonal elements appearing in non-increasing order such that ${\sigma }_1\ge \cdots \ge {\sigma }_n\ge 0$ [40].

2.2.2. Truncated Total Least Squares Method

The Tikhonov regularization method only considers the measurement errors on the right-hand side b, but the system errors also exist in the computed coefficient matrix of the system equation. In fact, whether the system errors come from the coefficient matrix or the measurement noise, they are inevitable in a real BLT experiment, so He et al. proposed the truncated total least squares (TTLS) method combined with a practical parameter-choice scheme termed as improved generalized cross-validation (IGCV) to cope with the BLT reconstruction problem including both the measurement noise and the system errors [27].

The TTLS method takes its basis in an SVD of the augmented matrix $(Ab)=\overline{U}\overline{\Sigma }{\overline{V}}^T={\displaystyle \sum _{i=1}^{n+1}{\overline{u}}_i}{\overline{\sigma }}_i{\overline{v}}_i^T$, $\overline{U}\in {\text{R}}^{m\times (n+1)}$, $\overline{V}\in {\text{R}}^{(n+1)\times (n+1)}$, where ${\overline{U}}^T\overline{U}={\overline{V}}^T\overline{V}=I_{n+1}$, and $\Sigma =diag({\overline{\sigma }}_1,\mathrm{...},{\overline{\sigma }}_{n+1})$ have singular values ${\overline{\sigma }}_1\ge \cdots \ge {\overline{\sigma }}_{n+1}\ge 0$. The matrix $\overline{V}$ is partitioned such that $\overline{V}=\left(\begin{array}{cc}{\overline{V}}_{11}& {\overline{V}}_{12}\\ {\overline{V}}_{21}& {\overline{V}}_{22}\end{array}\right)$, ${\overline{V}}_{11}\in {\text{R}}^{n\times k}$, and ${\overline{V}}_{22}\in {\text{R}}^{1\times (n+1-k)}$, where k is the truncation parameter. Then, the TTLS solution is given by $x_{\text{TTLS}}=-{\overline{V}}_{12}{\overline{V}}_{22}^{†}=-\frac{{\overline{V}}_{12}{\overline{V}}_{22}^T}{{\Vert {\overline{V}}_{22}\Vert }_2^2}$. Based on the modified GCV (MGCV) criterion for TTLS $G=\frac{{\Vert Ax-b\Vert }^2}{{(m-{\text{enp}}_k)}^2}$, where ${\text{enp}}_k$is the sum of the filter factors of TTLS, the truncation parameter k is obtained by the IGCV scheme.

2.2.3. Incomplete Variables Truncated Conjugate Gradient (IVTCG) Method

In view of the insufficient measurement, the source sparse characteristic and the oversmoothing of l₂-norm regularization in BLT, linear system Eq. (4) is covert to the following l₁-norm regularization problem based on the CS theory:

where ${\Vert x\Vert }_1={\displaystyle {\sum }_i{x}_i}$ is the l₁ norm of x. Introducing vectors u and v, and the substitution $x=u-v$, Eq. (9) can be reformulated as the following convex quadratic program with nonnegative constrained conditions:

where $u\ge 0$, $v\ge 0$, $z=[u{v]}^T$, $c=\lambda 1_{2N}+[-q{q]}^T$, $1_{2N}={\left[1,1,\cdots ,1\right]}^T$, $q=A^{T}x$, and $B=[A{-A]}^T[A-A]$. Then, the incomplete variables truncated conjugate gradient (IVTCG) algorithm was presented by He et al. to solve Eq. (10) [25]. First, we define the working set ${\Gamma }_k$ by the violation of the Karush-Kuhn-Tucker (KKT) optimal conditions at the kth iteration,

where $z^k$ is the current iteration point. Then, we define the other two working sets $I_k$ and $J_k$ as follows:

where $\delta >0$, $N_s>0$and $N_{\max }>N_s$are three constants, ${\widehat{J}}_k={\Gamma }_k\backslash I_k$ and ${\widehat{I}}_k=\left\{i\in \left\{1,\cdots ,2N\right\}|\left[{\left(z^k\right)}_i>0,{\left(z^k\right)}_i/{\left(\nabla F(z^k)\right)}_i>\delta \right]\right\}$. The variables with the set will be updated $I_k$and $J_k$ at the current iteration.

Finally, we need to decide the search direction. For the search direction $d_{I_k}^k$of $I_k$, it can be obtained by solving the sub-problem

For $J_k$, the search direction $d_{J_k}^k\text{=}-w_{J_k}$. The other entries of $d^k$ are set to 0. After determining the search direction, the next iteration point is derived by the backtracking method.

In the IVTCG method, the parameters $N_s$and $N_{\max }$are set as $N_s=\lfloor \raisebox{1ex}{$M$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\rfloor$and $N_{\max }=N_s+\lfloor \raisebox{1ex}{$N_s$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\rfloor$. The time complexity of the IVTCG algorithm is $O(M{m}^2)+O(MN)$.

2.2.4. Truncated Newton interior-point method (TNIPM)

Using the truncated Newton interior-point method (TNIPM), the object function (9) is transformed to a convex quadratic problem with inequality constraints:

Then, a logarithmic barrier function for the constraints $\left|x_i\right|\le u_i$ is added to the minimization function (14), and Eq. (14) is converted to the following problem:

where parameter $t\in (0,\infty )$.

In order to compute the search direction, a preconditioned conjugate gradient method is adopted to compute the Newton system

Then, we can obtain a feasible point by using the initial value and the search direction. The implementation of TNIPM in BLT reconstruction is presented in [24].

2.2.5. Primal-dual interior-point method (PDIP)

As an efficient method for solving the optimization problems, interior-point methods have been applied in many areas in the past twenty years. In [26], the primal-dual interior-point method was adopted to solve the BLT inverse problem. Equation (4) is transformed to the following l₁-norm minimization problem:

and then converted to the standard primal and dual forms in linear programming:

By introducing a logarithmic barrier term, (P) is converted to the following logarithmic barrier problems:

where θ is the barrier parameter. Then, we can obtain the Karush-Kuhn-Tucker conditions for ${\text{P}}_{\theta }$:

where $e={\left(1,\mathrm{...},1\right)}^T$, the matrix X and S are the $n\times n$ diagonal matrix whose diagonal entries are precisely the components of x and s respectively.

Hence, a new feasible point and the Newton direction $\left(\Delta x,\Delta y,\Delta s\right)$ at the current value $\left(x^k,y^k,s^k\right)$can be obtained by solving the following equation system:

The general framework of PDIPM is presented in [26].

2.2.6. Weighted iterative shrinkage/thresholding algorithm

In some cases, the l₁-norm regularization methods are often less sparse than l_p-norm (0 < p < 1) regularization methods. In the paper [28], we proposed a weighted iterative shrinkage/thresholding algorithm (WISTA) to solve the l_p-norm minimization problem (5) (0 < p < 1). Because we cannot solve the l_p-norm minimization problem directly, we reformed it by the weighted l₁-norm:

In order to obtain an ideal solution of (23), we used the iterative form as follows:

where the weight vector ω is calculated from the last iterate $x^i$, soft(z, q) is the soft-threshold function with a threshold q, and $\text{soft(}x\text{, }q\text{)=sign(}x\text{)max(}\left|x\right|\text{-}q\text{, 0)}$. For the implementation of WISTA, we set two stopping conditions. The first is that the quantity of descending in the objective function is less than the stopping threshold. The second is that the algorithm reached the maximum number of iterations.

The final algorithm of WISTA is formulated in [28].

3. Experiments and results

In this section, a series of numerical experiments and an in vivo experiment were conducted to benchmark the performance of the six investigated reconstruction algorithms. We used two indicators to evaluate the reconstruction results: the reconstruction time (Time) of these algorithms and the distance (Loc_Err) between the actual source position and the reconstruction position. We defined the reconstruction position as the point with the maximum density of the reconstructed source. The regularization parameter of the Tikhonov method was chosen by the classical L-curve method and the parameter-choice scheme of the TTLS method was the improved generalized cross validation (IGCV). The regularization parameters of TNIPM and IVTCG used in reconstruction were manually optimized and they were set as 1e-6 in this paper. All experiments were performed in MATLAB on a personal computer with 2.66 GHz Intel Core 2 Quad CPU and 4 GB RAM.

3.1. Numerical Experiments

3.1.1. Reconstruction in a single source

In this part, a cylindrical heterogeneous phantom was employed to demonstrate the features of the investigated algorithms. The numerical phantom consisted of muscle, lungs, heart, bone and liver, with the dimensions of 30 mm in height and 20 mm in diameter, as shown in Fig. 1 [39]. The optical parameters of the five tissues were specified as listed in Table 1 . A solid sphere source with a 1 mm radius was centered at (3, 5, 0) inside the right lung. Meanwhile, 3874 nodes and 17763 tetrahedral elements were used in the inverse reconstruction. In order to compare the performance of the six algorithms among various conditions, we designed five experiments for the single source reconstruction.

 

Fig. 1 Heterogeneous phantom with a single light source composed of muscle, lungs, heart, bone, liver and the source in the right lung.

Download Full Size | PDF

Table 1. Optical parameters of the heterogeneous phantom [39]

View Table | View all tables in this article

1. Reconstruction using different permissible source regions

As a priori knowledge, the permissible source region (PSR) is commonly used in the inverse reconstruction of BLT. Different PSR strategies may lead to different reconstruction results. In Fig. 2 , we can see the surface flux distribution of the heterogeneous phantom in the coronal, sagittal and translucency views respectively and six red lines indicate the PSR. We set four different PSRs to verify its influence on the inverse reconstruction of the six algorithms. The PSRs were defined as follows: PS1 = {(x, y, z)| the left lung}; PS2 = {(x, y, z)| 0<x, y<7, −5< z <5}; PS3 = {(x, y, z)| −5<x, y<8, −5< z <5}; and PS4 = {(x, y, z)| the whole phantom}. The volume ratios of the four PSR were 3.26%, 5.20%, 17.93% and 100% respectively. The reconstruction results are listed in Figs. 3 ­6. Since TTLS could only resolve the overdetermined equation, it did not provide results for PS4. In Fig. 6(b), because the reconstruction time of PDIP in PS4 was too long so as to make the time curves almost indistinguishable in the Cartesian coordinate system, we used the log plot to depict the reconstruction time of various methods.

 

Fig. 2 The surface flux distribution of the heterogeneous phantom from different views. (a)-(b) The coronal and sagittal views, respectively. (c) The translucency view of (a).

Download Full Size | PDF

 

Fig. 3 Axial views of the BLT reconstruction results of l₂ regularization methods using different PSRs at z = 0mm. (a)-(c) Results of TTLS; (d)-(g) Results of Tikhonov.

Download Full Size | PDF

We can draw the conclusion that the smaller the PSR, the better the reconstruction results obtained. From Figs. 3–5, we found that the reconstruction results of l_p (0 < p ≤ 1) regularization methods (TNIPM, IVTCG, PDIP, WISTA) were better than those of l₂ regularization methods (TTLS, Tikhonov). From Fig. 3, the reconstruction results of l₂ regularization methods were scattered around the true source. On the contrary, from Fig. 4 and Fig. 5, the reconstruction results of l_p (0 < p ≤ 1) regularization methods were more accurate and the reconstruction sources were smaller than those of l₂, except for the results of PDIP.

 

Fig. 4 Axial views of the BLT reconstruction results of l₁ regularization methods using different PSRs at z = 0mm. (a)-(d) Results of TNIPM; (e)-(h) Results of IVTCG; (i)-(l) Results of PDIP.

Download Full Size | PDF

 

Fig. 5 Axial views of the BLT reconstruction results of WISTA using different PSRs at z = 0mm.

Download Full Size | PDF

From Fig. 6, the smallest location error was obtained by PDIP for PS1 and the second-smallest location error was obtained by IVTCG for PS2. As the permissible source region became larger, the location error of PDIP increased. The worst result of IVTCG was obtained for PS1 and it may be due to the fact that IVTCG couldn’t develop its capacity of sparse reconstruction in a small PSR. The reconstruction results of TNIPM and WISTA remained stable with different PSRs. The reconstruction times of WISTA and IVTCG were smaller than the remaining algorithms and PDIP was the most time consuming algorithm. On average, the reconstruction locations of IVTCG, TNIPM and WISTA outperformed the other algorithms.

 

Fig. 6 Performance metrics for the six algorithms using different PSRs. (a) The distance errors of the BLT reconstruction results; (b) The reconstruction time of various methods.

Download Full Size | PDF

2. Reconstruction at different measurement noise levels

In this experiment, random noise at different levels was added to the synthetic data to evaluate the stability and robustness of the six algorithms using the permissible source region of PS1. The reconstruction results under different noise levels (0%-50%) of the Gaussian noise added to the synthetic data are shown in Fig. 7 and Table 2 [41]. The noise was generated by a MATLAB function randn and each algorithm was tested with a different noise sample. Tikhonov, TNIPM, IVTCG and WISTA were stable at all noise levels. The location error of TTLS increased a little bit at the noise level of 10%, and then it was kept stable until the noise level of 40%. The reconstruction results of PDIP behaved similarly. At the noise level of 50%, both the TTLS and PDIP methods became unstable, so we ran them several times and just showed the best results in Fig. 8(a) and Fig. 8(e) respectively. In Fig. 7(b), TNIPM has a larger computation cost and IVTCG methods are fast but produce larger localization errors.

 

Fig. 7 Performance metrics of various algorithms at different noise levels. (a) The distance errors of the BLT reconstruction results; (b) The reconstruction time of various methods.

Download Full Size | PDF

Table 2. Reconstruction results at different measurement noise levels

View Table | View all tables in this article

 

Fig. 8 Axial views of the BLT reconstruction results of six regularization methods at the 50% measurement noise level at z = 0mm. (a)-(f) are the results of TTLS, Tikhonov, TNIPM, IVTCG, PDIP and WISTA respectively.

Download Full Size | PDF

3. Reconstruction using different optical parameters

In order to investigate the influence of optical property on inverse reconstruction, four experiments with different deviations of absorption and scattering coefficients were conducted by adding 10% perturbation to the absorption and scattering coefficients of different organs. The corresponding results are shown in Table 3 and Fig. 9 . Therein, the symbols were defined as: P1 = {10% deviation of ${\mu }_a$ and ${\mu }_s$}; P2 = {-10% deviation of ${\mu }_a$ and ${\mu }_s$}; P3 = {10% deviation of ${\mu }_a$ and −10% deviation of ${\mu }_s$}; and P4 = {-10% deviation of ${\mu }_a$and 10% deviation of ${\mu }_s$}. From these results, we found that the results of P4 were better than those of the others and the best reconstruction location was obtained by IVTCG in P4. The reason should be that the optical parameters exhibit the characteristics of high scattering and low absorption in P4, which conforms to the application scope of the diffusion approximation model. From Fig. 9, IVTCG achieved the smallest location error and WISTA had the least amount of computation time. The optical property changes had little impact on the inverse reconstruction of WISTA and TNIPM. The l_p (0 < p ≤ 1) regularization methods (TNIPM, IVTCG, PDIP, WISTA) still obtained better reconstruction sources than the l₂ regularization methods (TTLS, Tikhonov).

Table 3. Reconstruction results using different optical parameters

View Table | View all tables in this article

 

Fig. 9 Performance metrics for the six algorithms using different optical parameters. (a) The distance errors of the BLT reconstruction results; (b) The reconstruction time of various methods.

Download Full Size | PDF

4. Reconstruction in tissue specificity

In this part, the influence of tissue specificity on the inverse reconstruction was observed in four experiments. In the first case of N1, the absorption, scattering and anisotropy coefficients of all organs were set the same as that of muscle to simulate a 3D homogeneous model. In the cases of N2-N4, we simulated three heterogeneous models (see Fig. 10 ). In the case of N2, the absorption, scattering and anisotropy coefficients of the lungs, heart and liver were set the same as those of muscle, thus the model was divided into two parts: muscle and bone. In the case of N3, the absorption, scattering and anisotropy coefficients of the lungs and heart were set the same, that is, the model was divided into four parts. The case N4 was the PS1 in the PSR experiment and there were five tissues in the model. The reconstruction results are presented in Table 4 and Fig. 11 . From these results, we found that the reconstruction results of the case N4 were better than the other three cases for l_p (0 < p ≤ 1) regularization methods (TNIPM, PDIP, WISTA), but this could not be observed for the l₂ regularization algorithms (TTLS, Tikhonov) and IVTCG. It may be that there were more optical parameters in tissues in N4, so it was more approachable in reality than the other three cases. As shown in Fig. 11(a), the PDIP’s performance was the best in the four cases and the reconstruction results of WISTA ranked second. In Fig. 11(b), WISTA still used the minimal amount of time, followed by IVTCG and Tikhonov.

 

Fig. 10 Tissue specificity models (various colors are for various segmented tissues).

Download Full Size | PDF

Table 4. Reconstruction results for tissue specificity

View Table | View all tables in this article

 

Fig. 11 Performance metrics for various algorithms of different tissue specificity. (a) The location errors of BLT reconstruction; (b) The reconstruction time of various methods.

Download Full Size | PDF

5. Reconstruction at different source locations

In this experiment, we carried out three experiments to observe the influence of the source locations on the inverse reconstruction. In case 1, the source location was (4,6,0) and the distance between the source center and the surface of phantom was 2.789mm. In case 2, that is PS1, the source location was (3,5,0) and the distance between the source center and the surface of phantom was 4.169mm. In case 3, the source location was (2,4,0) and the distance between the source center and the surface of phantom was 5.527mm. The reconstruction results are listed in Table 5 and the error bar can be seen in Fig. 12 . From the results, we observed that the smaller the depth was, the better the reconstruction results obtained. The standard deviation of WISTA was smallest indicating that it’s robust to source locations. Except for PDIP, the l_p (0 < p ≤ 1) regularization methods (TNIPM, IVTCG, WISTA) obtained smaller standard deviations than l₂ regularization methods (TTLS, Tikhonov).

Table 5. Reconstruction results at different source positions

View Table | View all tables in this article

 

Fig. 12 Error bar chart of the Loc_Err at different source positions.

Download Full Size | PDF

3.1.2. Reconstruction of double sources

In this part, three experiments were carried out to verify the performance of these algorithms on the double source resolution. In experiment 1, two sources were located in the left lung with their centers at (3, 5, 1.25) and (3, 5, −1.25), and with a distance of 2.5mm. In experiment 2, two sources were located in the left lung with their centers at (3, 5, 1. 5) and (3, 5, −1.5), and with a distance of 3mm. In experiment 3, two sources were located in the left lung with their centers at (3, 5, 2) and (3, 5, −2), and with a distance of 4mm. The right lung was specified as the PSR in the three experiments. The reconstruction results are presented in Figs. 13 –15 and Table 6 . From these results, we found that TTLS only reconstructs one source in case 1 and case 3, as shown in Fig. 13(a) and Fig. 13(c). In Fig. 13(f), the Tikhonov’s reconstruction of the two sources was almost overlapping. On the contrary, l_p (0 < p ≤ 1) regularization methods (TNIPM, IVTCG, PDIP, WISTA) obtained much better results. Comparing Fig. 14 with Fig. 15, we found that the reconstruction sources of WISTA were smaller than that of the three l₁ regularization methods (TNIPM, IVTCG, PDIP), which confirm that WISTA was sparser than the l₁ regularization methods.

 

Fig. 13 BLT reconstruction results for l₂ regularization methods in a double source case in a 3D view. (a)-(c) Results of TTLS; (d)-(f) Results of Tikhonov.

Download Full Size | PDF

Table 6. Reconstruction results in double sources

View Table | View all tables in this article

 

Fig. 14 BLT reconstruction results of l₁ regularization methods in a double source case in a 3D view. (a)-(c) Results of TNIPM; (d)-(f) Results of IVTCG; (g)-(i) Results of PDIP.

Download Full Size | PDF

 

Fig. 15 BLT reconstruction results of l_p (0 < p < 1) regularization methods in a double source case in a 3D view. (a)-(c) Results of WISTA.

Download Full Size | PDF

3.2. In vivo mouse experiment

The in vivo experiment was presented here to further verify the performance of these algorithms for the application of living animal studies. In the mouse experiment, an implanted source filled with luminescent liquid was placed into the abdomen of the living mouse. The source was 2.7 mm in diameter and 5.1 mm in length. The experimental details and the optical parameters of each organ were presented in [26]. In this experiment, the whole mouse was specified as the PSR, which made the weight matrix A in Eq. (5) undetermined. The reconstructed results are presented in Fig. 16 and Table 7 . The TTLS did not provide results because it could not resolve the undetermined equation. From Fig. 16(c), we found that the reconstruction location of TNIPM was on the bottom left corner of the mouse model and we concluded that its reconstruction failed. From Table 7, we can see that WISTA obtained the best reconstruction results and the run time of PDIP was the longest.

 

Fig. 16 BLT reconstruction results of the in vivo mouse experiment. (a) the 3-D view of the segmented micro-CT slices of the imaged mouse with a luminescent source; (b)-(f) the reconstruction results of Tikhonov, TNIPM, IVTCG, PDIP and WISTA.

Download Full Size | PDF

Table 7. Reconstruction results of the in vivo mouse experiment

View Table | View all tables in this article

4. Discussion and Conclusions

In this paper, we provided a comprehensive assessment of six regularization algorithms developed by our group. In order to investigate the performance of these algorithms in the BLT inverse reconstruction, we carried out five single source experiments, a double source experiment and an in vivo experiment. Because the inverse reconstruction equation is not well posed, the reconstruction results were influenced by various factors, such as the permissible source region, measurement noise, optical properties, tissue specificity and source positions. In this paper, we observed all of the five influential factors which affected the results of BLT reconstruction in the single source experiments. In the double source case, the ability of the six algorithms in resolving double sources was investigated. Finally, the feasibility and practicability of these algorithms were further validated by the in vivo experiment.

The extensive experiments showed that, under a wide range of conditions, there isn’t always a winner that achieves the best performance for all of the investigated cases. Various results for the six reconstruction algorithms under different conditions were obtained and some interesting conclusions could be addressed. Consistent with our anticipation, the reconstruction results became better as the PSR decreased, especially for the l₂ regularization methods (TTLS, Tikhonov), which could be mainly attributed to the fact that the weight matrix was overdetermined and the sparsity characteristics of the inverse problem didn’t pan out. In the noise experiment, four algorithms (Tikhonov, TNIPM, IVTCG, WISTA) were robust to noise. The PDIP became unstable at a noise level of 50% and TTLS became unstable at a noise level of 40%. The optical property variation experiments suggested that we could get better reconstruction results if the mathematical model used to describe the light propagation in the biological tissue fits the optical characteristics of the tissues. The tissue specificity experiments demonstrated the advantage of the heterogeneous model over the homogeneous model, which has been verified by many studies. But for l₂ regularization methods (TTLS, Tikhonov), especially the Tikhonov method, the reconstruction results didn’t improve significantly, so we should use the l_p (0 < p ≤ 1) regularization algorithms in the reconstruction of the heterogeneous model. In the last single source experiments, WISAT was the most robust algorithm at different source positions than the other five algorithms.

In the reconstruction of double sources, TTLS had two cases and Tikhonov had one case where they couldn’t separate the two sources, which contrasted sharply with the reconstruction results of the l_p (0 < p ≤ 1) regularization algorithms (TNIPM, IVTCG, PDIP, WISTA). It is due to the fact that l₂ regularization is easy in generating smooth solutions, which leads to multi-pseudo sources that can’t distinguish the two sources.

In the in vivo mouse experiment, except for the TNIPM method, other algorithms performed normally. As in the double source experiment, WISTA still obtained the best performance because it could maintain its preferable sparsity characteristics. We saw that the reconstruction results of the in vivo mouse experiment were unsatisfactory compared to the simulation experiment. Actually, many factors influenced the reconstruction results of the in vivo experiment, such as the signal acquisition, image registration of the CT data and optical data, the quality of mesh dissection, the accuracy of the optical propagation model, the inverse reconstruction algorithm and so on.

For most experiments, the l_p (0 < p ≤ 1) regularization algorithms (TNIPM, IVTCG, PDIP, WISTA) achieved better reconstruction results than l₂ regularization (TTLS and Tikhonov). In the four l_p (0 < p ≤ 1) regularization algorithms, IVTCG and WISTA were generally more effective than the other two methods.

In order to evaluate the l_p-regularization-based BLT reconstruction algorithms, we limited the review to six algorithms (TTLS, Tikhonov, TNIPM, IVTCG, PDIP, WISTA) that were developed by our group. This is mainly because our group has been engaged in the research of bioluminescence tomography for years, we have the source code of the six algorithms and we can show their best performance and give them a fair comparison. Some interesting conclusions were obtained in this paper which hope to provide useful information for the researcher in selecting and designing algorithms for BLT reconstruction.