## Abstract

Through restoration of the light source information in small animals *in vivo*, optical molecular imaging, such as fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT), can depict biological and physiological changes observed using molecular probes. *A*
*priori* information plays an indispensable role in tomographic reconstruction. As a type of *a priori* information, the sparsity characteristic of the light source has not been sufficiently considered to date. In this paper, we introduce a compressed sensing method to develop a new tomographic algorithm for spectrally-resolved bioluminescence tomography. This method uses the nature of the source sparsity to improve the reconstruction quality with a regularization implementation. Based on verification of the *inverse crime*, the proposed algorithm is validated with Monte Carlo-based synthetic data and the popular Tikhonov regularization method. Testing with different noise levels and single/multiple source settings at different depths demonstrates the improved performance of this algorithm. Experimental reconstruction with a mouse-shaped phantom further shows the potential of the proposed algorithm.

©2009 Optical Society of America

## 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*
_{2}-norm of the signals. The advantage of the *l*
_{2} 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*
_{1} 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*
_{1} 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*
_{1} 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*
_{1} 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*
_{1} 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*
_{1} 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.4399*n*
^{-2} + 0.7099*n*
^{-1} + 0.6681 + 0.0636*n* [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*

^{1}(Ω) is given considering Eqs. (1) and (2) for a specified wavelength

*λ*:

_{i}$$+{\int}_{\partial \Omega}\genfrac{}{}{0.1ex}{}{1}{2{A}_{n}\left(\mathbf{r}\right)}\Phi \left(\mathbf{r},{\lambda}_{i}\right)\psi \left(\mathbf{r}\right)d\mathbf{r}={\int}_{\Omega}S\left(\mathbf{r},{\lambda}_{i}\right)\psi \left(\mathbf{r}\right)d\mathbf{r}$$

where ∀Ψ(**r**)∈*H*
^{1}(Ω), *H*
^{1}(Ω) 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*(

*λ*) and

_{i}*B*(

*λ*) 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

_{i}*F*(

*λ*) can be flexibly selected depending on the choice of

_{i}*S*(

*λ*). Generally, we may select discretized elements or points as the unknown variables. Here, we use the point-based mode. When

_{i}*K*(

*λ*),

_{i}*C*(

*λ*) and

_{i}*B*(

*λ*) are summed as

_{i}*M*(

*λ*), we have:

_{i}Here, we consider the linear relationship between the unknown source variable *S*(*λ _{i}*) and the flux density Φ

*(*

^{b}*λ*) at the measurable boundary discretized points. 𝓕(

_{i}*λ*) can be adjustably obtained regarding the whole domain, or

_{i}*a priori*or a

*posteriori*permissible source region as unknown source region [19]. When the energy percentage of the wavelength

*λ*is

_{i}*γ*, we get:

_{i}where

Total energy *S* at all the wavelengths will be ∑_{i=1}
^{I}*γ _{i}*

*S*(

*λ*), where

_{i}*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*

_{2}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*

_{2}norm penalty term to the

*l*

_{2}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*

_{1}regularization is a more natural choice for this ill posed inverse problem. It follows the model:

where ∥*S*∥_{1} = ∑* _{i}* ∣

*S*∣ denotes the

_{i}*l*

_{1}of the vector

*S*. Since this functional is non-differentiable, we can use a differentiable approximation defined as ∥

*S*∥

_{1}≈ ∑

_{i=1}

^{n}*F*(

_{ε}*S*(

*i*)) [40] defined as

where *ε* is a small positive number.

#### 2.4. Algorithm

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*
_{0} for the source distribution should be given and the initial searching direction *d*
_{0} 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*> 0

_{k}where *ρ _{k}* = 1/(

*y*

^{T}_{k}*s*),

_{k}*V*=

_{k}*I*-

*ρ*(

_{k}*y*),

_{k}s^{T}_{k}*I*is the identity matrix. Since the inverse Hessian matrix is usually dense, the memory and time requirements for processing

*H*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

_{k}*m*iterations. Given an initial inverse Hessian approximation

*H*

_{k+1}

^{0}, the updated matrix

*H*is obtained

_{k}$$+{\rho}_{k-m}\left({V}_{k-1}^{T}\dots {V}_{k-m+1}^{T}\right){S}_{k-m}{S}_{k-m}^{T}\left({V}_{k-m+1}\dots {V}_{k-1}\right)$$

$$\text{}\vdots $$

$$+{\rho}_{k-2}{\left({V}_{k-1}\right)}^{T}{S}_{k-2}{S}_{k-1}^{T}\left({V}_{k-1}\right)$$

$$+{\rho}_{k-1}{S}_{k-1}{S}_{k-1}^{T}$$

## 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 15*mm* 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.0*mm* in width. The average element diameter in tetrahedra-based discretization was also 1.0*mm*. Figures 1(a) and 1(b) show the discretized meshes. In reconstruction, three wavelengths (600*nm*, 650*nm* and 700*nm*) 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^{7} 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).

When generating the synthetic data, a cubic source with a width of 1*mm* 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}*) 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*

_{Hex}*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.

### 3.1.1. Single Source Cases with the Homogeneous Media

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*
_{1} and *l*
_{2} 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^{6} at 600*nm*. 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*
_{2} 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.5*mm* error in depth localization. Similar localization information is obtained when the *l*
_{1} 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^{4}, 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^{6} photons. Even if the *l*
_{2} 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*
_{1} regularization is used, a similar reconstruction, as shown in Fig. 4(f) as that with 10^{6} 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*
_{2} regu-larization, the *l*
_{1} method shows improvement especially when the source is at a deeper position and high noise exists in the measured data.

### 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*
_{1} 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^{6} photons are tracked at 600*nm* 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*
_{2} and *l*
_{1} 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^{4}, 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*
_{2}- and *l*
_{1}-based reconstructions show similar source localization with those using 10^{6} photons. Note that these reconstructions are similar to those in the single source case with 10^{6} 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.

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^{6} photons are tracked. We cannot distinguish the source position accurately without regularization and with the *l*
_{2} 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*
_{1} 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^{4} 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*
_{2}-based reconstruction, or without regularization (Fig. 4(d)). Surprisingly, two sources can be distinguished in the *l*
_{1}-based reconstruction. However, the localization errors become bigger than those in the 10^{6} photon case.

### 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^{6} photons are used to generate the measured data at 600*nm*, 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 the*l*
_{2} method. When the *l*
_{1} 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^{4}, the reconstructed results are shown in Figs. 7(d) to 7(f). Note that the depth localization is improved with the *l*
_{1} 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.

### 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*
_{1}-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^{4} photons are tracked at 600*nm*. 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*
_{2} and *l*
_{1} 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*
_{2} and *l*
_{1} methods, there is an artifact in the *l*
_{2}-based reconstruction, which is shown in Fig. 8(b). To obtain good source depth and intensity reconstruction, more sophisticated *l*
_{1}-based reconstruction algorithms should be developed.

#### 3.2. Experimental data reconstructions

In order to test the spectrally-resolved *l*
_{1}-based BLT reconstruction algorithm, a commercially available solid mouse shaped homogeneous phantom was used. This phantom was fabricated with a polyester resin, TiO_{2} and Disperse Red. Table 2 shows the optical properties (*μ _{a}* and (

*μ′*) 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 500

_{s}*nm*to 700

*nm*and the peak was at about 567

*nm*) 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–575

*nm*) and DsRed (575–650

*nm*) emission filters were used to acquire the spectrally-solved measured data. In the reconstruction, the optical properties at 560

*nm*were used for the GFP filter-based data. The averaged

*μ*and

_{a}*μ′*from 580

_{s}*nm*to 660

*nm*(0.013

*mm*

^{-1}and 1.68

*mm*

^{-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.0

*mm*and included 1929 nodes and 7766 elements. The measured data was manually registered using the simultaneously obtained photograph in Amira.

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*
_{2} regularization-based reconstruction was performed, the reconstructed results are shown in Fig. 10(c). Due to the smoothness function of the *l*
_{2} 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*
_{1} regularization, and the reconstructed position (114.7,131.7,2.9) is similar with that based on *l*
_{2} method. However, the reconstructed results are compact, which shows that the BLT reconstruction is significantly improved when sparse *a priori* information is used.

## 4. Conclusion

In this paper, a spectrally-resolved *l*
_{1} 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*
_{1} 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*
_{2} reg-ularization. Particularly the case when the sources exist at deeper positions and the measured data contains high noise, the *l*
_{1}-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*
_{2} regularization improves performance, its smooth characteristic is not apropriate for the BLT problem. Since the *l*
_{1} 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*
_{1} 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*
_{1} 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.

## Acknowledgments

We gratefully thank Dr. Richard Taschereau for his useful discussion about parallel Monte Carlo method. This work is supported by the NIBIB R01-EB001458, NIH/NCI 2U24 CA092865 cooperative agreement, DE-FC02-02ER63520 grants, NSF DMS 0312222, ONR N00014-06-1-0345, NSF CCF-0528583 and the Project for the National Basic Research Program of China (973) under Grant No. 2006CB705700.

## References and links

**1. **R. Weissleder and U. Mahmood, “Molecular imaing,” Radiology **219**, 316–333 (2001). [PubMed]

**2. **V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weisslder, “Looking and listening to light: the evolution of whole body photonic imaging,” Nat. Biotechnol. **23**, 313–320 (2005). [CrossRef] [PubMed]

**3. **J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery **7**, 591–607 (2008). [CrossRef]

**4. **C. H. Contag and M. H. Bachmann, “Advances in bioluminescence imaging of gene expression,” Annu. Rev. Biomed. Eng. **4**, 235–260 (2002). [CrossRef] [PubMed]

**5. **S. Bhaumik and S. S. Gambhir, “Optical imaging of renilla luciferase reporter gene expression in living mice,” Proc. Natl. Acad. Sci. USA **99**, 377–382 (2002). [CrossRef]

**6. **T. F. Massoud and S. S. Gambhir, “Molecular imaging in living subjects: seeing fundamental biological processes in a new light,” Genes Dev. **17**, 545–580 (2003). [CrossRef] [PubMed]

**7. **G. Wang, E. A. Hoffman, G. McLennan, L. V. Wang, M. Suter, and J. F. Meinel, “Development of the first bioluminescence CT scanner,” Radiology **566**, 229 (2003).

**8. **G. Wang, Y. Li, and M. Jiang, “Uniqueness theorems in bioluminescence tomography,” Med. Phys. **31**, 2289–2299 (2004). [CrossRef] [PubMed]

**9. **G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. **50**, 4225–4241 (2005). [CrossRef] [PubMed]

**10. **Y. Lv, J. Tian, W. Cong, and G. Wang, “Experimental study on bioluminescence tomography with multimodality fusion,” Int. J. Biomed. Img. **1**, 86741 (2007).

**11. **C. Kuo, O. Coquoz, D. G. Stearns, and B. W. Rice, “Diffuse luminescence imaging tomography of in vivo bioluminescent markers using multi-spetral data,” in *Proceedings of the 3rd International Meeting of the Society* (MIT Press,2004), p. 227.

**12. **A. J. Chaudhari, F. Darvas, J. R. Bading, R. A. Moats, P. S. Conti, D. J. Smith, S. R. Cherry, and R. M. Leahy, “Hyperspectral and multispectral bioluminescence optical tomography for small animal imaging,” Phys. Med. Biol. **50**, 5421–5441 (2005). [CrossRef] [PubMed]

**13. **H. Dehghani, S. C. Davis, S. Jiang, B. W. Pogue, K. D. Paulsen, and M. S. Patterson, “Spectrally resolved bioluminescence optical tomography,” Opt. Lett. **31**, 365–367 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=ol-31-3-365. [CrossRef] [PubMed]

**14. **W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. V. Wang, E. A. Hoffman, G. McLennan, P. B. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express **13**, 6756–6771 (2005), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-13-18-6756. [CrossRef] [PubMed]

**15. **G. Wang, H. Shen, W. Cong, S. Zhao, and G. Wei, “Temperature-modulated bioluminescence tomography,” Opt. Express **14**, 7852–7871 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-17-7852. [CrossRef] [PubMed]

**16. **A. D. Klose, V. Ntziachristos, and A. H. Hielscher, “The inverse source problem based on the radiative transfer equation in optical molecular imaging,” J. Comput. Phys. **202**, 323–345 (2005). [CrossRef]

**17. **A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. **50**, R1–R43 (2005). [CrossRef]

**18. **Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express **14**, 8211–8223 (2006), http://www.opticsinfobase.org/abstract.cfm?URI=oe-14-18-8211. [CrossRef] [PubMed]

**19. **Y. Lv, J. Tian, W. Cong, G. Wang, W. Yang, C. Qin, and M. Xu, “Spectrally resolved bioluminescence tomography with adaptive finite element analysis: methodology and simulation,” Phys. Med. Biol. **52**, 4497–4512 (2007). [CrossRef] [PubMed]

**20. **X. Gu, Q. Zhang, L. Larcom, and H. Jiang, “Three-dimensional bioluminescence tomography with model-based reconstruction,” Opt. Express **12**, 3996–4000 (2004), http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-12-17-3996. [CrossRef] [PubMed]

**21. **S. Holder, *Electrical Impedance Tomography*, (Institute of Physics Publishing, Bristol and Philadelphia, 2005).

**22. **C. Kuo, O. Coquoz, T. L. Troy, H. Xu, and B. W. Rice, “Three-dimensional reconstruction of *in vivo* biolumines-cent sources based on multispectral imaging,” J. Biomed. Opt. **12**, 024007 (2007). [CrossRef] [PubMed]

**23. **W. Cong, K. Durairaj, L. V. Wang, and G. Wang, “A Born-type approximation method for bioluminescence tomography,” Med. Phys. **33**, 679–686 (2006). [CrossRef] [PubMed]

**24. **D. Donoho and I. Johnstone, “Ideal spatial adaption via wavelet shrinkage,” Biometrika **81**, 425–455 (1994). [CrossRef]

**25. **S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process. **41**, 3397–3415 (2002). [CrossRef]

**26. **L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” J. Phys. D **60**, 259–268 (1992). [CrossRef]

**27. **E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from incomplete and inaccurate measurements,” Commun. Pur. Appl. Math. **59**, 1207–1223 (2006). [CrossRef]

**28. **E. J. Candès, “Compressive sampling,” in *Proceedings of the International Congress of Mathematicians, Madrid, Spain***3**, 1433–1452 (2006).

**29. **C. Kuo, O. Coquoz, T. Troy, D. Zwarg, and B. Rice, “Bioluminescent tomography for *in vivo* localization and quantification of luminescent sources from a multiple-view imaging system,” Mol. Img. **4**, 370 (2005).

**30. **G. Wang, H. Shen, K. Durairaj, X. Qian, and W. Cong, “The first bioluminescence tomography system for simultaneous acquisition of multiview and multispectral data,” Int. J. Biomed. Img. **2006:Article ID 58601**, (2006).

**31. **M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: Boundary and source conditions,” Med. Phys. **22**, 1779–1792 (1995). [CrossRef] [PubMed]

**32. **S. S. Rao, *The finite element method in engineering*, (Butterworth-Heinemann, Boston, 1999).

**33. **R. Roy and E. M. Sevick-Muraca, “Active constrained truncated newton method for simple-bound optical tomography,” J. Opt. Soc. Am. A **17**, 1627–1641 (2000), http://www.opticsinfobase.org/abstract.cfm?URI=josaa-17-9-1627. [CrossRef]

**34. **S. J. Benson and J. Morè, “A limited-memory variable-metric algorithm for bound-constrained minimization,” Technical Report ANL/MCS-P909-0901, Mathematics and Computer Science Division, Argonne National Lab-
oratory (2001).

**35. **H. Li, J. Tian, F. Zhu, W. Cong, L. V. Wang, E. A. Hoffman, and G. Wang, “A mouse optical simulation envi-roment (MOSE) to investigate bioluminescent phenomena in the living mouse with the Monte Carlo Method,” Acad. Radiol. **11**, 1029–1038 (2004). [CrossRef] [PubMed]

**36. **S. A. Prahl, M. J. C. van Gemert, and A. J. Welch, “Determining the optical properties of turbid mediaby using the addingdoubling method,” Appl. Opt. **32**, 559–568 (1993), http://www.opticsinfobase.org/abstract.cfm?URI=ao-32-4-559. [CrossRef] [PubMed]

**37. **
Michael Lustig, *Sparse MRI*, PhD thesis, Stanford University, 2008.

**38. **D. Donoho, “For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution,” Commun. Pur. Appl. Math. **59**, 797–829 (2006). [CrossRef]

**39. **A. D. Klose and B. Beattie, “Bioluminescence tomography with SP_{3} equations,” in *Biomedical Optics Topical Meeting*, 2008, http://www.opticsinfobase.org/abstract.cfm?URI=BIOMED-2008-BMC8.

**40. **J. Cai, S. Osher, and Z. Shen, “Convergence of the Linearized Bregman Iteration for *l*_{1}-Norm Minimization,” CAM report 08–52, 2008.