1. Introduction

In vivo small animal optical imaging has become an important tool of biological discovery and preclinical applications [1][2][3]. When mouse models are labeled using optical molecular probes, the probes acting as light sources, reflect corresponding biological information through the emission of visible or near infrared (NIR) light photons. Optical molecular imaging equipment is used to detect the photon distribution over the surface of the small animal to non-invasively investigate these models [4]. In recent years, planar optical molecular imaging, and more specifically bioluminescence imaging, has been extensively applied in tumorigenesis studies, cancer diagnosis, metastasis detection, drug discovery and development, and gene therapies given its convenience and ease operation [5][6]. The technology that is capable to acquire three dimensional information of the light sources will become a next generation instrument for optical molecular imaging. Bioluminescence tomography (BLT) is one such instrument being developed for this purpose [7].

An indispensable parameter for bioluminescence tomography is a priori information, which can be used to localize the optical sources. Theoretically, the source uniqueness proof gives us an important reference [8]. Practically, the richer the a priori information we apply, the further improvements BLT reconstruction can yield. Currently, three types of a priori information are verified and extensively applied in reconstruction algorithms. These include anatomical information [9][10], spectrally-resolved measurements [11][12][9][13], and the distribution of surface photons [14]. Anatomical information is used to assign relevant optical properties to organs. Spectrally resolved data considers the source spectrum and the tissue absorption and scattering characteristics. The use of these a priori information significantly improves source reconstruction. The temperature dependent source spectral shift has recently led to a temperature-modulated bioluminescence tomography method which uses a focused ultrasound array [15]. In principle, this should belong to spectrally-resolved a priori information. The a priori permissible source region is defined by the surface photon distribution and improves the reconstruction by constraining the permissible source position [14]. Overall, it is necessary to define additional a priori parameters for BLT reconstruction.

The diffusion approximation is extensively used in BLT reconstruction despite the fact that higher approximations of the radiative transfer equation lead to improved reconstructed results in some situations [16]. The finite element method (FEM), analytical formulations and Born approximation theory have been applied in combination with the diffusion equation [17]. The FEM has become popular due to its ability to process complex heterogeneous geometries. The adaptive strategy has also been developed to further improve the reconstruction based on the FEM [18][19]. In BLT, although nonlinear optimization strategies [20] used in diffuse optical tomography (DOT) and expectation maximization (EM) algorithms [9] similar with that in positron emission tomography (PET) are used, a linear least square (LS) problem is easily obtained because of the linear nature of the BLT problem [14]. Meanwhile, the inverse crime needs to be carefully considered especially when new algorithms are evaluated using synthetic data [21].

BLT reconstruction is an ill-posed problem. Inhibiting noise in measured data and reducing the ill-posedness is necessary to obtain BLT reconstructions. Regularization is a useful method for such problems. Currently, the weighted least square method is used to reduce the measured noise effects [22][14]. The Tikhonov regularization is a popular method and is extensively applied in BLT reconstruction [23][19]. Mathematically, the Tikhonov method is aiming to stabilize the inverse of an ill-conditioned operator by minimizing a trade-off between a loss function and the l ₂-norm of the signals. The advantage of the l ₂ norm is that the associated optimization problem can be efficiently solved using a classic quadratic minimization algorithm. The disadvantage is that the solution obtained is often smoothed everywhere, resulting the loss of high frequency structures of the original signal, especially in the case of noise. Over the past several years, the l ₁ norm regularization has been investigated in the signal and image processing fields, such as wavelet thresholding denoising [24], basis pursuit [25], and total variation for edge preserving reconstruction [26]. Moreover, a new sampling theory related to l ₁ minimization, known as compressed sensing (or compressed sampling) provides a strong theoretical foundation for sparse approximations [27][28]. More accurately, this theory allows an exact reconstruction from a greatly reduced number of measurements through the use of convex programming. In other words, if the real signals or images are sparse on some basis and the measurement operator and sparsity basis satisfy certain coherent conditions, then the original sparse signal can be reconstructed with a greatly reduced sampling rate. In BLT, the unknown sources are contained in the reconstructed domain (such as a mouse). Non-invasive measurements only acquire the surface distribution of photons emitted by bioluminescence sources. When using small elements (such as tetrahedron or hexahedron) to discretize the whole domain, the number of the surface discretized points is significantly fewer than that of the volumetric discretized points. The un-dersampling is inevitable for BLT reconstruction. Compared with single view measurements, multi-view data acquisition improves the BLT reconstruction to a certain degree [29][30], but it limits the high throughput ability of optical imaging. Single view measurements need to be further investigated for improved BLT reconstruction.

Fortunately, when we use optical probes to observe the specified biological process of interest, the domain of the light source is relatively small and sparse compared with the entire reconstruction domain, in this case the mouse body. Here, by a combination of this a priori information and compressed sensing theory, a novel spectrally-resolved bioluminescence reconstruction algorithm is proposed. Specifically, based on the diffusion approximation model, the linear relationship between the spectrally-resolved measured data and the unknown source distribution is established by using the FEM. The l ₁ norm as a regularization term is combined into the BLT least squares problem, realizing the compressed sensing method. In order to reduce memory and time cost, a limited memory variable metric optimization method is used to solve the bound-constrained BLT problem. In numerical verifications, the inverse crime is demonstrated for different synthetic data sets from different finite element meshes and different simulation methods, showing that the Monte Carlo method is necessary for accurate simulation tests. Furthermore, BLT reconstructions with different noise levels and different source depths demonstrate the usefulness of the compressing sensing method-based l ₁ norm regularization, especially for sources located deep within tissues and having high noise. Finally, the proposed algorithm is further tested by experimental reconstruction. In the next section, we present the spectrally-resolved BLT framework based on l ₁ regularization. In the third section, we evaluate the performance of the proposed algorithm with various source settings. In the final section, we discuss the relevant issues and conclusions. To the best of the authors knowledge, this contribution represents the first time that bioluminescence source sparse characteristic is used to improve BLT reconstruction with the compressed sensing method.

2. Model

2.1. Photon diffusion model

The available bioluminescence probes typically emit photons in the range of 400–800nm. The diffusion models as the P ₁ approximation to the radiative transfer equation have been extensively applied in bioluminescence imaging. As the optical properties of biological tissues change depending on the wavelength, the diffusion model is also a function of the wavelength. Assuming the bioluminescence source intensity is stable when photons are collected, the steady-state diffusion equation can be used to depict the photon propagation in tissues:

where Ω and λ is the domain and the wavelength respectively; Φ(r, λ) denotes the photon flux density; S(r,λ) is the source energy density; μ_a(r,λ) is the absorption coefficient; D(r, λ)=1/(3(μ_a(r, λ)+(1 -g)μ_s(r,λ))) is the optical diffusion coefficient, μ_s(r, λ) the scattering coefficient, and g is the anisotropy parameter. On the boundary ∂Ω, the Robin boundary condition is used to depict the refractive index mismatch between the external medium n′ and Ω:

where v is the unit outer normal on ∂Ω.A(r;n,n′) can be approximately represented as:

where n′ is close to 1.0 when the mouse is in air; R(r) can be approximated by R(r) ≈ 1.4399n ^-2 + 0.7099n ^-1 + 0.6681 + 0.0636n [31]. When practical measurements are performed with a set of bandpass filters, the measured quantity is the outgoing flux density Q(r,λ) on the discretized wavelength λ_i, which is:

2.2. Linear Relationship Establishment

Based on the finite element theory [32], the weak solution of the flux density Φ(r, λ_i) ∈H ¹ (Ω) is given considering Eqs. (1) and (2) for a specified wavelength λ_i:

where ∀Ψ(r)∈H ¹(Ω), H ¹(Ω) is the Sobolev space, and Ψ(r) is an arbitrary piece-wise test function. In the numerical finite element computation, the domain Ω needs to be discretized into a group of small elements τ. Correspondingly in three dimensional computations, Ψ(r) is disretized as shape functions regarding the element τ. Tetrahedra and hexahedra are usually used as τ. Regardless of the specified element and shape function, we get the following finite element-based matrix form when the source term is unknown:

where λ_i, K(λ_i), C(λ_i) and B(λ_i) are called the mass, stiff and boundary matrix respectively. These are obtained from the first, second and third term in Eq. 5’s right side [14]. Note that F(λ_i) can be flexibly selected depending on the choice of S(λ_i). Generally, we may select discretized elements or points as the unknown variables. Here, we use the point-based mode. When K(λ_i), C(λ_i) and B(λ_i) are summed as M(λ_i), we have:

Here, we consider the linear relationship between the unknown source variable S(λ_i) and the flux density Φ^b(λ_i) at the measurable boundary discretized points. 𝓕(λ_i) can be adjustably obtained regarding the whole domain, or a priori or a posteriori permissible source region as unknown source region [19]. When the energy percentage of the wavelength λ_i is γ_i, we get:

where

Total energy S at all the wavelengths will be ∑_i=1 ^I γ_i S(λ_i), where I is the total wavelength number.

2.3. Regularization

In optical tomographic imaging, the physical meaning of the various parameters and constraining minimization problem by optimization methods have significant impact on object reconstruction [33]. Therefore, for Eq. 8, we get the following constrained minimization problem for the measured signal Φ^meas which corresponds to Φ^b:

where S^sup is a known upper bound vector and ∥ · ∥ denotes l ₂ norm of a vector. If we ignore the constraints, this also corresponds to solving the normal equation

When the spectrum of the operator 퓐^T퓐 is unbounded or ill-conditioned, the inverse of this equation can cause severe numerical instabilities. A standard procedure is to integrate a priori information in the solution, called regularization. For example, the simplest Tikhonov regular-ization consists of adding a l ₂ norm penalty term to the l ₂ loss functional, i. e.

where δ is a positive number called the regularization parameter which is used to balance the fidelity term and the regularization term. The related gradient of the objective function of Eq. 11 is written as:

This quadratic functional can be efficiently solved by a large range of convex programming technique.

However, this quadratic stabilizer intends to recover a smoothed version of S independent of the data structure. It is often incapable of recovering local singularities or discontinuities presented in the object in the case of noise. We are now considering a non-quadratic l_p norm penalty where 1 ≤ p < 2, since the functional ceases to be convex if p < 1. The main advantage of the non quadratic norm is to promote the sparsity of the solution. Roughly speaking, when p → 1, large components of S are less penalized when compared to the quadratic norm. However, the sum of small components are more penalized, thus leading to a sparse solution. In BLT, an important a priori information is the sparsity of the light source, therefore l ₁ regularization is a more natural choice for this ill posed inverse problem. It follows the model:

where ∥S∥₁ = ∑_i ∣S_i∣ denotes the l ₁ of the vector S. Since this functional is non-differentiable, we can use a differentiable approximation defined as ∥S∥₁ ≈ ∑_i=1 ⁿ F_ε (S(i)) [40] defined as

where ε is a small positive number.

2.4. Algorithm

 

Download Full Size | PDF

Through minimizing the objective function Θ(S), BLT reconstruction can be obtained. Θ(S) is a typical bound-constrained regularization-based least square problem. For the constraint problem, an active-set strategy which includes several types of Hessian matrix based optimization algorithms is adopted to obtain a desirable reconstruction [14][19]. Although this least square problem easily obtains the Hessian matrix, it requires a significant amount of memory during the optimization procedure, especially for large-scale problems. In addition, when computing the search direction, it is necessary to invert the Hessian matrix, which is time-consuming and severely affects the speed of BLT reconstruction. One solution is to use quasi-Newton methods. Generally, these build up an approximate Hessian matrix by using gradients and iteration algorithms. This approximate matrix is obtained in real-time by vector-vector multiplications and is easy to invert, saving memory and time requirements. Here, the limited memory variable metric bound constrained quasi-Newton method (BLMVM) [34] is used for BLT reconstruction. The detailed algorithm is shown in Algorithm 1.

Specifically, an initial guess S ₀ for the source distribution should be given and the initial searching direction d ₀ is also provided. Here, the operator 퓣_퓓 is defined as

where d(j),S(j) denotes the j-th element of d and S respectively. When the step size α_k is determined, the iterative solution at the next step S _k+1 can be calculated through the projection operator 풫 onto the box constraint, defined as

During the minimization procedure, the approximation H _k+1 of the inverse Hessian matrix at the next step is updated when y^T_k s_k > 0

where ρ_k = 1/(y^T_k s_k),V_k = I -ρ_k(y_ks^T_k), I is the identity matrix. Since the inverse Hessian matrix is usually dense, the memory and time requirements for processing H_k is prohibitive especially for large scale problems. In the BLMVM algorithm, a limited memory BFGS matrix is obtained by the vector pairs from the last m iterations. Given an initial inverse Hessian approximation H _k+1 ⁰, the updated matrix H_k is obtained

3. Results

3.1. Simulation verifications

Much attention should always be given to the inverse crime when new algorithms are verified using synthetic data. Monte Carlo methods can simulate the photon propagation better given the ability to incorporate Poisson noise in the simulation. In addition, the same discretized modes used in the forward simulation and inverse reconstruction will significantly affect the evaluation. A cubic domain with a width of 15mm was used to confirm these effects. The synthetic data was obtained through three types of modes, which are hexahedra- and tetrahedra-based FEM, and Monte Carlo method. The diffusion approximation equation was used in the FEM-based simulation. The same discretization with hexahedra-based simulation was used in BLT reconstruction and its element size was 1.0mm in width. The average element diameter in tetrahedra-based discretization was also 1.0mm. Figures 1(a) and 1(b) show the discretized meshes. In reconstruction, three wavelengths (600nm, 650nm and 700nm) were used to obtain spectrally-resolved measurements. We refer to the literature [9] to obtain the corresponding optical properties as listed in Table 1. Photon attenuation is approximately an exponential function of the effective attenuation coefficient ${\mu }_{\mathrm{eff}}\left({\mu }_{\mathrm{eff}}=\sqrt{3{\mu }_a\left({\mu }_S^{\prime }+{\mu }_a\right)}\right).$ In order to preserve the noise effect for all the wavelengths, we sampled 10⁷ photons for 600nm and half the number of these photons were used at two other wavelengths. Monte Carlo simulation is severely time-consuming. To accelerate the simulation, MPI-based parallel code based on the Molecular Optical Simulation Environment (MOSE) [35] was developed in order to perform spectrally-resolved simulations. Because the parallel program only records the information of the photons emitted through the boundary, we can consider that the cubic domain used in the MC method is not discretized, as shown in Fig. 1(c).

Table 1. Optical property at three wavelengths for cubic phantom in simulation verifications

View Table | View all tables in this article

When generating the synthetic data, a cubic source with a width of 1mm was placed at the center of the cubic domain. The source intensity at every wavelength was set to “1.0”. In the reconstruction, we used the synthetic data on the top surface, while the additional noise was not considered. Additionly, regularization methods were not used. Based on three different types of synthetic data, the reconstructed results are shown in Figs. 1(d) to 1(f). When using the hex-ahedral mesh, the same mathematical model and discretized mesh were used in synthetic data generation and reconstruction. Although there are some artifacts in the reconstructed results, the source position is well localized, as seen in Fig. 1(d). However, when the tetrahedral-based synthetic data is used, the reconstructed results become inaccurate. Also, in Fig. 1(e), it is observed that there are some reconstructed results below the actual source position. MC-based reconstructed results in Fig. 1(f) are similar with those based on the tetrahedral mesh. The difference is that the reconstruction around the actual source position becomes increasingly obscure and its position is proximal to the top surface. Actually, the use of different simulation methods results in different levels of noise in the synthetic data, affecting the reconstructed results. Therefore, the synthetic data results are compared in terms of the three different methods and the results are shown in Fig. 2(b). We can find the Hex- and Tet-based synthetic data is almost the same. The maximal relative error (RE) (RE = (Φ_Tet -Φ_Hex)/Φ_Hex) between them is only 10.4%. However, these errors introduce significant effects during reconstruction, showing the ill-posedness of the BLT problem and the necessity of regularization. Furthermore, we can see there are large errors between the Hex- and MC-based data especially when the discretized points are far from the center. It is clear that this produces the poor reconstructed results. From another aspect, it is necessary to use the MC-based synthetic data for testing the proposed algorithm due to its precise simulation and the inverse crime problem. Note that for convenient comparison the MC-based data is normalized using the Hex-based data based on its respective maximal flux density.

 

Fig. 1. Verification to the inverse crime problem. The discretizations of the cubic domain in Figs. (a), (b) and (c) were used to generate the synthetic data using the finite element method (Figs. (a) and (b)) and Monte Carlo method (Fig. (c)). Figures (d), (e) and (f) are the reconstructed results respectively when the real source central position is at (0.0,0.0,0.0). The synthetic data on the top surface is used in reconstruction.

Download Full Size | PDF

 

Fig. 2. Quantitative comparison between HEX-, TET- and MC-based synthetic data at 650nm

Download Full Size | PDF

3.1.1. Single Source Cases with the Homogeneous Media

 

Fig. 3. BLT reconstructions when the real source central position is at (0.0,0.0,0.0) and 10⁶ photons are tracked to generate the synthetic data at 600nm. Figures (a), (b) and (c) are the photons distribution at 600nm, 650nm and 700nm. Figures (d), (e) and (f) are the reconstructed results without regularization method and with l ₂ and l ₁ methods.

Download Full Size | PDF

With the MC-based synthetic data, the proposed algorithm is verified. For each method, we present our results with an optimal parameter δ chosen from a series of values. Generally, based on the reconstruction experience, the ranges of l ₁ and l ₂ parameters are from 10^-7 to 10^-4 and from 10^-6 to 10^-3 respectively. The regularization parameters become larger along with the noise increase. In the single source case, we use the same settings with those used in the inverse crime evaluations and reduce the simulated photon number to 10⁶ at 600nm. Figs. 3(a) to 3(c) show the photon distribution on the top surface of the cubic domain at three wavelengths. It is obvious that they are different because of the effect of the optical properties at different wavelengths. When regularization methods are not used in the reconstruction, we get similar reconstructed results as in Fig. 1(f), and which are shown in Fig. 3(d). We cannot accurately localize the source position. Figure 3(e) shows the reconstructed results when the l ₂ regulariza-tion method is used. The center position of the reconstructed source is at (0.0,0.0,1.5). Due to the effect of the noise on the source depth information contained in the synthetic data, there is an 1.5mm error in depth localization. Similar localization information is obtained when the l ₁ regularization is used. Both of the regularization-based BLT reconstructions show good source localization when compared with and without regularization.

When we reduce the tracking photons number to 10⁴, figures 4(a) to 4(c) show the photon distribution on the top surface. Because few photons are emitted through the boundary, high Poisson noise exists in the synthetic data. It is difficult to distinguish the difference between the three wavelengths. When no regularization method is used, we obtain a degraded reconstruction, which is shown in Fig. 4(d), compared with that obtained by 10⁶ photons. Even if the l ₂ regularization is used, we cannot always accurately localize the source position using the reconstructed results (Fig. 4(e)) no matter how the regularized parameter is adjusted. When the l ₁ regularization is used, a similar reconstruction, as shown in Fig. 4(f) as that with 10⁶ photons is obtained. Due to the higher noise level, the center of the reconstructed source is at (0.0,0.0,2.0) and the localization errors further are increased further. However, compared with the l ₂ regu-larization, the l ₁ method shows improvement especially when the source is at a deeper position and high noise exists in the measured data.

 

Fig. 4. BLT reconstructions when the real source central position is at (0.0,0.0,0.0) and 10⁴ photons are tracked to generate the synthetic data at 600nm. Figures (a), (b) and (c) are the photons distribution at 600nm, 650nm and 700nm. Figures (d), (e) and (f) are the reconstructed results without regularization method and with l ₂ and l ₁ methods.

Download Full Size | PDF

3.1.2. Dual source cases with the homogeneous media

In this simulated reconstruction, dual source settings are considered in order to evaluate the l ₁ regularization method. Both sources have the same settings as those used in the single source cases and are placed at (-3.0,0.0,3.0) and (3.0,0.0,3.0) respectively. First, 10⁶ photons are tracked at 600nm for each source. Since the sources are close to the top surface and many photons can be captured, we obtain good reconstructed results without regularization, which is shown in Fig. 5(a). The center positions of the reconstructed sources are (-2.5,0.0,3.5) and (2.5,0.0,3.5). Note that even if we set the regularization parameters with the l ₂ and l ₁ methods, similar results (Figs. 5(b) and 5(c)) are obtained with those without regularization, illustrating the robust nature of the regularization methods. When the photon number reduces to 10⁴, the reconstructed results are shown in Figs. 5(d) to 5(f). Obviously, without regularization, we cannot obtain accurate source localization due to the high noise in the synthetic data. The l ₂- and l ₁-based reconstructions show similar source localization with those using 10⁶ photons. Note that these reconstructions are similar to those in the single source case with 10⁶ photons. In other words, more photons are required for BLT reconstruction with multiple sources compared with single sources. This is true even if these sources are closer to the measured surface than the latter case.

 

Fig. 5. Dual source BLT reconstructions when the real source central positions is at (-3.0,0.0,3.0) and (3.0,0.0,3.0). Figures (a), (b) and (c) are the corresponding reconstructed results without regularization and with l ₂ and l ₁ methods when 10⁶ photons are tracked at 600nm. Figures (d), (e) and (f) are the counterparts corresponding to (a), (b) and (c) when 10⁴ photons are tracked at 600nm.

Download Full Size | PDF

When the two sources are moved to (-3.0,0.0,0.0) and (3.0,0.0,0.0), Figures 6(a) to 3(f) show the reconstructed results when 10⁶ photons are tracked. We cannot distinguish the source position accurately without regularization and with the l ₂ method, although the lower values in the latter BLT reconstruction show that there are two sources. In contrast, two sources can be distinguished from the reconstruction with the l ₁ method despite the fact that the localization errors (the reconstructed center positions are (-2.0,0.0,2.5) and (2.0,0.0,2.5)), are bigger than those in the single source case (Fig. 3(f)). Based on the synthetic data with 10⁴ photons, the reconstruction results are shown in Figs. 6(d) to 4(f). Due to the higher noise, we can’t distinguish two sources (Fig. 4(e)) in the l ₂-based reconstruction, or without regularization (Fig. 4(d)). Surprisingly, two sources can be distinguished in the l ₁-based reconstruction. However, the localization errors become bigger than those in the 10⁶ photon case.

 

Fig. 6. Dual source BLT reconstructions when the real source central positions is at (-3.0,0.0,0.0) and (3.0,0.0,0.0). Figures (a), (b) and (c) are the corresponding reconstructed results without regularization and with l ₂ and l ₁ methods when 10⁶ photons are tracked at 600nm. Figures (d), (e) and (f) are the counterparts corresponding to (a), (b) and (c) when 10⁴ photons are tracked at 600nm.

Download Full Size | PDF

3.1.3. Dual source cases with the heterogeneous media

Furthermore, heterogeneous BLT reconstructions with dual source settings are performed. Two sources are placed at (-3.0,0.0,0.0) and (3.0,0.0,0.0) respectively. The heterogeneous characteristics of the domain are realized by placing a 5×5×5 cube within the homogeneous domain. The center of this cube is the same as that of the source at (-3.0,0.0,0.0). The absorption and reduced scattering coefficients at three wavelengths are set to 0.038 and 1.82, 0.015 and 1.73, and 0.004 and 1.57 respectively. When 10⁶ photons are used to generate the measured data at 600nm, the reconstructed results are shown in Figs. 7(a) to 7(c). They are similar with the BLT reconstructions in the homogeneous domain and two sources can not be distinguished without regularization and with thel ₂ method. When the l ₁ regularization is used, the reconstructed results have little difference between using heterogeneous and homogeneous domains, as shown in Figs. 6(c) and 7(c). Due to the heterogeneous media characteristics, the reconstructed positions become (-1.5,0.45,2.55) and (1.5,-0.45,2.55) respectively. Although the depth localization is similar in heterogeneous and homogeneous domains, the reconstructed precisions at X and Y directions become worse. When the photon number is reduced to 10⁴, the reconstructed results are shown in Figs. 7(d) to 7(f). Note that the depth localization is improved with the l ₁ regularization due to the heterogeneous media characteristics, and the reconstructed positions are (-1.35,-0.6,2.7) and (1.8,-0.6,3.45). With the reduced number of photons, the heterogeneous characteristics improve the reconstruction precision.

 

Fig. 7. Dual source BLT reconstructions with the heterogeneous media when the real source central positions is at (-3.0,0.0,0.0) and (3.0,0.0,0.0). Figures (a), (b) and (c) are the corresponding reconstructed results without regularization and with l ₂ and l ₁ methods when 10⁶ photons are tracked at 600nm. Figures (d), (e) and (f) are the counterparts corresponding to (a), (b) and (c) when 10⁴ photons are tracked at 600nm.

Download Full Size | PDF

3.1.4. Multiple source cases

With respect to different source intensities, three sources with different depths are set in the homogeneous domain to test the l ₁-based reconstruction method. Their positions and intensities are (-2.0,2.0,4.0) and 1.0, (0.0,0.0,0.0) and 5.0, and (3.0,-3.0,2.0) and 3.0 respectively. In this simulation, 10⁴ photons are tracked at 600nm. The reconstructed results are shown in Figs. 8(a) to 8(c). Without regularization methods, the three sources cannot be distinguished (Fig. 8(a)). When the l ₂ and l ₁ methods are used, the three sources can be distinctly distinguished. However, overall there is a coupling between the source depth and intensity, and the source depth localizations become worse, while the source intensities cannot be reconstructed accurately. When comparing the reconstructed results based on the l ₂ and l ₁ methods, there is an artifact in the l ₂-based reconstruction, which is shown in Fig. 8(b). To obtain good source depth and intensity reconstruction, more sophisticated l ₁-based reconstruction algorithms should be developed.

3.2. Experimental data reconstructions

In order to test the spectrally-resolved l ₁-based BLT reconstruction algorithm, a commercially available solid mouse shaped homogeneous phantom was used. This phantom was fabricated with a polyester resin, TiO₂ and Disperse Red. Table 2 shows the optical properties (μ_a and (μ′_s) at six wavelengths measured with the inverse adding doubling method [36]. To imitate the bioluminescence source, an optical fiber coupled to a green LED was embedded within the phantom. The emission spectrum of the LED ( wavelength range was from 500nm to 700nm and the peak was at about 567nm) was similar to that of a bioluminescence source. More detailed information about this phantom can be obtained elsewhere in [22]. To acquire the shape of this phantom and localize the source position, an Imtek microCAT system (Siemens Preclinical Solutions, Knoxville, TN) was used. Since the diameter of the optical fiber was only 200μm, the source could be considered as a point source. GFP (515–575nm) and DsRed (575–650nm) emission filters were used to acquire the spectrally-solved measured data. In the reconstruction, the optical properties at 560nm were used for the GFP filter-based data. The averaged μ_a and μ′_s from 580nm to 660nm (0.013mm ^-1 and 1.68mm ^-1) was considered for the DrRed filter-based measurement. To avoid the effects of the curved surface in the measured data, the photon distribution was obtained from the bottom surface of the phantom in a Caliper IVIS 100 imaging system (Caliper Life Sciences Alameda, CA). Using the commercial software Amira 3.0 (Mercury Computer Systems, Inc. Chelmsford, MA), the tetrahedral-based finite element volumetric mesh shown in Fig. 10(a) for reconstruction was generated based on the CT images. With respect to the photons distribution, about 2/3 of the whole phantom was selected for mesh generation. The mesh had the average element diameter of 3.0mm and included 1929 nodes and 7766 elements. The measured data was manually registered using the simultaneously obtained photograph in Amira.

 

Fig. 8. Triple source BLT reconstructions with the homogeneous media when the real source central positions is at (-2.0,2.0,4.0) , (0.0,0.0,0.0), and (3.0,-3.0,2.0). Figures (a), (b) and (c) are the corresponding reconstructed results without regularization and with l ₂ and l ₁ methods when 10⁴ photons are tracked at 600nm.

Download Full Size | PDF

Table 2. Optical properties of Caliper mouse phantom at six wavelengths

View Table | View all tables in this article

The photon distributon for 2 minute acquisitions using two types of filters is shown in Fig. 9. Since the optical properties at two wavelength ranges are different, we can clearly observe the difference in the photon distribution. Using the CT images, the actual source position was localized at (114.5,131.0,3.0). When the BLT reconstruction without regularization was performed, the reconstructed results are shown in Fig. 10(b), indicating a distributed source. The difference between the corresponding maximal values in several distributed regions is small. When the maximal reconstructed values are used to decide the reconstructed position, its center is localized at (111.7,132.6,2.7). Overall, there are large reconstruction errors especially along X-axis direction. When the l ₂ regularization-based reconstruction was performed, the reconstructed results are shown in Fig. 10(c). Due to the smoothness function of the l ₂ regularization, the whole reconstructed region is almost filled by the values close to the maximal. However, the center position of the reconstructed results is easily localized at (115.1,131.7,2.4). The errors at three axes are 0.6, 0.7, 0.6 compared with the actual position. Compared with the reconstruction without regularization, this localization is better. Figure 10(d) shows the reconstructed results with l ₁ regularization, and the reconstructed position (114.7,131.7,2.9) is similar with that based on l ₂ method. However, the reconstructed results are compact, which shows that the BLT reconstruction is significantly improved when sparse a priori information is used.

 

Fig. 9. Surface radiance images of the mouse-shaped phantom with embedded fiber optic source using GFP and DsRed emission filters.

Download Full Size | PDF

 

Fig. 10. Experimental BLT reconstructions with mouse-shaped phantom. Figure (a) are the volumetric mesh used in reconstruction. Figures (b), (d) and (d) show the reconstructed results without regularization and with l ₂ and l ₁ methods.

Download Full Size | PDF

4. Conclusion

In this paper, a spectrally-resolved l ₁ regularization based reconstruction algorithm is proposed. Based on the linear relationship between the unknown source variable and the boundary measurements, the spectrally-resolved information is used to improve BLT reconstruction. Sparse source characteristics are considered in a graceful way along with the l ₁ regularization method. The use of quasi-Newton optimization methods accelerates the BLT reconstruction. Simulation verifications with MC-based synthetic data show that the use of the sparse a priori information significantly improves the BLT reconstruction compared with the popular l ₂ reg-ularization. Particularly the case when the sources exist at deeper positions and the measured data contains high noise, the l ₁-based methods are necessary to obtain improved location reconstruction. Reconstruction of experimental data further shows the effectiveness of the proposed algorithm.

In BLT reconstruction, it is vital to solve the ill-posed and unique problems. Regularization and a priori information significantly improve the reconstruction. From the results obtained here, reconstruction without regularization is not stable. Although l ₂ regularization improves performance, its smooth characteristic is not apropriate for the BLT problem. Since the l ₁ reg-ularization considers sparse information, it is more suitable for BLT reconstruction especially when multiple sources exist. Note that to strictly meet the condition of compressed sensing theory, signal incoherence should be considered further. In practice, for a more general inverse problem, it is not simple to check the incoherence condition. Perfect matching between theory and practice is limited in some simple cases [37]. It is still standard to use a nonquadratic norm to promote sparsity, roughly speaking, for most large under-determined systems of linear equations, the minimal of l ₁ norm solution is also the sparsest solution [38]. Fortunately, the improved reconstructed results show the inherent characteristics of BLT problem can meet the compressed sensing theory to a certain extent.

Although in vivo mouse BLT reconstructions with the diffusion approximation theory obtain good results, the nature of photon propagation in biological tissues demonstrates that more precise mathematical models should be considered for the BLT problem when the biolumines-cence source is in or close to specific tissues. Several improved models have been proposed to improve the reconstruction quality [16][39]. Since BLT reconstruction is a linear inverse source problem in nature, as a regularization method, the l ₁ method can be used in improved precise model based reconstructions.

In conclusion, we have developed a spectrally-resolved compressed sensing based reconstruction method for BLT, obtained encouraging preliminary results in both numerical simulations and physical phantom experiments, and established that our proposed method is effective for BLT. In vivo mouse studies using the proposed method will be reported in the future.