Truncated Total Least Squares Method with a Practical Truncation Parameter Choice Scheme for Bioluminescence Tomography Inverse Problem

Abstract

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.

1. Introduction

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 𝑙₂-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.

2. Methodology

2.1. Diffusion Approximation and Boundary Condition

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]: −∇⋅(𝐷(𝐱)∇Φ(𝐱))+𝜇_𝑎(𝐱)Φ(𝐱)=𝑆(𝐱)(𝐱∈Ω), (1) where Ω∈𝑅³ is the bounded domain, Φ(𝐱) represents the photon flux density, and 𝑆(𝐱) denotes the energy density distribution of an internal bioluminescence source, 𝐷(𝐱)=1/(3(𝜇_𝑎(𝐱)+𝜇^′_𝑠(𝐱))) is the optical diffusion coefficient with 𝜇_𝑎(𝐱) being the optical absorption coefficient and 𝜇^′_𝑠(𝐱) 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]: Φ(𝐱)+2𝐴(𝐱;𝑛,𝑛^′)𝐷(𝐱)(𝑣(𝐱)⋅∇Φ(𝐱))=0(𝐱∈𝜕Ω), (2) where 𝜕Ω denotes the boundary, 𝑣(𝐱) represents the unit outer normal on 𝜕Ω,𝐴(𝐱;𝑛,𝑛^′)≈(1+𝑅(𝐱))/(1−𝑅(𝐱)) and 𝑅≈−1.4399𝑛⁻²+0.7099𝑛⁻¹+0.6681+0.0636𝑛, n is the ratio of optical reflective index of the inner tissue to that outside the boundary, and 𝑛^′ 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]:𝑄(𝐱)=−𝐷(𝐱)(𝑣(𝐱)⋅∇Φ(𝐱))=Φ(𝐱)2𝐴(𝐱;𝑛,𝑛^′)(𝐱∈𝜕Ω). (3)

2.2. The Model of BLT Reconstruction

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 𝑆 and the photon flux density Φ can be expressed as the following matrix form: 𝑀Φ=𝐹𝑆, (4) where Φ and 𝑆 are the collection of all the nodal values of the photon flux density and source density, 𝑀=𝐾+𝐶+𝐵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 Φ=𝑀⁻¹𝐹𝑆. In fact, only partial photon on the boundary can be acquired in the BLT experiment, therefore, Φ can be partitioned into the measurable boundary data Φ^𝑚 and other immeasurable values Φ^𝑖, and thus the reconstruction of the bioluminescent source is to identify the unknown vector 𝑆 from the photon flux density Φ^𝑚. 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 Φ^𝑚 and the source energy density distribution 𝑆^𝑝 in the light source permissible region, that is, 𝐴𝑆^𝑝=Φ^𝑚, (5) where the coefficient matrix 𝐴 is ill-conditioned and can cause severe numerical instabilities in the solution. Therefore, it cannot be directly solved using a simple least squares method.

2.3. Regularization

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: min‖𝐴𝑆^𝑝−Φ^𝑚‖²+𝜆‖𝑆^𝑝‖², (6) where 𝜆>0 is a properly chosen regularization parameter. As a function of the regularization parameter, the solution of (5) is given by 𝑆^𝑝=(𝐴^𝑇𝐴+𝜆𝐼)⁻¹𝐴^𝑇Φ^𝑚. (7) However, the reconstructions with Tikhonov regularization method assume that the coefficient matrix 𝐴 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: 𝐴min_,Φ^𝑚‖(𝐴,Φ^𝑚)−(𝐴,Φ^𝑚)‖_𝐹subjectto𝐴𝑆^𝑝=Φ^𝑚, (8) where ‖⋅‖_𝐹 denotes the Frobenius norm, 𝐴 and Φ^𝑚 are the error versions of 𝐴 and Φ^𝑚, respectively, and (𝐴,Φ^𝑚) is the augmented matrix that combines matrix 𝐴 and vector Φ^𝑚 by using Φ^𝑚 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 (𝐴,Φ^𝑚), 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.(1)Compute the SVD of the augmented matrix (𝐴,Φ^𝑚)(𝐴,Φ^𝑚)=𝑈Σ𝑉^𝑇=^𝑛+1_𝑖=1𝑢_𝑖𝜎_𝑖𝑣_𝑖^𝑇, (9) where 𝑈 is m×(n+1), 𝑉 is (𝑛+1)×(𝑛+1), and U^𝑇U=𝑉^𝑇𝑉=𝐼_𝑛+1, and Σ is an (𝑛+1)×(𝑛+1) diagonal matrix with the singular values 𝜎₁>𝜎₂>⋯>𝜎_𝑛+1 on the diagonal.(2)Select a truncation parameter 𝑘≤min(𝑛,rank(𝐴,Φ^𝑚)).(3)Partition the matrix 𝑉∈R^{(𝑛+1)×(𝑛+1)} such that 𝑉𝑉=¹¹𝑉₂₁𝑉₁₂𝑉₂₂,𝑉₁₁∈R^𝑛×𝑘,𝑉₂₂∈R^{1×(𝑛+1−𝑘)}. (10)(4)Then the TTLS solution is given by 𝑆^𝑝_TTLS=−𝑉₁₂𝑉^†₂₂𝑉=−¹²𝑉^𝑇₂₂‖𝑉₂₂‖²₂. (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 𝑘 must be carefully determined.

2.4. Choice of the Truncation Parameter

MGCV criterion proposed by Sima in [17] makes use of the filter factor formulation of the TTLS solution proved in [11]: 𝑆^𝑝_TTLS=^𝑛_𝑖=1𝑓_𝑖𝑢^𝑇_𝑖Φ^𝑚𝜎_𝑖𝑣_𝑖, (12) where the filter factor values 𝑓_𝑖=^𝑘_𝑗=1𝑣²_𝑛+1,𝑗‖𝑉^𝑘₂₂‖²(𝜎²_𝑖𝜎²_𝑖−𝜎²_𝑗),𝑖=1,…,𝑛. (13)

The property used for choosing the truncation parameter 𝑘 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 𝑘, the filter factors with indices 𝑖=1,…,𝑘 stay close to 1 and the filter factors with indices 𝑖=𝑘+1,…,𝑛 stay close to 0; when the truncation level gradually increases to a certain critical value, the filter factors with indices nearby 𝑘 increase dramatically. It implies a way to identify the value of 𝑘 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: 𝐺=‖𝐴𝑆^𝑝_reg−Φ^𝑚‖²(trace(𝐼−𝐴𝐴^†))², (14) where 𝐴^† presents the pseudoinverse of 𝐴. With filter factors, the denominator can be computed by means of the following expression: trace(𝐼−𝐴𝐴^†)=𝑚−(𝑛−𝑝)−^𝑝_𝑖=1𝑓_𝑖, (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_𝑘 as the effective number of parameters: enp_𝑘=^𝑛_𝑖=1𝑓_𝑖=^𝑛𝑘_𝑖=1𝑗=1𝑣²_𝑛+1,𝑗‖𝑉^𝑘₂₂‖²(𝜎²_𝑖𝜎²_𝑖−𝜎²_𝑗). (16)

According to the properties of filter factors mentioned above, for a 𝑘 above a certain critical value, the filter factors for TTLS solutions with indices nearby 𝑘 are larger than 1. A fact can be derived that enp_𝑘 is greater than 𝑘 when 𝑘 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 𝐺=‖𝐴𝑆^𝑝_TTLS−Φ^𝑚‖²(𝑚−enp_𝑘)². (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.

Step 1. Use the MGCV criterion to get an initial truncation parameter 𝑘 and compute 𝑘_max, where 𝑘_max≤𝑛 is the maximum 𝑘 such that enp₁≤enp₂≤⋯≤enp_𝑘max≤𝑚; at the same time an array of GCV function values 𝐺(𝑖) is obtained, 𝑖=[1,𝑘_max].

Step 2. For 𝑖∶𝑘∼𝑘_max, find the local minimum points that satisfy the conditions: 𝐺(𝑖−1)>𝐺(𝑖) and 𝐺(𝑖+1)>𝐺(𝑖).

Step 3. For all the local minimum points, compute the residual error ‖𝐴𝑆^′_TTLS,𝑖−Φ^𝑚‖, where 𝑆^′_TTLS,𝑖 is an approximation of the TTLS solution for a given truncation level 𝑖, that is, only the top 70% of the nodal values are kept for computation convenience, and then the final truncation parameter 𝑘=argmin‖𝐴𝑆^′_TTLS,𝑖−Φ^𝑚‖.

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

3. Experiments and Results

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‖𝐴𝑆^𝑝−Φ^𝑚‖/‖Φ^𝑚‖. Absolute error (AE) of the reconstructed source location is used to describe the accuracy of the reconstruction, which is defined by(𝑥_𝑖,𝑟−𝑥_𝑖)²+(𝑦_𝑖,𝑟−𝑦_𝑖)²+(𝑧_𝑖,𝑟−𝑧_𝑖)², where (𝑥_𝑖,𝑟,𝑦_𝑖,𝑟,𝑧_𝑖,𝑟) is the reconstructed center of each source and (𝑥_𝑖,𝑦_𝑖,𝑧_𝑖) 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 ∫𝑆(𝐱)𝑑𝐱 of the source intensity over its support [25].

3.1. Numerical Simulation Verifications

In the numerical simulation, a 30 mm diameter and 30 mm 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.5 mm diameter with 1 nW/mm³energy density are located in the left lung and the centers are 𝑆₁=(−9mm,−1.5mm,15mm) and 𝑆₂=(−9mm,1.5mm,15mm), respectively. The power of each source is 0.5236 nW. In the following single source case, only the source centered at 𝑆₁ is considered.

Table 1 

Optical parameters of the heterogeneous phantom.

(a)

(b)

Figure 1 

(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={(𝑥,𝑦,𝑧)∶8<(𝑥²+𝑦²)^1/2<12,13.5<𝑧<16.5} [25], where (𝑥,𝑦,𝑧) 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 Φ^𝑚 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 𝐴, 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 𝑘 for TTLS. It should be noticed that the parameter k identified by MGCV is 64, whereas the optimal parameter 𝑘 obtained by IGCV is 78. It is because ‖𝐴𝑆^′_TTLS,78−Φ^𝑚‖ is 0.0038 and ‖𝐴𝑆^′_TTLS,64−Φ^𝑚‖ 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 2 and 3.

Table 2 

Quantitative results in single source case.

Table 3 

Quantitative results in double source case.

(a)

(b)

(c)

Figure 2 

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 selection for TTLS with IGCV.

In single source test, we found that all the methods under consideration can detect the source with the same center location 𝑆^𝑅₁=(−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 𝐴, but TTLS combined with IGCV performs best in all quantitative indices under different noise or error levels.

(a)

(b)

Figure 3 

Reconstructed results under measurement noise level of 10% and Gaussian system error level of 1% with TTLS + IGCV: (a) 𝑥-𝑦 view of at 𝑧=15mm plane, (b) 𝑦-𝑧 view at 𝑥 =−9mm plane; the white circle indicates the real source.

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 𝑆^𝑅₁=(−9.20mm,−1.62mm,14.12mm) and 𝑆^𝑅₂ = (−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 𝑆^𝑅₁ and 𝑆^𝑅₂, 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 𝐴, 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

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 𝑥-𝑦 views at 𝑧=15mm plane of the results by Tikhonov + GCV, TTLS + MGCV, and TTLS + IGCV; (d), (e), and (f) are the corresponding 𝑦-𝑧 views at 𝑥=9.5mm plane of the reconstruction results, respectively; the white circle indicates the real source.

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 10 mm in diameter is utilized, which is expressed as PR^′={(𝑥,𝑦,𝑧)∣((𝑥+7.5)²+𝑦²+(𝑧−15)²)^1/2<5,(𝑥,𝑦,𝑧)∈Leftlung}. 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 𝐴 in this section.

Table 4 

Quantitative results for ball shape permission region𝑃𝑅^′ in double source case.

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 𝑆₁ is up to 1.2 mm. 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.

3.2. Physical Experiment Verifications

A physical experiment was carried out to further investigate the performance of the proposed method. A cylindrical phantom of 45 mm height and 22.5 mm radius was designed to evaluate different methods. The phantom shown in Figure 5(a) was made from nylon, and one small hole of 2.95 mm radius and 21 mm 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.15 mL, thus a cylindrical source with a 2.95 mm radius and 5.4 mm 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 660 nm are 0.91 mm^-1 and 0.0138 mm^-1, respectively.

(a)

(b)

(c)

Figure 5 

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).

(a)

(b)

(c)

(d)

(e)

Figure 6 

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 phantom after mapping from the CCD camera.

According to the photon flux density distribution on the phantom surface, the source permissible region is set as PR^′′={(𝑥,𝑦,𝑧)∶7<(𝑥²+𝑦²)^1/2<13,𝑥>−3,19.7<𝑧<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.76 mm away from the actual source center.

Table 5 

Reconstruction results in physical phantom experiment.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7 

Reconstructed results in physical experiment: (a)–(c) are the 𝑥-𝑦 views at 𝑧=28mm plane of the reconstructed results using Tikhonov + GCV, TTLS + MGCV, and TTLS + IGCV method, respectively; (d)–(f) are the 𝑦-𝑧 views at 𝑥=9.5mm plane of the corresponding results, where the white contours indicates the real sources.

4. Discussion and Conclusion

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 𝐴 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 𝐴. 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 𝐴, 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.

Acknowledgments

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.