Novel l<sub>2,1</sub>-norm optimization method for fluorescence molecular tomography reconstruction

Shixin Jiang; Jie Liu; Yu An; Guanglei Zhang; Jinzuo Ye; Yamin Mao; Kunshan He; Chongwei Chi; Jie Tian

doi:10.1364/BOE.7.002342

1. Introduction

In recent years, fluorescence molecular imaging (FMI) has attracted great attention due to the increasing availability of fluorescent proteins, dyes, and probes, which enable the noninvasive study of gene expression, protein function, interactions, and a large number of cellular processes [1, 2]. Based on FMI, fluorescence molecular tomography (FMT) can three-dimensionally (3-D) visualize in vivo fluorescence bio-distribution in small animal studies, which has been widely used for tumor diagnosis and drug discovery [3–6]. FMT reconstruction is an inverse problem, where the fluorescent source is obtained by the system matrix and measurement data sets. Therefore, how to precisely and quickly solve the FMT problem is one of the most challenging problems in FMT studies [7].

FMT reconstruction is an ill-posed problem since only the photon distribution on the surface is measurable. Although the problem can be mitigated by increasing the measurement data sets, such as increasing the number of excitation angles, it is still difficult to obtain satisfactory results because of the data interpolation errors and charge-coupled device (CCD) measurement errors caused by the shot noise of the CCD camera [8–10]. In consequence, FMT reconstruction results are unstable and sensitive to noise. Different methods have been designed to achieve a meaningful approximate solution. Regularization is typically used to tackle the inverse problem. Amongst the traditional regularization methods, l₂-norm regularization (also known as Tikhonov regularization) is the most commonly used method due to its convenient implementation [11, 12]. The primary benefit of using Tikhonov regularization is the simplicity of the optimization problem involved, which can be efficiently solved by standard minimization tools, such as Newton’s method and conjugate gradient method [13–15]. However, the image obtained by Tikhonov regularization turns out to be over-smoothed when we penalize large values in this method. As a result, sharp boundaries of the reconstructed images in general are difficult to obtain via Tikhonov regularization [16].

Since fluorescent sources are usually sparsely distributed and small in size compared with the entire reconstruction domain [17], sparsity regularization methods for FMT reconstruction have been studied extensively [18, 19]. To overcome the over-smooth limitation of the l₂-norm regularization, many recent methods have been proposed to exploit the sparsity of the reconstructed image [20–22], such as l₁ and smooth-l₀ norms. While these methods have been proved superior in overcoming the over-smoothing limitation of the traditional l₂-norm regularization, the inverse problem has been made more difficult to solve due to the non-convexity and non-smoothness of such norms [23].

While continuing efforts are made to develop efficient reconstruction methods for FMT, many challenging problems remain to be solved. Measurement noises are unavoidable in experiments. The standard sparsity based methods are not robust enough in presence of measurement errors. As a result, it is often difficult to distinguish true signals from the relatively high intensity noises. However, such limitations existing in standard sparsity based methods can be overcome by an advanced sparsity technique called structured sparsity [24]. The advantage of structured sparsity based methods is that it can still perform well when the measurement data sets are very limited. Moreover, structured sparsity based methods are also more robust to noise compared to the standard sparsity methods, since the standard sparsity based methods only exploit the sparseness of the reconstructed image, whereas the structured sparsity utilizes correlations or structures of the reconstructed image [25, 26].

As mentioned before, fluorescent signals usually have a similar property of sparse signals based on the fact that the domains of the fluorescent sources are small and sparse compared with the entire reconstruction domain in FMT, because the size of early-stage tumors tagged with the fluorescent probes are small. Moreover, fluorescent sources possess some special features in their structure. The domains of the fluorescent sources are clustered together. The region where the fluorescent sources are clustered has nonzero values while the rest of the regions where other tissues are located all have zero value. Therefore, it is reasonable to assume that the fluorescent sources have a group structure [27, 28]. In this paper, we demonstrate a fast and effective reconstruction method based on group structured sparsity. Considering the group structure of the fluorescent region, the mixed l_2,1-norm was used [29]. The l_2,1-norm is the sum of the l₂-norms of the blocks of coefficients associated with each component. However, due to non-differentiability of the mixed l_2,1-norm, to accurately and effectively solve the l_2,1-norm model, Nesterov’s method [30, 31] was used. Because the l_2,1-norm combines the advantages of L1-regularization and L2-regularization, the reconstruction of the l_2,1-norm is neither over-smoothed nor over-shrunk. Thus, l_2,1-norm is an effective reconstruction strategy for FMT. Comprehensive numerical simulations and in vivo experiments show that our proposed method is more accurate, efficient, and robust for fluorescence reconstruction compared to the Tikhonov-L2 method and IS_L1 method.

This paper is organized as follows. The forward diffusion approximation model and the l_2,1-norm optimization model of FMT are presented in section 2. In section 3, numerical heterogeneous phantom experiments were conducted to evaluate the performance of our proposed method. Then, an in vivo bead-implanted mouse experiment was conducted to demonstrate the feasibility of l_2,1-norm in an in vivo application. Finally, we discuss the results and conclude the results in Section 4.

2. Method

2.1 Photon propagation model

For steady-state FMT with point excitation sources, the photon propagation model in highly scattering media can be described by the following coupled diffusion equation [32, 33]:

{\begin{array}{l} - \nabla [D_{x} (r) \nabla Φ_{x} (r)] + μ_{a x} (r) Φ_{x} (r) = Θ δ (r - r_{l}) & (r \in Ω) \\ - \nabla [D_{m} (r) \nabla Φ_{m} (r)] + μ_{a m} (r) Φ_{m} (r) = Φ_{x} (r) η μ_{a f} (r) & (r \in Ω) \end{array}

where the subscripts

x

and

m

denote the excitation and emission wavelength, respectively;

Ω

is the imaging domain of the problem;

Φ_{x, m}

is the photon flux density;

μ_{a x, a m}

is the absorption coefficient;

D_{x, m} = {[3 μ_{a x, a m} + 3 μ_{s x, s m} (1 - g)]}^{- 1}

is the diffusion coefficient;

μ_{s x, s m}

is the scattering coefficient,

g

is the anisotropy parameter;

η μ_{a f}

denotes the fluorescent field which is to be reconstructed and

Θ

denotes the amplitude of the point source. To solve the equations, the Robin-type boundary conditions are imposed on boundary

\partial Ω

of the imaging domain

Ω

:

\begin{matrix} Φ_{x, m} (r) + 2 q D_{x, m} (r) [\bar{v} (r) \cdot Φ_{x, m} (r)] = 0 & (r \in \partial Ω) \end{matrix}

where

\bar{v}

is the unit outward norm vector on the boundary.

q

is a constant which depends on the optical reflective index mismatch on the boundary [34].

Based on the finite element theory, Eqs. (1) and (2) can be transformed into the following matrix-form equations:

K_{x} Φ_{x} = S_{x}

K_{m} Φ_{m} = F X

where

X

denotes the unknown fluorescent source vector to be reconstructed.

K_{x}

denotes the system matrix during excitation and

K_{m}

denotes the system matrix during emission. They are used to calculate the system weight matrix A.

S_{x}

is the excitation source distribution.

F

is obtained by discretizing the fluorescent yield distribution. Based on Eqs. (3) and (4), the FMT problem can be formulated as the following linear matrix equation:

Φ = A X

where

Φ

denotes the measurements of FMT, and

A

denotes the system weight matrix. X denotes the intensity of the fluorescence distribution in biological tissues [35]. Therefore, solving the FMT inverse problem is aimed at recovering the fluorescent distribution X in the aforementioned linear matrix equation.

2.2 Reconstruction based on l_2,1-norm optimization

As mentioned above, FMT is usually an ill-posed problem, which means that the dimension of the null space of matric A is not zero. That is, the solution of the problem is not unique in this situation. Even though the FMT problem may become less ill-posed when we capture more fluorescence measurement data sets are captured by the cooled CCD camera, it can also remain ill-conditioned. Therefore, the errors in the FMT problem may be large, which will affect the accuracy of the reconstruction results. Hence, Eq. (5) needs to be regularized in order to achieve an accurate and robust solution. We consider that a priori information encompasses the sparsity of the fluorescent sources. Although there are some standard structural methods (e.g. l₁ and smooth-l₀ norms), it is very difficult to solve the inverse problem due to the non-convexity and non-smoothness of such norms [23]. To overcome the limitation, structured sparsity theories have been used [24–26]. In addition to the sparsity of the structure, we exploit another priori information that the fluorescent sources have the group structured sparsity. Since the fluorescent sources are clustered together, we suppose that they have a group structure, where the components in the same group are all zeros or nonzeros. Considering the group structure of the fluorescent region, l_2,1-norm optimization was adopted in this paper. The optimization function with respect to Eq. (5) is formulated as follows:

\min_{X} E (X) = \frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + λ {‖ X ‖}_{2, 1}

where λ is the regularization parameter used to balance the two terms, and

{‖ X ‖}_{2, 1} = \sum_{i = 1}^{n} {‖ x_{i} ‖}_{2}

is the l_2,1-norm of matrix X.

Optimizing l_2,1-norm based problems in general is not an easy task, due to non-differentiability inherent in the mixed l_2,1-norm, where the objective function in Eq. (6) is not smooth. Thus, we propose to solve the equivalent convex smooth reformulations of Eq. (6). First, we introduce an additional variable $t = {[t_{1}, t_{2}, \dots, t_{3}]}^{T}$ , where $t_{i}$ acts as the upper-bound of $x_{i}$ . Because the $\frac{1}{2} {‖ A X - Φ ‖}_{2}^{2}$ term is a convex smooth loss function, then, the object function is equivalent to the following constrained convex smooth optimization problem:

\min_{(t, X) \in D} \frac{1}{2} {‖ A X - Φ ‖}_{2}^{2} + λ \sum_{i = 1}^{n} t_{i}

where

D = {(t, X) | ‖ x_{i} ‖ \leq t_{i}, \forall i = 1, 2, \dots, n}

.

We employ Nesterov’s method to solve Eq. (7). Nesterov’s method has a much faster convergence rate than the traditional methods such as subgradient descent and gradient descent. Nesterov’s method is based on two sequences ${s_{i}}$ and ${x_{i}}$ in which ${s_{i}}$ is the sequence of search points and ${x_{i}}$ is the sequence of approximate solutions. The search point $s_{i}$ is the affine combination of $x_{i - 1}$ and $x_{i}$ given by

s_{i} = x_{i} + α_{i} (x_{i} - x_{i - 1})

where

α_{i}

is the combination coefficient. The approximate solution

x_{i + 1}

is computed as a “gradient” step of

s_{i}

by

x_{i + 1} = π_{G} (s_{i} - \frac{1}{β_{i}} \nabla f (s_{i}))

where

f (x) = \frac{1}{2} {‖ A x - Φ ‖}_{2}^{2}

,

\nabla f (x) = A^{T} (A x - Φ)

,

\frac{1}{β_{i}}

is the stepsize, and

β_{i}

is determined by the line search depending on the Armijo-Goldstein rule so that

β_{i}

should be appropriate for

s_{i}

.

π_{G} (v)

, which is the Euclidean projection of

v

onto convex set G:

π_{G} (v) = \min_{x \in G} \frac{1}{2} {‖ x - v ‖}^{2}

The choice of

β

is key for the convergence analysis of Nesterov’s method. For a given

β > 0

, we define

f_{β, x} (y) = f (x) + 〈 \nabla f (x), y - x 〉 + \frac{β}{2} {‖ y - x ‖}_{2}^{2}

which is the tangent line of

f (\cdot)

at x, regularized by the square distance between x and y. It is easy to confirm that the quadratic function

f_{β, x} (y)

is strongly convex, and minimizing

f_{β, x} (y)

shows that domain G is equivalent to the problem of Euclidean projections onto G, which is,

π_{G} (x - \frac{1}{β} \nabla f (x)) = \arg \min_{y \in G} f_{β, x} (y)

We say that

β

is “appropriate” for

x

, if the following inequation holds:

f (π_{G} ((x - \frac{1}{γ_{0}} \nabla f (x))) \leq f_{β, x} (π_{G} ((x - \frac{1}{β} \nabla f (x)))

More detailed descriptions can be found in [36].

For a sequence of ${s_{i}}$ and ${x_{i}}$ , $s_{i}$ and $x_{i + 1}$ are computed recursively according to Eqs. (8) and (9) until the optimal solution $x^{*}$ is obtained. The illustration of Nesterov’s method is shown in Fig. 1, and the flowchart of the proposed method is illustrated in Table 1.

Fig. 1 The illustration of Nesterov’s method. We set $x_{1} = x_{0}$ , so $s_{1} = x_{1}$ . We assume that $x_{6} = x_{7} = x *$ , and $x *$ is an optimal solution.

Download Full Size | PDF

Table 1. Nesterov’s Method For L_2,1-Norm Optimization

View Table | View all tables in this article

3. Results

In this section, both numerical simulation experiments and in vivo mouse experiments are conducted to analyze the robustness, accuracy and efficiency of our proposed method. All of the reconstructions were carried out on our desktop computer with 4 GB RAM and 3.60 GHz Intel Core i7-4790 CPU.

In order to quantitatively evaluate the performance of the reconstruction results, the position error (PE) and contrast-to-noise ratio (CNR) are introduced in this paper. PE is defined to analyze the location error between the barycenter of the reconstruction region and real fluorescent region. The definition of PE is given by

P E = {‖ P_{r} - P_{0} ‖}_{2}

where

P_{0}

is the barycenter of the real fluorescent region and

P_{r}

is the barycenter of the reconstructed region.

CNR is introduced to indicate whether the reconstructed region could be clearly distinguished from the background.

C N R = \frac{μ_{V O I} - μ_{V O B}}{\sqrt{ω_{V O I} σ_{V O I}^{2} + ω_{V O B} σ_{V O B}^{2}}}

where

ω_{V O I}

,

ω_{V O B}

are the weight factors of the volume of interest (VOI) and volume of background (VOB) relative to the entire image volume, while

μ_{V O I}

,

μ_{V O B}

are the mean values of VOI and VOB, and

σ_{V O I}^{}

,

σ_{V O B}^{}

are the respective standard deviations. In this paper, the VOI is defined according to a threshold of 30% of the maximum fluorescent intensity.

3.1 Convergence analysis

To evaluate the convergence of our proposed method, a one-source mouse-mimicking heterogeneous phantom with a 20 mm height and 20 mm diameter was utilized, as shown in Fig. 2(a). Four kinds of materials were used to represent lungs (L), heart (H), bone (B) and muscle (M). The optical parameters for each kind of material are presented in Table 2, while the excitation and emission wavelength corresponding to the optical parameters are 780nm and 830nm, respectively [37]. The black dots in Fig. 2(b) denote different locations of the excitation light sources. A spherical fluorescent source with a 2 mm diameter centered in the z = 0 plane was placed in the right lungs, as shown in Fig. 2(b). The fluorescent yield of the sources was set to be 0.02 mm⁻¹. For each excitation point source, the emitted fluorescent measurements were captured from the opposite side of the heterogeneous cylindrical phantom within a 160° field of view (FOV), as shown in Fig. 2(b). To reconstruct the fluorescent sources, this phantom was discretized into 5657 nodes and 32670 tetrahedral elements.

Fig. 2 The one-source mouse-mimicking heterogeneous phantom, (a) is the 3D view of the fluorescent sources. (b) is the cross-section in the z = 0 plane. The black dots in (b) denote the excitation point sources. Fluorescence was collected from the opposite cylindrical surface within 160° FOV for each point source.

Download Full Size | PDF

Table 2. The Optical Parameters Of The Heterogeneous Phantom

View Table | View all tables in this article

Usually, the tomographic imaging quality is typically sensitive to the selection of the regularization parameter λ. Hence, λ with different orders of magnitude ranging from 10⁻¹ to 10⁻¹⁰ were tested. For each parameter λ, different iteration numbers, ranging from 50 to 500 were conducted to test the convergence of our proposed method. The reconstruction results in cross-sectional views are demonstrated in Figs. 3(a)-3(j). Furthermore, the corresponding numerical results are shown in Fig. 3(k). It can be observed that the errors of the reconstruction results are slightly changed when the iteration number is greater than 300. The errors are reduced with the decrease of the parameter. When parameter λ<10⁻⁶, the distribution shape of the source becomes deviated from the original shape. We therefore conclude from above that our proposed method is convergent when the iteration number is greater than 300, and the reconstruction results are satisfactory when parameter λ = 10⁻⁵.

Fig. 3 Reconstruction results with different regularization parameters: (a)-(j) are cross-sectional views in the z = 0 plane when iteration number is 300, and λ = 10⁻¹, 10⁻², 10⁻³, 10⁻⁴, 10⁻⁵, 10⁻⁶, 10⁻⁷, 10⁻⁸, 10⁻⁹ and 10⁻¹⁰ respectively. The red circles in the slice images represent the real locations of the fluorescent sources. (k) numerical analysis for the reconstruction errors with different iteration numbers and parameters.

Download Full Size | PDF

3.2 Results of the heterogeneous phantom experiment

To evaluate the performance of the proposed method, a two-source mouse-mimicking heterogeneous phantom was utilized, as shown in Fig. 4(a). Four kinds of materials were used to represent lungs (L), heart (H), bone (B) and muscle (M). The optical parameters for each kind of material are presented in Table 2. Similarly, to reconstruct the fluorescent sources, this phantom was discretized into 5657 nodes and 32670 tetrahedral elements.

Fig. 4 The two-source mouse-mimicking heterogeneous phantom, (a) is the 3D view of the fluorescent sources. (b) is the cross-section in the z = 0 plane. The black dots in (b) denote the excitation point sources. Fluorescence was collected from the opposite cylindrical surface within 160° FOV for each point source.

Download Full Size | PDF

To evaluate our proposed method, two other classical reconstruction methods were chosen to compared to reconstruct the same data set. One was the Tikhonov based method with the L2-norm (Tikhonov_L2) [38]. The other one was the iterated shrinkage based method with the L1-norm (IS_L1) [35]. To ensure there were sufficient convergences of the three methods, all the iteration numbers were set to 400, and the parameter of our proposed method was set to 10⁻⁵. The parameter of the Tikhonov_L2 method was set to 10⁻² and the IS_L1 method was set to 10⁻³, which were optimal choices according to the analysis above and reference [35], [38]. The reconstruction results of the three methods are shown in Fig. 5. The first row in Fig. 5 shows the 3-D reconstruction results corresponding to Tikhonov-L2, IS_L1 and l_2,1-norm respectively, while the second row illustrates the cross-sectional views corresponding to the three methods. The red spheres in the 3-D view and red circles in the cross-sectional views denoted the real positions and regions of the fluorescent sources. All of the cross-sectional views in Fig. 5 are shown in the z = 0 slice.

Fig. 5 Reconstruction results from the Tikhonov_L2 method (a), the IS_L1 method (b), and the proposed method (c) for two fluorescent sources and four measurement data sets. The first row illustrates the 3-D views of the reconstruction results. The second row illustrates the slice images in the $z = 0$ plane. The red circles in the slice images represent the real locations of the fluorescent source.

Download Full Size | PDF

The quantitative analysis of the performance of the three methods is listed in Table 3. Compared to the Tikhonov_L2 method, our proposed method results in much smaller PE values. The CNRs of our proposed method are higher than the Tikhonov_L2 method and the IS_L1 method, which implies that our proposed method is more effective in distinguishing the fluorescent signal from the background. The time consumptions of the three different methods show that our proposed method was faster than the other two methods. In addition, the fluorescence reconstruction results show that our proposed method achieves higher contrasts, accuracy and efficiency.

Table 3. Quantitative Analysis Of The Three Methods In Numerical Phantom Experiments

View Table | View all tables in this article

In order to assess the reconstruction performance in different organs, we varied the source location. One was located in the right lung, and the other one was located in the muscle, as shown in Fig. 6. The diameters of the two sources were both set as 2 mm, and the optical parameters for each kind of material are listed in Table 2. To reconstruct the fluorescent sources, this phantom was discretized into 5693 nodes and 34017 tetrahedral elements.

Fig. 6 The two-source mouse-mimicking heterogeneous phantom, source 1 located in the right lung and source 2 located in the muscle, (a) is the 3D view of the fluorescent sources. (b) is the cross-section in the z = 0 plane. The black dots in (b) denote the excitation point sources. Fluorescence was collected from the opposite cylindrical surface within 160° FOV for each point source.

Download Full Size | PDF

To evaluate the reconstruction performance, the proposed method was also compared with the Tikhonov_L2 method and the IS_L1 method. To ensure convergence of the three methods, the iteration number was set to 400 for all three methods. The parameter of our proposed method was set to 10⁻⁵. The parameter of the Tikhonov_L2 method and the IS_L1 method was set to 10⁻² and 10⁻³ respectively, which were the optimal choices according to the above analysis and references [35], [38]. The reconstruction results of the three methods are shown in Fig. 7. The first row in Fig. 7 shows the 3-D reconstruction results corresponding to Tikhonov-L2, IS_L1 and l_2,1-norm respectively, while the second row illustrates the cross-sectional views corresponding to the three methods. The red spheres in the 3-D view and red circles in the cross-sectional views denote the real positions and regions of the fluorescent sources. All of the cross-sectional views in Fig. 7 are shown in the z = 0 slice.

Fig. 7 Reconstruction results from the Tikhonov_L2 method (a), the IS_L1 method (b), and the proposed method (c) for two fluorescent sources and four measurement data sets. The first row illustrates the 3-D views of the reconstruction results. The second row illustrates the slice images in the z = 0 plane. The red circles in the slice images represent the real locations of the fluorescent source.

Download Full Size | PDF

The quantitative analysis of the performance of the three methods is listed in Table 4. Since different organs have different optical parameters, by using the IS_L1 method, we can only accurately reconstruct the one source located in the muscle. The other source results in over-shrinkage. Although we could obtain both sources with Tikhonov_L2 method, the PEs are larger than those of our proposed method. Hence, compared to the Tikhonov_L2 and the IS_L1 method, our proposed method gives better results in particular when the sources are located in different organs.

Table 4. Quantitative Analysis Of The Three Methods In Numerical Phantom Experiments with Different Source Locations

View Table | View all tables in this article

To further evaluate the robustness of the reconstruction performance, we conducted the test of noise intensity for our proposed method. Because noise is inevitable in the experiments, it is usually difficult to distinguish true signals from noise especially when the noise is large. As a result, the tomographic imaging quality is sensitive to noise. In this experiment, we used the heterogeneous phantom whose sources were both in the right lung. We set 4 groups of different noise intensities to simulate the worst case scenario. Due to the sensitivity to noise, the robustness of the FMT reconstruction is critical for preclinical in vivo application. Here, the measurement data sets were corrupted by 5%, 10%, 15%, and 20% Gaussian noise respectively.

The reconstruction results of the three methods with different noise intensities are shown in Fig. 8, and the quantitative analysis of the performance for the three methods is listed in Table 5. The columns in Fig. 8 denote the reconstruction results corresponding to different noise intensities. The first, third and fifth rows in Fig. 8 illustrated the 3-D reconstruction results corresponding to Tikhonov_L2, IS_L1 and l_2,1-norm, respectively. While the second, fourth and sixth rows show the cross-sectional views corresponding to the three methods. The red circles in the cross-sectional views denote the real positions and regions of the fluorescent sources. All of the cross-sectional views in Fig. 8 are shown in the z = 0 slice.

Fig. 8 Reconstructed results of heterogeneous numerical phantom experiments with different noise intensities. The columns denote the reconstructed results corresponding to different noise intensities (5%, 10%, 15%, 20%). The first and second rows are the reconstruction results obtained using the Tikhonov_L2 method. The third and fourth rows correspond to the IS_L1 method. The fifth and sixth rows correspond to the l_2,1-norm method. _3-D denotes the 3-D views of the reconstruction results, and _CV denotes the cross-sectional views of the results. The red spheres in the 3-D views and the red circles in the cross-sectional views show the real positions of the fluorescent sources.

Download Full Size | PDF

Table 5. Quantitative analysis of the three methods in numerical phantom experiments with different noise intensities.

View Table | View all tables in this article

Based on the experimental results presented above, we arrive at the following conclusions. (1) When the intensity of the noise is increased, our proposed method is more robust to noise. Even when the measurement data sets were corrupted by 20% Gaussian noise, our proposed method could achieve effective results compared to those obtained from the Tikhonov_L2 and IS_L1 methods. (2) For the same data sets and noise intensity, our proposed method achieves better performance in PEs for both targets. Besides, the fluorescence reconstructed by our proposed method has a higher contrast than the other two methods. Unlike the L1-regularization and L2-regularization method, the l_2,1-norm regularization method acted out with a stronger convergence property and fewer reconstruction artifacts. (3) For the same data sets, Tikhonov_L2 method takes approximately 100 seconds to complete the reconstruction, IS_L1 can considerably decrease the computation time to approximately 50 seconds. Due to the relatively simple equivalent form of l_2,1-norm regularization and accelerated techniques, it takes only 5 seconds to obtain satisfactory results using l_2,1-norm method. This shows that our proposed method produces significant improvement in computation time and efficiency.

3.3 Practical application of the in vivo mouse experiment

To evaluate the feasibility of our proposed method in in vivo applications of FMT, an in vivo mouse experiment was conducted. In this experiment, an adult Kunming mouse (Laboratory Animal Center, Peking University, China) was imaged based on a dual-modality micro-CT and FMT imaging system exploited previously by our group [39, 40]. The schematic illustration of the dual-modality micro-CT and FMT imaging system is presented in Fig. 9. The system mainly consists of a continuous wave laser (with a center wavelength of 671 nm), along with an ultrasensitive cooled CCD camera (PIXIS 1024BR, Princeton Instruments, USA), an x-ray generator (UltraBright, Oxford Instruments, USA), an x-ray detector (CMOS Flat-panel Detector, Hamamatsu, Japan), a set of optical lenses and a rotating stage.

Fig. 9 The schematic illustration of the dual-modality micro-CT and FMT imaging system.

Download Full Size | PDF

Before x-ray and optical data acquisition, the mouse was anesthetized first, and then the fluorescent bead (3mm in diameter) with the cy5.5 solution was implanted into the interspaces in the neighborhood of the liver stereotactically. The cy5.5 solution has an extinction coefficient of 0.019 ${mm}^{- 1} μ M^{- 1}$ and a quantum efficiency of 0.23, the peak excitation wavelength of the solution is 671 nm and the emission wavelength is 710 nm [41], The main process of the in vivo experiments can be summarized as follows. Firstly, the optical data were acquired. The detector FOV was set to be 160° and eight projections were adopted. The mouse was scanned using micro-CT after the acquisition of the optical data. Since the fluorescent bead was wrapped in a plastic material, it could be detected by the micro-CT system to locate the real fluorescent region easily. Then, the structural data of the mouse were reconstructed by the Feldkamp-Davis-Kress method [42]. As shown in Fig. 10, the green square marks the position of the fluorescent bead at the coordinates (43.20mm, 46.00mm, 6.40mm). In addition, the major organs including muscle, lungs, heart, liver and kidneys were segmented to construct the heterogeneous mouse model. The optical properties of the mouse organs were calculated according to Ref [43], which are listed in Table 6. Furthermore, in order to depict the photon distribution on the surface of the heterogeneous mouse torso, the fusion of the mesh and the fluorescence data were carried out via a 3-D surface flux reconstruction algorithm [44]. The heterogeneous mouse model and the fusion are shown in Fig. 11. Finally, the heterogeneous mouse model was discretized into a volumetric mesh with 5302 nodes and 29414 tetrahedral elements for the reconstruction of FMT.

Fig. 10 Anatomical structure of the mouse. (a) 3-D visualization of the mouse. (b) Transverse view. (c) Sagittal view. (d) Coronal view. The green square marker in (b), (c), (d) illustrates the location of the fluorescent bead.

Download Full Size | PDF

Table 6. Optical Properties Of The In Vivo Mouse Model

View Table | View all tables in this article

Fig. 11 Heterogeneous mouse model for the in vivo experiments. (a)Heterogeneous mouse torso for imaging reconstructions, including muscle, kidneys, liver, lungs and heart. (b) Surface view of the mouse torso with the frontal view measurements mapped on it.

Download Full Size | PDF

After the above procedures, the reconstruction of FMT was performed by using three different reconstruction methods (Tikhonov_L2, IS_L1 and l_2,1-norm). The 3-D reconstruction results of the three methods are shown in Fig. 12. The cross-sectional views of the results are demonstrated in Fig. 13. Quantitative analysis comparing PE, CNR and computation time of the three methods is shown in Table 7.

Fig. 12 3-D views of the results utilizing the Tikhonov_L2 method (a), the IS_L1 method (b), and the proposed method (c). The blue plane in the figure is the z = 6.4 mm slice.

Download Full Size | PDF

Fig. 13 Cross-sectional views of the in vivo region reconstruction results. (a), (d), and (e) are the transverse views of the reconstruction results via the Tikhonov_L2, IS_L1 and proposed method, respectively. The results are compared with the corresponding micro-CT slices. The cross on the red line denotes the real position of the fluorescent bead.

Download Full Size | PDF

Table 7. Optical Properties Of The In Vivo Mouse Model

View Table | View all tables in this article

The experimental results indicate that the reconstruction results achieved by our proposed l_2,1-norm method are better than the other two contrast methods. It can be observed that the fluorescent source reconstructed by the Tikhonov_L2 method is not accurately localized with a location error of 1.50 mm and it is widely scattered. The IS_L1 method and our proposed method are both able to obtain satisfactory reconstruction results. However, the results of our proposed method are better than those results of the IS_L1 method since the fluorescent source reconstructed by our proposed method has not only a higher accuracy, but also a higher contrast to the background. The reconstruction time of our proposed method is 7.45 seconds which is one order of magnitude less than the other two contrasting methods, verifying the absolute advantages of our proposed method in efficiency. This further demonstrates the improvement of our method on the efficiency and computation cost of the problem, and implies that our proposed method has potential to achieve more reliable fluorescent targets in biological applications.

4. Conclusion

In this study, a novel method based on the l_2,1-norm method has been proposed to reconstruct the internal fluorescent sources for FMT. It is well known that the ill-posedness of the FMT inverse problem can cause large errors in reconstruction. A variety of reconstruction algorithms have been proposed, such as the Gauss–Newton method, conjugate gradient method, interior-point method, etc [13–15]. to reconstruct the fluorescent distributions. To solve the ill-posedness problem of FMT, the sparsity of the fluorescent region was used as a priori knowledge in this paper. Besides, we exploited the property of fluorescent sources that fluorescent sources are clustered and that the image domain is zero everywhere except for the rigion where fluorescent sources are located. Based on this property, we assumed in our problem formulation that the fluorescent sources had a group structure [27, 28]. Hence, l_2,1-norm optimization based on group structured sparsity was utilized to solve the problem in this work. To accurately and effectively solve the l_2,1-norm model, Nesterov’s method was employed. Because the l_2,1-norm method combines the advantages of L1-regularization and L2-regularization. The reconstruction of the l_2,1-norm method was neither over-smoothed nor over-shrunk. Thus, l_2,1-norm is an effective reconstruction strategy for FMT, and it can reconstruct satisfactory fluorescent targets inside biological tissues.

To evaluate the performance of the proposed method, both numerical simulation experiments and in vivo mouse experiments are conducted. For comparison, the Tikhonov_L2 method and IS_L1 method are also used on the same data set to obtain reconstruction results. Based on the simulations and in vivo experiments, we arrive at the following conclusions. (1) Given the same number of measurement data sets, our proposed method gives higher reconstruction accuracy than the Tikhonov_L2 method and IS_L1 method. (2) Our method can enhance the robustness to noise and prevent artifacts in the background. (3) Our method proved to be more efficient compared to traditional methods. The computation time of our proposed method was about one order of magnitude less than that of traditional methods.

Although the proposed l_2,1-norm method achieved promising results as demonstrated in this paper, there are still some challenging problems in FMT, and other performances of the proposed method, such as spatial resolution and depth sensitivity, should also be tested systematically. It should be noted that the l_2,1-norm method is based on cluster assumption, but in cases where nanoparticles are injected systemically, the cluster assumption is no longer valid. Therefore, more effective and universal reconstruction methods need be developed in the future. It is believed that FMT will provide more potentials for earlier detection and characterization of the disease and evaluation of the treatment with rapid development of 3-D reconstruction algorithms [45].

Acknowledgments

This paper is supported by the National Natural Science Foundation of China under Grant Nos. 81571836, 81227901, 61231004 and 61501462, the Fundamental Research Funds for the Central Universities under Grant No. 2015JBM026, the China Postdoctoral Science Foundation under Grant Nos. 2015M570032 and 2015T80155, the Beijing Natural Science Foundation under Grant Nos. 7164270, 20140491524, and the Open Research Project under Grant 20150103 from SKLMCCS. The authors thank the International Innovation Team of CAS.

References and links

1. W. R. Zipfel, R. M. Williams, and W. W. Webb, “Nonlinear magic: multiphoton microscopy in the biosciences,” Nat. Biotechnol. 21(11), 1369–1377 (2003). [CrossRef] [PubMed]

2. D. J. Stephens and V. J. Allan, “Light microscopy techniques for live cell imaging,” Science 300(5616), 82–86 (2003). [CrossRef] [PubMed]

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

4. C. Qin, S. Zhu, and J. Tian, “New optical molecular imaging systems,” Curr. Pharm. Biotechnol. 11(6), 620 (2010). [CrossRef] [PubMed]

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

6. C. C. Leng and J. Tian, “Mathematical method in optical molecular imaging,” Sci. China Life Sci. 58, 1–13 (2015). [PubMed]

7. V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods 7(8), 603–614 (2010). [CrossRef] [PubMed]

8. X. Song, D. Wang, N. Chen, J. Bai, and H. Wang, “Reconstruction for free-space fluorescence tomography using a novel hybrid adaptive finite element algorithm,” Opt. Express 15(26), 18300–18317 (2007). [CrossRef] [PubMed]

9. J. Ye, Y. Du, Y. An, C. Chi, and J. Tian, “Reconstruction of fluorescence molecular tomography via a nonmonotone spectral projected gradient pursuit method,” J. Biomed. Opt. 19(12), 126013 (2014). [CrossRef] [PubMed]

10. Y. An, J. Liu, G. Zhang, J. Ye, Y. Du, Y. Mao, C. Chi, and J. Tian, “A Novel Region Reconstruction Method for Fluorescence Molecular Tomography,” IEEE Trans. Biomed. Eng. 62(7), 1818–1826 (2015). [CrossRef] [PubMed]

11. D. Wang, X. Song, and J. Bai, “Adaptive-mesh-based algorithm for fluorescence molecular tomography using an analytical solution,” Opt. Express 15(15), 9722–9730 (2007). [CrossRef] [PubMed]

12. W. Bangerth and A. Joshi, “Adaptive finite element methods for the solution of inverse problems in optical tomography,” Inverse Probl. 24(3), 034011 (2008). [CrossRef]

13. R. Noumeir, G. E. Mailloux, H. Mallouche, and R. Lemieux, “Comparison between ML-EM and modified Newton algorithms for SPECT image reconstruction,” in OE/LASE'93: Optics, Electro-Optics, & Laser Applications in Science& Engineering, (Academic, 1993).

14. T. Tarvainen, M. Vauhkonen, and S. R. Arridge, “Gauss–Newton reconstruction method for optical tomography using the finite element solution of the radiative transfer equation,” J. Quant. Spectrosc. Ra. 43(4), 252–261 (2010).

15. W. Zhu, Y. Wang, Y. Yao, J. Chang, H. L. Graber, and R. L. Barbour, “Iterative total least-squares image reconstruction algorithm for optical tomography by the conjugate gradient method,” J. Opt. Soc. Am. A 14(4), 799–807 (1997). [CrossRef] [PubMed]

16. N. Cao, A. Nehorai, and M. Jacobs, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007). [CrossRef] [PubMed]

17. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). [CrossRef] [PubMed]

18. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]

19. Y. Lu, X. Zhang, A. Douraghy, D. Stout, J. Tian, T. F. Chan, and A. F. Chatziioannou, “Source reconstruction for spectrally-resolved bioluminescence tomography with sparse a priori information,” Opt. Express 17(10), 8062–8080 (2009). [CrossRef] [PubMed]

20. M. Süzen, A. Giannoula, and T. Durduran, “Compressed sensing in diffuse optical tomography,” Opt. Express 18(23), 23676–23690 (2010). [CrossRef] [PubMed]

21. J. Prakash, C. B. Shaw, R. Manjappa, R. Kanhirodan, and P. K. Yalavarthy, “Sparse Recovery Methods Hold Promise for Diffuse Optical Tomographic Image Reconstruction,” IEEE J. Sel. Top. Quantum Electron. 20(2), 74–82 (2014). [CrossRef]

22. S. Okawa, Y. Hoshi, and Y. Yamada, “Improvement of image quality of time-domain diffuse optical tomography with l sparsity regularization,” Biomed. Opt. Express 2(12), 3334–3348 (2011). [CrossRef] [PubMed]

23. C. Chen, F. Tian, H. Liu, and J. Huang, “Diffuse optical tomography enhanced by clustered sparsity for functional brain imaging,” IEEE Trans. Med. Imaging 33(12), 2323–2331 (2014). [CrossRef] [PubMed]

24. J. Huang, T. Zhang, and D. Metaxas, “Learning with structured sparsity,” J. Mach. Learn. Res. 13(7), 3371–3412 (2011).

25. R. Baraniuk, V. Cevher, M. Duarte, and C. Hegde, “Model-based compressive sensing,” IEEE Trans. Inf. Theory 56(4), 1982–2001 (2010). [CrossRef]

26. V. Cevher, P. Indyk, C. Hegde, and R. Baraniuk, “Recovery of clustered sparse signals from compressive measurements,” presented at the Int. Conf. Sampl. Theory Appl., Marseille, France (2009).

27. M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” J. R. Stat. Soc. Ser. A Stat. Soc. 68(1), 49–67 (2006). [CrossRef]

28. J. Huang and T. Zhang, “The Benefit Of Group Sparsity,” Ann. Stat. 38(4), 1978–2004 (2010). [CrossRef] [PubMed]

29. L. Jacob, G. Obozinski, and J. P. Vert, “Group lasso with overlap and graph lasso,” in International Conference on Machine Learning (ACM, 2009), pp. 433–440. [CrossRef]

30. A. Nemirovski, Efficient methods in convex programming (Academic, 2005).

31. Y. Nesterov, Introductory lectures on convex optimization (Springer Science & Business Media, 2004).

32. J. H. Lee, A. Joshi, and E. M. Sevick-Muraca, “Fully adaptive finite element based tomography using tetrahedral dual-meshing for fluorescence enhanced optical imaging in tissue,” Opt. Express 15(11), 6955–6975 (2007). [CrossRef] [PubMed]

33. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. Wang, E. Hoffman, G. McLennan, P. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [CrossRef] [PubMed]

34. A. Joshi, W. Bangerth, and E. Sevick-Muraca, “Adaptive finite element based tomography for fluorescence optical imaging in tissue,” Opt. Express 12(22), 5402–5417 (2004). [CrossRef] [PubMed]

35. D. Han, J. Tian, S. Zhu, J. Feng, C. Qin, B. Zhang, and X. Yang, “A fast reconstruction algorithm for fluorescence molecular tomography with sparsity regularization,” Opt. Express 18(8), 8630–8646 (2010). [CrossRef] [PubMed]

36. J. Liu, S. Ji, and J. Ye, “Multi-Task Feature Learning Via Efficient L2,1-Norm Minimization,” Eprint Arxiv339–348 (2009).

37. A. B. Milstein, J. J. Stott, S. Oh, D. A. Boas, R. P. Millane, C. A. Bouman, and K. J. Webb, “Fluorescence optical diffusion tomography using multiple-frequency data,” J. Opt. Soc. Am. A 21(6), 1035–1049 (2004). [CrossRef] [PubMed]

38. W. Bangerth, “A Framework for the Adaptive Finite Element Solution of Large-Scale Inverse Problems,” SIAM J. Sci. Comput. 30(6), 2965–2989 (2008). [CrossRef]

39. S. Zhu, J. Tian, G. Yan, C. Qin, and J. Feng, “Cone beam micro-CT system for small animal imaging and performance evaluation,” Int. J. Biomed. Imaging 2009, 960573 (2009). [CrossRef] [PubMed]

40. Ping Wu, Yifang Hu, Kun Wang, and Jie Tian, “Bioluminescence Tomography by an Iterative Reweighted (l)2-Norm Optimization,” IEEE Trans. Biomed. Eng. 61(1), 189–196 (2014). [CrossRef] [PubMed]

41. F. Gao, H. Zhao, Y. Tanikawa, and Y. Yamada, “A linear, featured-data scheme for image reconstruction in time-domain fluorescence molecular tomography,” Opt. Express 14(16), 7109–7124 (2006). [CrossRef] [PubMed]

42. G. Yan, J. Tian, S. Zhu, Y. Dai, and C. Qin, “Fast cone-beam CT image reconstruction using GPU hardware,” J. XRay Sci. Technol. 16(4), 225–234 (2008).

43. 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(17), 4225–4241 (2005). [CrossRef] [PubMed]

44. X. Chen, X. Gao, D. Chen, X. Ma, X. Zhao, M. Shen, X. Li, X. Qu, J. Liang, J. Ripoll, and J. Tian, “3D reconstruction of light flux distribution on arbitrary surfaces from 2D multi-photographic images,” Opt. Express 18(19), 19876–19893 (2010). [CrossRef] [PubMed]

45. C. Chi, Y. Du, J. Ye, D. Kou, J. Qiu, J. Wang, J. Tian, and X. Chen, “Intraoperative imaging-guided cancer surgery: from current fluorescence molecular imaging methods to future multi-modality imaging technology,” Theranostics 4(11), 1072–1084 (2014). [CrossRef] [PubMed]

Input: $f (x)$ , $\nabla f (x)$ , $x_{0}$ , $β_{0}$
Output: $x$
1.	Initialize $x_{1} = x_{0} = 0$ , $β_{0} = 1$ , $t_{- 1} = 0$ , $t_{0} = 1$
2.	for $i = 1$ to $n$ do
3.	set $α_{i} = \frac{t_{i - 2} - 1}{t_{i - 1}}$ , $s_{i} = x_{i} + α_{i} (x_{i} - x_{i - 1})$
4.	for $j = 0$ to $k$ do
5.	set $β = 2^{j} β_{i - 1}$
6.	compute $x_{i + 1} = π_{G} (s_{i} - \frac{1}{β} \nabla f (s_{i}))$
7.	if $f (x_{i + 1}) \leq f_{β, x} (x_{i + 1})$
8.	$β_{i} = β$ , break
9.	end if
10.	end for
11.	set $t_{i} = (1 + \sqrt{1 + 4 t_{i - 1}^{2}}) / 2$
12.	end for
13.	Set $x = x_{i}$

Material	$μ_{a x}$ (m⁻¹)	${μ^{'}}_{s x}$ (m⁻¹)	$μ_{a m}$ (m⁻¹)	${μ^{'}}_{s m}$ (m⁻¹)
Muscle	5.2	1080	6.8	1030
Lungs	13.3	1970	20.3	1950
Heart	8.3	1010	10.4	990
Bone	2.4	1750	3.5	1610

Method	PE(mm)	CNR	Time(s)
Tikhonov_L2	S1:0.13 S2:0.18	30.09	101.07
IS_L1	S1:0.29 S2:0.56	41.32	50.25
l_2,1-norm	S1:0.13 S2:0.18	109.3	4.54

Method	PE(mm)	CNR	Time(s)
Tikhonov_L2	S1:0.65 S2:0.30	19.48	99.45
IS_L1	S1:- - S2:0.15	64.99	48.23
l_2,1-norm	S1:0.14 S2:0.16	89.69	4.37

Method	Quantitative Index	noise intensities
Method	Quantitative Index	5%	10%	15%	20%
Tikhonov _L2	PE(mm)	S1: 0.14 S2: 0.18	S1: 0.19 S2: 0.26	S1: 0.81 S2: 0.89	S1:– S2:–
	CNR	24.70	14.30	8.52	9.17
	Time(s)	97.96	96.55	95.85	92.54
IS_L1	PE(mm)	S1: 0.47 S2: 0.61	S1: 0.55 S2: 1.47	S1: 0.86 S2: 1.70	S1: 1.40 S2: 1.70
	CNR	53.54	44.17	32.02	39.48
	Time(s)	49.66	53.48	50.97	51.21
l_2,1-norm	PE(mm)	S1: 0.13 S2: 0.50	S1: 0.50 S2: 1.00	S1: 1.00 S2: 1.01	S1: 0.98 S2: 1.00
	CNR	84.67	83.47	121.24	102.93
	Time(s)	4.41	4.61	4.47	4.45

Novel l_2,1-norm optimization method for fluorescence molecular tomography reconstruction

Abstract

1. Introduction

2. Method

2.1 Photon propagation model

2.2 Reconstruction based on l_2,1-norm optimization

3. Results

3.1 Convergence analysis

3.2 Results of the heterogeneous phantom experiment

3.3 Practical application of the in vivo mouse experiment

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (13)

Tables (7)

Equations (15)

Biomedical Optics Express

Material	$μ_{a x} (m^{- 1})$	$μ_{s x} (m^{- 1})$	$μ_{a m} (m^{- 1})$	$μ_{a x} (m^{- 1})$
Liver	343.67	677	228.31	648
Lungs	191.82	2172	126.55	2124
Heart	57.45	962	38.29	905
Kidneys	64.45	2248	43.03	2109
Muscle	84.5	427.3	56.3	379.2

Method	PE(mm)	CNR	Time(s)
Tikhonov_L2	1.50	12.43	82.43
IS_L1	0.85	87.91	54.67
l_2,1-norm	0.54	191.22	7.45

Abstract

1. Introduction

2. Method

2.1 Photon propagation model

2.2 Reconstruction based on l2,1-norm optimization

3. Results

3.1 Convergence analysis

3.2 Results of the heterogeneous phantom experiment

3.3 Practical application of the in vivo mouse experiment

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (13)

Tables (7)

Equations (15)

Biomedical Optics Express

2.2 Reconstruction based on l_2,1-norm optimization