Home | About | Journals | Submit | Contact Us | Français |

**|**Int J Biomed Imaging**|**v.2010; 2010**|**PMC2874932

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Methodology
- 3. Experiments and Results
- 4. Discussion and Conclusion
- References

Authors

Related links

Int J Biomed Imaging. 2010; 2010: 291874.

Published online 2010 May 19. doi: 10.1155/2010/291874

PMCID: PMC2874932

Xiaowei He,^{1, 2}^{} Jimin Liang,^{1} Xiaochao Qu,^{1} Heyu Huang,^{1} Yanbin Hou,^{1} and Jie Tian^{1, 3}^{}^{*}

*Jie Tian: Email: gro.eeei@nait

Academic Editor: Guo Wei

Received 2009 September 29; Revised 2010 March 4; Accepted 2010 March 8.

Copyright © 2010 Xiaowei He 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.

This article has been cited by other articles in PMC.

In bioluminescence tomography (BLT), reconstruction of internal bioluminescent source distribution from the surface optical signals is an ill-posed inverse problem. In real BLT experiment, apart from the measurement noise, the system errors caused by geometry mismatch, numerical discretization, and optical modeling approximations are also inevitable, which may lead to large errors in the reconstruction results. Most regularization techniques such as Tikhonov method only consider measurement noise, whereas the influences of system errors have not been investigated. In this paper, the truncated total least squares method (TTLS) is introduced into BLT reconstruction, in which both system errors and measurement noise are taken into account. Based on the modified generalized cross validation (MGCV) criterion and residual error minimization, a practical parameter-choice scheme referred to as improved GCV (IGCV) is proposed for TTLS. Numerical simulations with different noise levels and physical experiments demonstrate the effectiveness and potential of TTLS combined with IGCV for solving the BLT inverse problem.

In recent years, molecular imaging has emerged as a promising tool in basic, preclinical and clinical research for monitoring a variety of molecular and cellular processes in living organisms [1–4]. As one of molecular imaging modality, bioluminescence tomography (BLT) has attracted much attention due to its exquisite sensitivity and cost effectiveness.

The key problem of BLT is to reconstruct the bioluminescent source distribution inside a biological tissue from the optical signals detected on the body surface, which is a highly ill-posed inverse problem. By using numerical method such as finite element method (FEM), the inverse problem of BLT can be formulated into a nonsquare matrix equation, where the coefficient matrix is typically ill-conditioned [5]. Hence overcoming the ill-posedness and seeking a stable solution of the matrix equation are the major issues of BLT inverse problem. For this purpose, the inverse problem is often transformed to a least squares problem incorporated with the regularization technique. Tikhonov regularization is the most widely-used method in BLT reconstruction [6–8]. It is aiming to stabilize the inverse of an ill-conditioned operator and minimize the effects of the inevitable error by minimizing a trade-off between the loss function and the *l*_{2}-norm [8]. However, previous studies based on Tikhonov regularization only consider noise in the measurement. In fact, some system errors also exist in the computed coefficient matrix of the system equation. These errors may take place in such aspects as FEM discretization, geometrical mismatch, optical parameters inaccuracy and model approximation, and so forth. System errors as well as measurement noise are inevitable in real BLT experiment, which may lead to large errors for the reconstruction results.

The total least squares (TLS) method is a generalization of the least squares approximation method when the data in both sides of the matrix equation are perturbed [9, 10]. Based on the TLS method, the truncated total least squares (TTLS) method is proposed for regularization of ill-conditioned linear systems [11]. It is inspired by truncated singular value decomposition (TSVD) which aims at limiting the contribution of noise by cutting off a certain number of terms in the singular value decomposition of coefficient matrix [12]. Truncation level plays the role of regularization parameter in truncation methods, which has great influence on the quality of the solution. As a result, determining an appropriate truncation level for TTLS is a critical step in the inverse procedure. Most existing parameter-choice schemes such as L-curve, discrepancy principle, generalized cross-validation (GCV), and zero crossing methods assume that the coefficient matrix is exactly known, that is, it is not contaminated by noises or errors [13–16]. In [17], a truncation level choice criterion named modified GCV (MGCV) is proposed for TTLS method; theoretical analysis and simulation tests show its potential for solving ill-posed linear system. However, it has been recognized that choice schemes of regularization parameter are mostly problem-dependent and practical parameter-choice scheme for BLT reconstruction deserves further study.

In this paper, the aim of our study is to extend the BLT reconstruction to the case including both the measurement noise and the system errors. For this purpose, TTLS method combined with a practical scheme termed as improved GCV (IGCV) is proposed to solve the BLT inverse problem. In the next section, our methodology of solving the inverse problem in BLT is described. In Section 3, we demonstrate the performance of the TTLS method combined with IGCV scheme in BLT reconstruction using numerical simulation and physical experiments in various source and noise level settings. Finally, we draw a conclusion and discuss the relevant issues.

In general, light propagation in living subjects is mainly hindered by both tissue scattering and absorption [7, 8]. Considering that bioluminescent photons belong to the near-infrared region where scattering predominates over absorption [3], the propagation of photon can be well modeled by the following steady-state diffusion equation [18]:

$$-\nabla \xb7\left(D\left(\mathbf{x}\right)\nabla \mathrm{\Phi}\left(\mathbf{x}\right)\right)+{\mu}_{a}\left(\mathbf{x}\right)\mathrm{\Phi}\left(\mathbf{x}\right)=S\left(\mathbf{x}\right)\hspace{1em}\left(\mathbf{x}\in \Omega \right),$$

(1)

where *Ω* *R*^{3} is the bounded domain, Φ(**x**) represents the photon flux density, and *S*(**x**) denotes the energy density distribution of an internal bioluminescence source, *D*(**x**) = 1/(3(*μ*_{a}(**x**) + *μ*_{s}′(**x**))) is the optical diffusion coefficient with *μ*_{a}(**x**) being the optical absorption coefficient and *μ*_{s}′(**x**) the reduced scattering coefficient, respectively.

Assuming that the BLT experiment is performed in a totally dark environment, the equation is subject to a Robin boundary condition [18]:

$$\mathrm{\Phi}\left(\mathbf{x}\right)+2A\left(\mathbf{x};n,{n}^{\prime}\right)D\left(\mathbf{x}\right)\left(v\left(\mathbf{x}\right)\xb7\nabla \mathrm{\Phi}\left(\mathbf{x}\right)\right)=0\hspace{1em}\left(\mathbf{x}\in \partial \Omega \right),$$

(2)

where *Ω* denotes the boundary, *v*(**x**) represents the unit outer normal on *Ω*, *A*(**x**; *n*, *n*′) ≈ (1 + *R*(**x**))/(1 − *R*(**x**)) and *R* ≈ −1.4399*n*^{−2} + 0.7099*n*^{−1} + 0.6681 + 0.0636*n*, *n *is the ratio of optical reflective index of the inner tissue to that outside the boundary, and *n*′ is close to 1.0 when the subject is in air [8]. In a bioluminescent imaging experiment, the measurable photon flux density on *Ω* can be calculated by the following outgoing radiation [18]:

$$Q\left(\mathbf{x}\right)=-D\left(\mathbf{x}\right)\left(v\left(\mathbf{x}\right)\xb7\nabla \mathrm{\Phi}\left(\mathbf{x}\right)\right)=\frac{\mathrm{\Phi}\left(\mathbf{x}\right)}{2A\left(\mathbf{x};n,{n}^{\prime}\right)}\hspace{1em}\left(\mathbf{x}\in \partial \Omega \right).$$

(3)

Based on (1), (2), and (3), the essence of the BLT reconstruction is to estimate the light source distribution inside the biological tissues from the measured flux on the surface, given the corresponding optical parameters of the tissues. In order to solve the BLT inverse problem, FEM was introduced to solve the diffusion equation in [8, 18–20] because of its capability to process volume with arbitrary geometries. After the discretization using FEM, the linear relation between the bioluminescence source intensity *S* and the photon flux density Φ can be expressed as the following matrix form:

$$M\mathrm{\Phi}=FS,$$

(4)

where Φ and *S* are the collection of all the nodal values of the photon flux density and source density, *M* = *K* + *C* + *B*is a positive-definite matrix, and *K*, *C,* and *B* are called the mass, stiff, and boundary matrix, respectively. The photon density Φ can be obtained from Φ = *M*^{−1}*F**S*. In fact, only partial photon on the boundary can be acquired in the BLT experiment, therefore, Φ can be partitioned into the measurable boundary data Φ^{m} and other immeasurable values Φ^{i}, and thus the reconstruction of the bioluminescent source is to identify the unknown vector *S* from the photon flux density Φ^{m}. According to the uniqueness theorem, the BLT solution is not unique in the general case [21]. Some prior information or constraints such as permissible area of source should be imposed on the unknown variables to obtain a meaningful reconstruction result. Considering the source permissible region, we can obtain the linear relation between the photon flux density Φ^{m} and the source energy density distribution *S*^{p} in the light source permissible region, that is,

$$A{S}^{p}={\mathrm{\Phi}}^{m},$$

(5)

where the coefficient matrix *A* is ill-conditioned and can cause severe numerical instabilities in the solution. Therefore, it cannot be directly solved using a simple least squares method.

In order to obtain a stable solution, regularization methods are typically used for solving inverse problems [8, 22, 23]. The commonly used Tikhonov regularization method approximately solves (5) by converting it into the following minimization problem:

$$\mathrm{min}\hspace{0.17em}{||A{S}^{p}-{\mathrm{\Phi}}^{m}||}^{2}+\lambda {||{S}^{p}||}^{2},$$

(6)

where *λ* > 0 is a properly chosen regularization parameter. As a function of the regularization parameter, the solution of (5) is given by

$${S}^{p}={\left({A}^{T}A+\lambda I\right)}^{-1}{A}^{T}{\mathrm{\Phi}}^{m}.$$

(7)

However, the reconstructions with Tikhonov regularization method assume that the coefficient matrix *A* is exactly known and noises only exist in the measurement. The regularization solutions computed by (7) do not take system errors into account.

As mentioned in the introduction, TLS method is designed for the case that both sides of the matrix equation are subject to errors. BLT inverse problem can be stated with TLS formulation as follows:

$$\underset{\tilde{A},{\tilde{\mathrm{\Phi}}}^{m}}{\mathrm{min}\hspace{0.17em}}{\hspace{0.17em}||\left(A,{\mathrm{\Phi}}^{m}\right)-\left(\tilde{A},{\tilde{\mathrm{\Phi}}}^{m}\right)||}_{F}\mathrm{\hspace{0.17em}\hspace{0.17em}}\hspace{1em}\text{subject}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{to}\mathrm{\hspace{0.17em}\hspace{0.17em}}\tilde{A}{S}^{p}={\tilde{\mathrm{\Phi}}}^{m},$$

(8)

where ||·||_{F} denotes the Frobenius norm, $\tilde{A}$ and ${\tilde{\Phi}}^{m}$ are the error versions of *A* and Φ^{m}, respectively, and (*A*, Φ^{m}) is the augmented matrix that combines matrix *A* and vector Φ^{m} by using Φ^{m} as the last column of the new matrix. Based on the TLS method, TTLS method is proposed by Fierro in [11] for regularization of ill-conditioned linear systems. In TTLS, the redundant information in (*A*, Φ^{m}), associated with the small singular values, is discarded and the original ill-conditioned problem is replaced with another appropriate and more well-conditioned problem.

The TTLS algorithm used in this paper can be summarized as follows.

- Compute the SVD of the augmented matrix (
*A*, Φ^{m})where $\tilde{U}$ is$$\left(A,{\mathrm{\Phi}}^{m}\right)=\tilde{U}\tilde{\mathrm{\Sigma}}{\tilde{V}}^{T}={\displaystyle \sum}_{i=1}^{n+1}{\stackrel{\u0305}{u}}_{i}{\stackrel{\u0305}{\sigma}}_{i}{{\stackrel{\u0305}{v}}_{i}}^{T},$$(9)*m*×(*n*+1), $\tilde{V}$ is (*n*+ 1) × (*n*+ 1), and ${\tilde{\text{U}}}^{T}\tilde{\text{U}}={\tilde{V}}^{T}\tilde{V}={I}_{n+1},$ and $\tilde{\Sigma}$ is an (*n*+ 1) × (*n*+ 1) diagonal matrix with the singular values ${\stackrel{\u0305}{\sigma}}_{1}>{\stackrel{\u0305}{\sigma}}_{2}>\cdots >{\stackrel{\u0305}{\sigma}}_{n+1}$ on the diagonal. - Select a truncation parameter
*k*≤ min(*n*, rank(*A*, Φ^{m})). - Partition the matrix $\tilde{V}\in {\text{R}}^{(n+1)\times (n+1)}$ such that$$\tilde{V}=\left(\begin{array}{cc}\begin{array}{c}{\tilde{V}}_{11}\\ {\tilde{V}}_{21}\end{array}& \begin{array}{c}{\tilde{V}}_{12}\\ {\tilde{V}}_{22}\end{array}\end{array}\right),\hspace{1em}{\tilde{V}}_{11}\in {\text{R}}^{n\times k},\mathrm{\hspace{0.17em}\hspace{0.17em}}{\tilde{V}}_{22}\in {\text{R}}^{1\times (n+1-k)}.$$(10)
- Then the TTLS solution is given by$${S}_{\text{TTLS}}^{p}=-{\tilde{V}}_{12}{\tilde{V}}_{22}^{\u2020}=-\frac{{\tilde{V}}_{12}{\tilde{V}}_{22}^{T}}{{||{\tilde{V}}_{22}||}_{2}^{2}}.$$(11)

In fact, the aim of TTLS regularization is to appropriately identify an optimal truncation level, and then to construct a truncated solution that can capture the essential features of the unknown true solution, without explicit knowledge about the true solution and even without a priori knowledge about the magnitude of the noise in the data. For this purpose, truncation level *k* must be carefully determined.

MGCV criterion proposed by Sima in [17] makes use of the filter factor formulation of the TTLS solution proved in [11]:

$${S}_{\text{TTLS}}^{p}={\displaystyle \sum}_{i=1}^{n}{f}_{i}\frac{{\stackrel{\u0305}{u}}_{i}^{T}{\mathrm{\Phi}}^{m}}{{\stackrel{\u0305}{\sigma}}_{i}}{\stackrel{\u0305}{v}}_{i},$$

(12)

where the filter factor values

$${f}_{i}={\displaystyle \sum}_{j=1}^{k}\frac{{v}_{n+1,j}^{2}}{{||{V}_{22}^{k}||}^{2}}\left(\frac{{\stackrel{\u0305}{\sigma}}_{i}^{2}}{{\stackrel{\u0305}{\sigma}}_{i}^{2}-{\sigma}_{j}^{2}}\right),\hspace{1em}i=1,\dots ,n.$$

(13)

The property used for choosing the truncation parameter *k* is that when the parameter is greater than a certain crucial value, the TTLS solution is very sensitive to the noise or errors. Specifically, for small truncation level *k*, the filter factors with indices *i* = 1,…, *k* stay close to 1 and the filter factors with indices *i* = *k* + 1,…, *n* stay close to 0; when the truncation level gradually increases to a certain critical value, the filter factors with indices nearby *k* increase dramatically. It implies a way to identify the value of *k* where the filter factors change their steady behavior into erratic growth behavior.

As for the regularization problem in BLT, the choice of regularization parameter with classical GCV is by means of minimizing the GCV function:

$$G=\frac{{||A{S}_{\text{reg}}^{p}-{\mathrm{\Phi}}^{m}||}^{2}}{{\left(\text{trace}\left(I-A{A}^{\u2020}\right)\right)}^{2}},$$

(14)

where *A*^{†} presents the pseudoinverse of *A*. With filter factors, the denominator can be computed by means of the following expression:

$$\text{trace}\left(I-A{A}^{\u2020}\right)=m-\left(n-p\right)-{\displaystyle \sum}_{i=1}^{p}{f}_{i},$$

(15)

where *p* is the rank of matrix ∑ with the singular values on the diagonal. We denote the sum of the filter factors of TTLS solution by enp_{k} as the effective number of parameters:

$$\text{en}{\text{p}}_{k}={\displaystyle \sum}_{i=1}^{n}{f}_{i}={\displaystyle \sum}_{i=1}^{n}\hspace{0.17em}{\displaystyle \sum}_{j=1}^{k}\frac{{v}_{n+1,j}^{2}}{{||{V}_{22}^{k}||}^{2}}\left(\frac{{\stackrel{\u0305}{\sigma}}_{i}^{2}}{{\stackrel{\u0305}{\sigma}}_{i}^{2}-{\sigma}_{j}^{2}}\right).$$

(16)

According to the properties of filter factors mentioned above, for a *k* above a certain critical value, the filter factors for TTLS solutions with indices nearby *k* are larger than 1. A fact can be derived that enp_{k} is greater than *k* when *k* reaches this critical value, which is used to modify the above classic GCV function to suit the TTLS case. And then the MGCV criterion for TTLS is obtained

$$G=\frac{{||A{S}_{\text{TTLS}}^{p}-{\mathrm{\Phi}}^{m}||}^{2}}{{\left(m-\text{en}{\text{p}}_{k}\right)}^{2}}.$$

(17)

However, the regularization parameter directly identified by (17) may be not optimal for the specific BLT reconstruction problem. Inspired by L-curve method, we propose a hybrid scheme that combines MGCV with the minimization of the corresponding residual norm for regularization parameter choice. The IGCV scheme for TTLS is summarized in the following steps.

Use the MGCV criterion to get an initial truncation parameter *k* and compute *k*_{max}, where *k*_{max} ≤ *n* is the maximum *k* such that enp_{1} ≤ enp_{2} ≤ ≤enp_{kmax} ≤ *m*; at the same time an array of GCV function values *G*(*i*) is obtained, *i* = [1, *k*_{max}].

For *i* : *k* ~ *k*_{max}, find the local minimum points that satisfy the conditions: *G*(*i* − 1) > *G*(*i*) and *G*(*i* + 1) > *G*(*i*).

For all the local minimum points, compute the residual error ||*A**S*_{TTLS,i}′ − Φ^{m}||, where *S*_{TTLS,i}′ is an approximation of the TTLS solution for a given truncation level *i*, that is, only the top 70% of the nodal values are kept for computation convenience, and then the final truncation parameter *k* = argmin ||*A**S*_{TTLS,i}′ − Φ^{m}||.

Thus, a proper truncation parameter *k* for TTLS is sought according the above IGCV scheme.

The experiments implemented in this section are to test the performance of TTLS combined with IGCV for BLT inverse problem. To demonstrate the effectiveness of the proposed scheme, we compare the following reconstruction algorithms: Tikhonov method with classical GCV (Tik-GCV), TTLS method with MGCV (TTLS-M), and TTLS method with the proposed IGCV (TTLS-I). The parameter-choice scheme of Tikhonov method is different from that of TTLS method because MGCV and IGCV are specially designed for TTLS. A similar scheme, namely, classical GCV, is adopted in Tikhonov method for comparison convenience. The qualities of the reconstruction are assessed by the following quantitative indices: relative residual error (RRE), reconstructed location error, and reconstructed source power. Here, RRE is used to depict the extent of the solution fitting the measured data and is defined as ||*A**S*^{p} − Φ^{m}||/||Φ^{m}||. Absolute error (AE) of the reconstructed source location is used to describe the accuracy of the reconstruction, which is defined by$\sqrt{{({x}_{i,r}-{x}_{i})}^{2}+{({y}_{i,r}-{y}_{i})}^{2}+{({z}_{i,r}-{z}_{i})}^{2}}$
, where (*x*_{i,r}, *y*_{i,r}, *z*_{i,r}) is the reconstructed center of each source and (*x*_{i}, *y*_{i}, *z*_{i}) the actual center. Considering the ill-posedness of the BLT inverse problem, it is difficult to discriminate the influence of small source of high density and large one of low density [24]. So we prefer reconstructed source power compared with the actual value to source density for evaluating the quality of the reconstruction results. And the source power is estimated by computing the integral ∫*S*(**x**)*d ***x** of the source intensity over its support [25].

In the numerical simulation, a 30mm diameter and 30mm high cylindrical mouse chest phantom is designed to evaluate the performance of the reconstruction method. The structure of the phantom is shown in Figure 1(a). The phantom is heterogeneous and the corresponding optical parameters are set as in Table 1 [25]. Two sphere sources of 0.5mm diameter with 1nW/mm^{3}energy density are located in the left lung and the centers are *S*_{1} = (−9mm, −1.5mm, 15mm) and *S*_{2} = (−9mm, 1.5mm, 15mm), respectively. The power of each source is 0.5236nW. In the following single source case, only the source centered at *S*_{1} is considered.

(a) A cross-section through two luminescent sources (S) in the left lung of a mouse phantom consisting of bone (B), heart (H), lungs (L), and tissue (T). (b) A 3D view of the permissible region.

In order to reduce the ill-posedness of the inverse problem, a priori information of the source permissible region (PR) is incorporated to our method, which is shown in Figure 1(b) as PR = {(*x*, *y*, *z*) : 8 < (*x*^{2} + *y*^{2})^{1/2} < 12, 13.5 < *z* < 16.5} [25], where (*x*, *y*, *z*) is the coordinates of the corresponding FEM mesh vertices.

Generally speaking, simulated data used in reconstruction algorithms for inverse problems often come from the numerical solution of the forward problem. To avoid the typical issue of *inverse crime*, we use different FEM discretization for the forward process and reconstruction algorithms. Specifically, the forward model contains 11997 mesh vertices corresponding to 66334 tetrahedral elements, whereas the reconstruction model is consisting of 5277 vertices and 27465 tetrahedral elements. In addition, we employ Lagrange-Quadratic interpolation function in the forward process owing to the observation that high-order interpolation function can improve the numerical accuracy of the forward solution [26, 27].

To comprehensively simulate the noise and system errors involved in real BLT experiment, the photon flux density Φ^{m} is added with Gaussian white noise, and the coefficient matrix A is added with a system errors matrix. Due to the complexity of error sources, it is difficult to have an exact mathematic model to describe the system errors accurately. Hence we adopted the commonly used Gaussian white noise [28–30] and exponential noise to simulate the errors in matrix *A*, respectively.

As discussed in Section 2, regularization parameter is the crucial factor that affects the quality of regularization solution to inverse problem. Figure 2 illustrates the determination of regularization parameters in single source case with measurement noise level of 10% and Gaussian system error level of 1%. Among them, Figure 2(c) shows the residual error values of all the local minimum points described in our improved scheme IGCV, which are used for the selection of an optimal truncation parameter *k* for TTLS. It should be noticed that the parameter* k* identified by MGCV is 64, whereas the optimal parameter *k* obtained by IGCV is 78. It is because ||*A**S*_{TTLS,78}′ − Φ^{m}|| is 0.0038 and ||*A**S*_{TTLS,64}′ − Φ^{m}|| is 0.0046, which indicate that 64 is not the optimal parameter value according to IGCV criterion. The determination of regularization parameter in double sources case is similar to that of single source case. For space limitation, we just provide the final regularization parameter obtained in various noise settings in Tables Tables2 and2 and and33.

Regularization parameter determination in single-source case under measurement noise level of 10% and Gaussian system error level of 1%: (a) GCV function curve for Tikhonov, (b) MGCV function curve for TTLS, (c) illustration of the truncation parameter **...**

In single source test, we found that all the methods under consideration can detect the source with the same center location *S*_{1}^{R} = (−9.20mm, −1.62mm, 14.12mm) in different noise levels, but the reconstructed source power varies with different reconstruction methods. Although the absolute error of the source location is 0.911mm, the reconstructed source center is the nearest node to the original location in the aforementioned FEM discretization. Figure 3 only shows the reconstruction results by our proposed method with measurement noise level of 10% and Gaussian system error level of 1%. The detailed quantitative reconstruction results for the single-source model in various noise settings are listed in Table 2. The optimal results are listed in bold. Based on the simulation results in single-source case, it is clear that all the reconstruction methods can estimate the source location with no matter Gaussian or exponential noise in matrix *A*, but TTLS combined with IGCV performs best in all quantitative indices under different noise or error levels.

Reconstructed results under measurement noise level of 10% and Gaussian system error level of 1% with TTLS + IGCV: (a) *x*–*y* view of at *z* = 15mm plane, (b) *y*–*z* view at *x* =−9mm plane; the white circle indicates the **...**

In the double sources case, both of the two sphere sources located in the left lung are tested. The final reconstruction results are listed in Table 3. Under all the noise conditions considered in this paper, the three methods can reconstruct the two sources at *S*_{1}^{R} = (−9.20mm, −1.62mm, 14.12mm) and *S*_{2}^{R} = (−9.42mm, 1.69mm, 14.94mm), which are 0.911mm and 0.467mm away from the actual ones, respectively. In fact, they are the nearest nodes to the original source locations under the FEM mesh used in our tests. However, with the increase of noise or error level, besides the optimal nodes *S*_{1}^{R} and *S*_{2}^{R}, some artifacts appear in the reconstruction results, which are illustrated in Figure 4. Simulation results in double sources case further show that although there are differences between the results of different noise pattern in matrix *A*, similar conclusions can be obtained. As shown in Table 3, the reconstruction results of TTLS combined with MGCV are comparable to that of TTLS combined with IGCV when noise level is low; whereas with the increase of noise or error, TTLS combined with IGCV outperforms the other methods in all quantitative indices.

Reconstructed results in double source case under measurement noise level of 20% and system error level of 5%. (a), (b), and (c) separately show the *x*–*y* views at *z* = 15mm plane of the results by Tikhonov + GCV, TTLS + MGCV, and TTLS + **...**

For BLT inverse problem, permission region is an effective way to regularize the solution by restricting the source distribution within a proper permissible region. In order to further test the proposed method, a ball shape permissible region of 10mm in diameter is utilized, which is expressed as PR′ = {(*x*, *y*, *z*)((*x* + 7.5)^{2} + *y*^{2} + (*z*−15)^{2})^{1/2} < 5, (*x*, *y*, *z*) Left lung}. The sources settings in this section are the same as the aforementioned double sources case. The source distribution in the ball permission region was reconstructed, and the results are summarized in Table 4. Considering that the different system error pattern has little effect on the reconstruction results in the foregoing simulations, we only add Gaussian noise to the system matrix *A* in this section.

It is shown in Table 4 that TTLS combined with IGCV still performs best under all the noise levels in terms of RRE, reconstructed power and source location. Compared with the results in Table 3, the location accuracy for ball shape permission region PR′ is lower. For example, the largest deviation of the reconstructed position of *S*_{1} is up to 1.2mm. It is clear that all reconstruction methods under consideration suffer from performance degradation with the relaxation of the permission region. However; the proposed method outperforms the other two methods and produces acceptable reconstruction results in our tests.

A physical experiment was carried out to further investigate the performance of the proposed method. A cylindrical phantom of 45mm height and 22.5mm radius was designed to evaluate different methods. The phantom shown in Figure 5(a) was made from nylon, and one small hole of 2.95mm radius and 21mm depth was drilled in the phantom to inject luminescent mixed solution used as the light source. In our physical experiment, the total volume of the mixed solution injected into the hole is 0.15mL, thus a cylindrical source with a 2.95mm radius and 5.4mm height is centered at (9.88mm, 1.5mm, 26.7mm), as shown in Figure 5(b). The optical parameters of the phantom were determined by a time-correlated single photon counting (TCSPC) system specifically constructed for the optical properties of the turbid medium [31]. The measured values of absorption and reduced scattering coefficients at the wavelength around 660nm are 0.91mm^{−1} and 0.0138mm^{−1}, respectively.

Physical phantom. (a) The homogeneous physical phantom; (b) The location of the single source in the phantom; (c) The cross-section of the phantom and the four directions of the CCD camera during data acquisition.

A scientific cooled back-illuminated CCD camera (PIXIS 2048B) is used to collect the outgoing photons from the phantom surface. The photon flux density from different angles can be acquired by rotating the stage under the phantom, as illustrated in Figure 5(c). Figures 6(a)–6(d) exhibits the four views of the cylindrical phantom obtained by the CCD camera, respectively. Because the data captured by CCD camera is planar, mapping it onto 3D surface of the cylindrical phantom must be accomplished before reconstruction, which will also bring some inevitable errors to the measured data [32]. The mapping result was shown in Figure 6(e).

The normalized surface measurement of the homogeneous phantom. (a), (b), (c), and (d) are left view, front view, right view, and back view of the cylindrical phantom on the CCD camera, respectively; (e) is the flux density on the surface of the cylindrical **...**

According to the photon flux density distribution on the phantom surface, the source permissible region is set as PR′′ = {(*x*, *y*, *z*) : 7 < (*x*^{2}+*y*^{2})^{1/2} < 13, *x* > −3, 19.7 < *z* < 33.7}. In the reconstruction process, the phantom model consists of 2734 vertices corresponding to13551 tetrahedral elements. The schemes for the selection of regularization parameters are identical to those in numerical simulations. The final reconstruction results and the corresponding regularization parameter are listed in Table 5. The 3D views of the reconstructed results using different methods are presented in Figure 7, which verified the feasibility and effectiveness of the proposed method. As is evident from the images in Figure 7 and the data in Table 5, TTLS combined with IGCV successfully reconstructed the luminescent source with the minimum distance of 1.76mm away from the actual source center.

Reconstructed results in physical experiment: (a)–(c) are the *x*–*y* views at *z* = 28mm plane of the reconstructed results using Tikhonov + GCV, TTLS + MGCV, and TTLS + IGCV method, respectively; (d)–(f) are the *y*– **...**

BLT reconstruction is a highly ill-posed inverse problem where small measurement noise and system errors in the input data can produce large changes in the results. In addition, bioluminescence signals are generally very weak, thus the noise or errors will significantly affect the reconstruction quality. Regularization technique has played an important role in solving BLT inverse problem. And most of the previous works assume that there is only measurement noise, which affects the right-hand side of the system equations. However, the computed coefficient matrix *A* in the model also has some errors, which may be caused by the calculation errors, the geometrical approximation, optical parameter inaccuracy, as well as the assumption of diffusion equation model itself. For example, the FEM discretization typically adds some errors to the matrix *A*. Hence, there is a need for seeking methods that can deal with the errors in both sides of the system equation. TTLS is a truncation regularization method that can take account of both system errors and measurement noise in the reconstruction process. This method depends on a parameter called truncation level; this single parameter has a significant influence on the regularization solutions. In this paper, IGCV, a practical scheme for determining the truncation parameter, is proposed to be combined with TTLS method for solving BLT inverse problem.

Simulations considering both system errors and measurement noise are conducted to investigate the performance of the proposed reconstruction method. Due to the lack of an accurate model to describe the system errors arising from multiple sources, commonly used Gaussian white noise and exponential noise are adopted to simulate the errors in matrix *A*, respectively. In addition, physical phantom experiments further test the proposed method.

Both the numerical simulations and physical experiments demonstrate the effectiveness of the proposed method. Tests with different noise levels show that TTLS with combined IGCV is able to produce much better reconstruction results than Tikhonov method, and TTLS combined with IGCV performs better than TTLS combined with MGCV, especially when both sides of the system equation are contaminated by measurement noise and system errors. Based on the experiments in this paper, we can draw a preliminary conclusion that TTLS combined with IGCV criterion is a potential reconstruction method for BLT inverse problem. Further investigation of the performance of the proposed method on animal experiments will be conducted in our future work.

This work is supported by the Program of the National Basic Research and Development Program of China (973) under Grant no. 2006CB705700, the Cheung Kong Scholars and Innovative Research Team in University (PCSIRT) under Grant no. IRT0645, the Chair Professors of Cheung Kong Scholars Program of Ministry of Education of China, CAS Hundred Talents Program, the National Natural Science Foundation of China under Grant no. 30873462, 60532050, 30900334, the Beijing Municipal Natural Science Foundation of China under Grants no. 4071003, the CAS Scientific Research Equipment Develop Program (YZ0642, YZ200766), the Natural Science Basic Research Plan in Shaanxi Province of China under Grant no. 2009JQ8018, and the Science Foundation of Northwest University under Grant no. 09NW34.

1. Willmann JK, van Bruggen N, Dinkelborg LM, Gambhir SS. Molecular imaging in drug development. *Nature Reviews Drug Discovery*. 2008;7:591–607. [PubMed]

2. Ntziachristos V, Ripoll J, Wang LV, Weissleder R. Looking and listening to light: the evolution of whole-body photonic imaging. *Nature Biotechnology*. 2005;23(3):313–320. [PubMed]

3. Tian J, Bai J, Yan X-P, Bao S, Li Y, Liang W, Yang X. Multimodality molecular imaging: improving image quality. *IEEE Engineering in Medicine and Biology Magazine*. 2008;27(5):48–57. [PubMed]

4. Contag CH, Bachmann MH. Advances in in vivo bioluminescence imaging of gene expression. *Annual Review of Biomedical Engineering*. 2002;4:235–260. [PubMed]

5. Schweiger M, Arridge SR, Hiraoka M, Delpy DT. The finite element method for the propagation of light in scattering media: boundary and source conditions. *Medical Physics*. 1995;22(11):1779–1792. [PubMed]

6. Tikhonov AN, Aresenin VY. *Solutions of Ill-Posed Problems*. Washington, DC, USA: V. H. Winston & Sons; 1977.

7. Wang G, Cong W, Durairaj K, et al. In vivo mouse studies with bioluminescence tomography. *Optics Express*. 2006;14(17):7801–7809. [PubMed]

8. Lu Y, Zhang X, Douraghy A, et al. Source reconstruction for spectrally-resolved bioluminescence tomography with sparse A priori information. *Optics Express*. 2009;17(10):8062–8080. [PMC free article] [PubMed]

9. Huffel SV, Vandewalle J. *The Total Least Squares Problem: Computational Aspects and Analysis*. Philadelphia, Pa, USA: Society for Industrial Mathematics; 1991. (Frontiers in Applied Mathematics).

10. Golub GH, Van Loan CF. An analysis of the total least squares problem. *SIAM Journal on Numerical Analysis*. 1997;17(6):883–893.

11. Fierro RD, Golub GH, Hansen PC, O’Leary DP. Regularization by truncated total least squares. *SIAM Journal of Scientific Computing*. 1997;18(4):1223–1241.

12. Hansen PC. Truncated singular value decomposition solutions to discrete ill-posed problems with ill-determined numerical rank. *SIAM Journal on Scientific and Statistical Computing*. 1990;11(3):503–518.

13. Hansen PC. Analysis of discrete ill-posed problems by means of the L-curve. *SIAM Review*. 1992;34(4):561–580.

14. Hansen PC, Kilmer ME, Kjeldsen RH. Exploiting residual information in the parameter choice for discrete ill-posed problems. *BIT Numerical Mathematics*. 2006;46(1):41–59.

15. Golub GH, Heath M, Wahba G. Generalized cross-validation as a method for choosing a good ridge parameter. *Technometrics*. 1979;21:215–223.

16. Moody JE. Moody. *Advances in Neural Information Processing Systems*. Vol. 4. Palo Alto, Calif, USA: Morgan Kaufmann; 1992. The effective number of parameters: an analysis of generalization and regularization in nonlinear learning systems; pp. 847–854.

17. Sima DM, Huffel SV. Level choice in truncated total least squares. *Computational Statistics and Data Analysis*. 2007;52(2):1103–1118.

18. Cong W, Wang G, Kumar D, et al. Practical reconstruction method for bioluminescence tomography. *Optics Express*. 2005;13(18):6756–6771. [PubMed]

19. Song X, Wang D, Chen N, Bai J, Wang H. Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm. *Optics Express*. 2007;15(26):18300–18317. [PubMed]

20. Song X, Yi J, Bai J. A parallel reconstruction scheme in fluorescence tomography basedon contrast of independent inversed absorption properties. *International Journal of Biomedical Imaging*. 2006;2006(26):7 pages. Article ID 70839. [PMC free article] [PubMed]

21. Wang G, Li Y, Jiang M. Uniqueness theorems in bioluminescence tomography. *Medical Physics*. 2004;31(8):2289–2299. [PubMed]

22. Cong AX, Wang G. Multispectral bioluminescence tomography: methodology and simulation. *International Journal of Biomedical Imaging*. 2006;2006:7 pages. Article ID 57614.

23. He X, Tian J, Wu Y, Hou Y, Ren N, Penga K. Study of four regularization methods for the inverse problem in bioluminescence tomography. In: Medical Imaging 2009: Biomedical Applications in Molecular, Structural, and Functional Imaging, vol. 7262; 2009; Lake Buena Vista, Fla, USA. p. 10 pages. *Proceedings of SPIE*. 7262B.

24. Dehghani H, Davis SC, Pogue BW. Spectrally resolved bioluminescence tomography using the reciprocity approach. *Medical Physics*. 2008;35(11):4863–4871. [PubMed]

25. Jiang M, Zhou T, Cheng J, Cong W, Wang G. Image reconstruction for bioluminescence tomography from partial measurement. *Optics Express*. 2007;15(18):11095–11116. [PubMed]

26. Hou Y, Tian J, Wu Y, Liang J. A new numerical method for BLT forward problem based on high-order finite elements. *Communications in Numerical Methods in Engineering*. 2009;25(6):667–681.

27. Han R, Liang J, Qu X, et al. A source reconstruction algorithm based on adaptive hp-FEM for bioluminescence tomography. *Optics Express*. 2009;15(18):11095–11116. [PubMed]

28. Zhu W, Wang Y, Yao Y, Chang J, Graber HL, Harbour RL. Iterative total least-squares image reconstruction algorithm for optical tomography by the conjugate gradient method. *Journal of the Optical Society of America A*. 1997;14(4):799–807. [PubMed]

29. Zhu W, Wang Y, Zhang J. Total least-squares reconstruction with wavelets for optical tomography. *Journal of the Optical Society of America A*. 1998;15(10):2639–2650. [PubMed]

30. Shou G, Xia L, Jiang M, Wei Q, Liu F, Crozier S. Truncated total least squares: a new regularization method for the solution of ECG inverse problems. *IEEE Transactions on Biomedical Engineering*. 2008;55(4):1327–1335. [PubMed]

31. Qin D, Zhao H, Tanikawa Y, Gao F. Experimental determination of optical properties in turbid medium by TCSPC technique. In: Progress in Biomedical Optics and Imaging, vol. 6434; 2007; San Jose, Calif, USA. p. 10 pages. *Proceedings of SPIE*. 64342E.

32. Ripoll J, Yessayan D, Zacharakis G, Ntziachristos V. Experimental determination of photon propagation in highly absorbing and scattering media. *Journal of the Optical Society of America A*. 2005;22(3):546–551. [PubMed]

Articles from International Journal of Biomedical Imaging are provided here courtesy of **Hindawi Publishing Corporation**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |