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

: Fluorescence molecular tomography (FMT) is a promising tomographic method in preclinical research, which enables noninvasive real-time three-dimensional (3-D) visualization for in vivo studies. The ill-posedness of the FMT reconstruction problem is one of the many challenges in the studies of FMT. In this paper, we propose a l 2,1 -norm optimization method using a priori information, mainly the structured sparsity of the fluorescent regions for FMT reconstruction. Compared to standard sparsity methods, the structured sparsity methods are often superior in reconstruction accuracy since the structured sparsity utilizes correlations or structures of the reconstructed image. To solve the problem effectively, the Nesterov’s method was used to accelerate the computation. To evaluate the performance of the proposed l 2,1 -norm method, numerical phantom experiments and in vivo mouse experiments are conducted. The results show that the proposed method not only achieves accurate and desirable fluorescent source reconstruction, but also demonstrates enhanced robustness to noise.


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][4][5][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][9][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 2 -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][14][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 2 -norm regularization, many recent methods have been proposed to exploit the sparsity of the reconstructed image [20][21][22], such as l 1 and smooth-l 0 norms. While these methods have been proved superior in overcoming the over-smoothing limitation of the traditional l 2 -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 2 -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,1norm 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.

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]: where the subscripts x and m denote the excitation and emission wavelength, respectively; Ω is the imaging domain of the problem; where 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: where X denotes the unknown fluorescent source vector to be reconstructed. x K denotes the system matrix during excitation and m K denotes the system matrix during emission. They are used to calculate the system weight matrix A. x S 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: 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.

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 1 and smooth-l 0 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][25][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: ( ) where λ is the regularization parameter used to balance the two terms, and Optimizing l 2,1 -norm based problems in general is not an easy task, due to nondifferentiability 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
is a convex smooth loss function, then, the object function is equivalent to the following constrained convex smooth optimization problem: where i α is the combination coefficient. The approximate solution 1 x i + is computed as a "gradient" step of i s by ( ) , which is the Euclidean projection of v onto convex set G: The choice of β is key for the convergence analysis of Nesterov's method. For a given 0 which is the tangent line of ( ) f ⋅ at x, regularized by the square distance between x and y. It is easy to confirm that the quadratic function , ( ) shows that domain G is equivalent to the problem of Euclidean projections onto G, which is, We say that β is "appropriate" for x , if the following inequation holds: More detailed descriptions can be found in [36]. (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.
We assume that 6 7 * x x x = = , and * x is an optimal solution.

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 where 0 P is the barycenter of the real fluorescent region and r P is the barycenter of the reconstructed region.
CNR is introduced to indicate whether the reconstructed region could be clearly distinguished from the background.
where VOI ω , VOB ω are the weight factors of the volume of interest (VOI) and volume of background (VOB) relative to the entire image volume, while VOI μ , VOB μ are the mean values of VOI and VOB, and VOI σ , VOB σ are the respective standard deviations. In this paper, the VOI is defined according to a threshold of 30% of the maximum fluorescent intensity.

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 −1 . 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.  Usually, the tomographic imaging quality is typically sensitive to the selection of the regularization parameter λ. Hence, λ with different orders of magnitude ranging from 10 −1 to 10 −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 −6 , 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 −5 .

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. 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 −5 . The parameter of the Tikhonov_L2 method was set to 10 −2 and the IS_L1 method was set to 10 −3 , 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.
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.   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. 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 −5 . The parameter of the Tikhonov_L2 method and the IS_L1 method was set to 10 −2 and 10 −3 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. 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. 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. 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.

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

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][14][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,1norm 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].