Open Access
Add to BibSonomy
Abstract
Fluorescent molecular tomography (FMT) is a highly sensitive and noninvasive imaging approach for providing three-dimensional distribution of fluorescent marker probes. However, owing to its light scattering effect and the ill-posedness of inverse problems, it is challenging to develop an efficient reconstruction algorithm that can achieve the exact location and morphology of the fluorescence source. In this study, therefore, in order to satisfy the need for early tumor detection and improve the sparsity of solution, we proposed a novel L1-L2 norm regularization via the forward-backward splitting method for enhancing the FMT reconstruction accuracy and the robustness. By fully considering the highly coherent nature of the system matrix of FMT, it operates by splitting the objective to be minimized into simpler functions, which are dealt with individually to obtain a sparser solution. An analytic solution of L1-L2 norm proximal operators and a forward-backward splitting algorithm were employed to efficiently solve the nonconvex L1-L2 norm minimization problem. Numerical simulations and an in-vivo glioma mouse model experiment were conducted to evaluate the performance of our algorithm. The comparative results of these experiments demonstrated that the proposed algorithm obtained superior reconstruction performance in terms of spatial location, dual-source resolution, and in-vivo practicability. It was believed that this study would promote the preclinical and clinical applications of FMT in early tumor detection.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
With the development of medical imaging technology, fluorescence molecular imaging (FMI) has widely applied in the study of drug delivery, protein function, preclinical, clinical, etc. [1–5]. However, FMI can only obtain the photon intensity of the object surface because of the absorption and scattering of light transmission. Therefore, fluorescence molecular tomography (FMT), which can accurately visualize and quantify the three-dimensional (3D) distribution of fluorescent sources in deep turbidities by capturing the surface fluorescence distribution of fluorescence sources using an ultra-sensitive charge-coupled device (CCD) camera, is favored by researchers [6–8].
However, owing to the strong scattering property of biological tissues and the limited boundary measurements with noise, FMT reconstruction is always an ill-posed problem [9,10]. Moreover, the mutual coherence of system matrix derived by finite element method (FEM) is always as high as above 90%, which also increase the complexity of FMT reconstruction [11]. Thus, the reconstruction is still a challenging problem [12–14]. Over the past few years, many reconstruction methods have been developed to reduce the ill-posedness [15–17]. Tikhonov regularization have widely used to overcome ill-posedness by adding $L_{2}$ norm [18]. The primary benefit of using $L_{2}$ norm is the simplicity of the optimization problem. However, it causes over-smoothing, leading to a blurring of the sources and making it difficult to obtain sharp boundaries [19].
FMT is usually utilized for early tumor imaging, and in this stage the tumor is always small and sparse in size compared to the entire body of the subject [20,21]. Mathematically speaking, the $L_{0}$ norm is the sparsest constrain. The reconstruction algorithm based on the $L_{0}$ norm regularization algorithm involves a problem of combinatory optimization, which makes it inefficient or unfeasible for practical applications. Thus, to improve the quality of the reconstruction, convex $L_{1}$ norm regularization has been adopted in FMT reconstruction. Several algorithms based on the $L_{1}$ norm have been proposed for optical molecular imaging and allow for more accurate results and less computation time, such as the iterated shrinkage algorithm with the $L_{1}$ norm [22], fast iterative soft-thresholding algorithm [23] , $L_{1}$ norm regularization algorithm based on the split Bregman algorithm [24,25], gradient-based techniques, and the incomplete variable truncated conjugate gradient algorithm. Moreover, the algorithm based on the $L_{1}$ norm can obtain satisfactory results and overcome the over-smoothing cause by the $L_{2}$ norm.
Recently, some novel regularization algorithm base on non-convex $L_{p}$ norm $(0<p<1)$ have been proposed. For example, He et al. proposed the half thresholding pursuit algorithm, with significant improvements in achieving the accuracy and sparse for the reconstruction of FMT [26]. In addition, Guo et al. proposed a non-convex sparse regularization approach (nCSRA) framework, which have a better multiple-source resolution capability in terms of position resolution, intensity resolution and size resolution [27]. Compared with the $L_{1}$ norm, the $L_{p}$ norm $(0<p<1)$ is a better approximation of the $L_{0}$ norm. Thus, an algorithm based on $L_{p}$ norm can always obtain a sparser solution compared with $L_{1}$ norm algorithms, thus ensuring stability at the same time.
Recently, it was demonstrated in a series of papers that the difference of the $L_{1}$ norm and $L_{2}$ norm ($L_{1}$-$L_{2}$ norm), outperforms $L_{1}$ norm and $L_{p}$ norm in terms of promoting sparsity when system matrix is highly coherent [28]. Inspired by this, $L_{1}$-$L_{2}$ norm regulation via difference of convex algorithm ($L_{1}$-$L_{2}$ via DCA) algorithm, has been used in FMT reconstruction having a great improvement in the sparsity of the solution [29]. However, the $L_{1}$-$L_{2}$ via DCA algorithm decomposes the objective function as the difference of two convex functions, convert to $L_{1}$ norm minimization problem, which may not have analytical solutions. Thus, the robustness of $L_{1}$-$L_{2}$ via DCA algorithm is therefore poor, especially for multiple sources resolution. Furthermore, because the $L_{1}$-$L_{2}$ via DCA algorithm for $L_{1}$-$L_{2}$ amounts to solving an $L_{1}$ norm minimization problem iteratively as a subproblem, it is much slower than $L_{1}$ norm minimization [30].In recently years, forward-backward splitting (FBS) method has been proposed to solve the nonlinear and nonconvex least-squares problem with tight frame sparsity regularization in quantitative photoacoustic tomography [31] and Low-dose CT [32]. These research prove the convergence and efficiency of FBS, and provide a proximal operator of the $L_{1}$-$L_{2}$ metric.
To satisfy the need for early tumor detection in FMT, a novel $L_{1}$-$L_{2}$ norm regularization algorithm via forward-backward splitting ($L_{1}$-$L_{2}$ via FBS) algorithm was proposed to recover the distribution of small fluorescent sources in orthotopic glioma mouse models. First, $L_{1}$-$L_{2}$ via FBS algorithm employ proximity operators and gradient descent to obtain the search point. It is worth mentioning that we derive an analytical solution for the proximal operator of the $L_{1}$-$L_{2}$ norm. Second, they operate by splitting the objective to be minimized into simpler functions that are dealt with individually to get a sparser solution. Thus, the main advantages of $L_{1}$-$L_{2}$ via FBS algorithm is full considering the highly coherent of the system matrix of FMT and the handling of non-differentiable objectives and non-convex constraints due to the $L_{1}$-$L_{2}$ norm while maintaining the simplicity of gradient descent methods and the sparsity of results. To validate the performance of the proposed algorithm, groups of comparison experiments were designed using numerical simulations and in-vivo imaging experiments for the $L_{1}$-$L_{2}$ via DCA algorithm [29], the incomplete variables truncated conjugate gradient algorithm (IVTCG) based on $L_{1}$ norm [33], and $L_{1/2}$ norm regulation based on Iterative shrinkage-thresholding algorithm (ISTA-$L_{1/2}$) [27]. It has experimentally shown that the proposed $L_{1}$-$L_{2}$ via FBS algorithm performance better in spatial location, dual-source resolution, robustness, and in-vivo practicability.
2. Materials and methods
2.1 Photon propagation model
In biological bodies, photon propagation within the near-infrared spectral band has a highly scattering feature. For steady-state FMT with point excitation sources, the following coupled diffusion equations (DE) have been commonly used to model the forward problem of FMT, which is defined as [34–36]:
To solve these equations, the Robin-type boundary conditions are added on the boundary $\partial \Omega$ of the domain $\Omega$ [37]:
where $\textbf {n}$ represents the outward normal vector to the surface. $q$ is a constant depending on the optical reflective index mismatch at the boundary.Based on FEM method, we can obtain the following linear model:
where $\textbf {A}$ is an $M \times N$ weighting matrix, and $\textbf {x}$ is an $N \times 1$ vector denotes the unknown internal distribution of the probes, $\pmb {\Phi }$ is an $M \times 1$ vector represents the measurements of surface photon distribution.2.2 Reconstruction Based on the forward-backward splitting (FBS) algorithm
In order to alleviate the ill-posedness of the inverse problem, the $L_{1}$-$L_{2}$ norm regularization was adopted in FMT inverse problem, writing as:
Obviously, $f(\textbf {x})$ is differentiable convex function; however $g(\textbf {x})$ is no-convex function and nondifferentiable, Eq. (5) cannot be minimized using simple gradient descent methods. For this, a forward-backward splitting (FBS) method was adopted to address each term in Eq. (5).
In the first stage, for improving the computational speed, the search point ${\textbf {s}_{k}}$, was determined by forward splitting, which can be obtained as follows:
where ${f{'}}(\textbf {x}) = \textbf {A}(\textbf {A}\textbf {x} - \pmb {\Phi })$ denotes the gradient of $f(\textbf {x})$ at the point $\textbf {x}$; and $t_{k}$ is the step size. and utilize the Lipschitz condition to get the value of ${t_k}$ [38].In the second stage, to ensure global convergence, an approximate solution $\textbf {x}$, and without straying too far from a starting point ${{\textbf {s}}_k}$. It can be through the proximal operator of $\textrm {g}(\textbf {x})$ obtained:
For any satisfying Eq. (9), we can obtain ${{\textbf {s}}_k} = \sigma {\textbf {q}} + (1 - \frac {\sigma }{{||{\textbf {x}}|{|_2}}}){\textbf {x}}$. And through a simple calculation, we have:
Therefore, the ${{\textbf {x}}^*}$ can be determined among all ${\textbf {x}}$ satisfying Eq. (10) with the largest norm. Without loss of generality ${{\textbf {s}}_k}^i$ is a non-increasing vector, ${{\textbf {s}}_k}^1 \ge {{\textbf {s}}_k}^2 \ge \cdots \ge {{\textbf {s}}_k}^i \ge 0$, then we discuss the optimal solution ${{\textbf {x}}^*}$:
1) When ${{\textbf {s}}_k}^1 > \sigma$, then ${{\textbf {s}}_k}^1 - \sigma {\textbf {q}_1} > 0$. For $i = 1$, the right-hand side (RHS) of Eq. (10) is positive, so the left-hand side (LHS) is positive as well; owing to the non-negative constraint on the unknown vector, we obtain $(1 - \frac {\sigma }{{||{\textbf {x}}|{|_2}}}) > 0$. For this $i$, which ${{\textbf {s}}_k}^i \le \sigma$, we have ${\textbf {x}}_i^* = 0$; Otherwise for this $i$, the LHS of Eq. (10) is positive, while the RHS is nonpositive. For this $i$ which ${{\textbf {s}}_k}^i > \sigma$, we have ${{\textbf {q}}_i} = 1$. Thus, for RHS of Eq. (10):${{{\textbf {s}}_k}} - \sigma {\textbf {q}} = {soft}({{\textbf {s}}_k},\sigma )$; ${soft}({{\textbf {s}}_k},\sigma )$ is the proximal operator of ${L_1}$ norm regularized, appoint $\textbf {z} = {soft}({{\textbf {s}}_k},\sigma )$, and defined as Eq. (13).
Do a simple operation on Eq. (10):
Based on the above statement, the pseudo code of the main steps of $L_{1}$-$L_{2}$ via FBS algorithm is given in Algorithm 1.
Algorithm 1. $L_{1}$-$L_{2}$ via FBS algorithm
3. Experiments setting
3.1 Numerical simulations
A series of numerical simulations based on the head of the common digital mouse model were conducted to evaluate the performance of $L_{1}$-$L_{2}$ via FBS algorithm, as shown in Fig. 1(a). The mouse head was segmented into three organs: the muscle, skull, and brain. Their related optical properties list in Table 1, with an excited light wave length is $750nm$ and an emission wavelength is $820nm$ [10]. Figure 1(b) shows the distribution of four excitation source, which were located one transport mean free path beneath the surface on the plane of $z=18mm$.
Fig. 1. (a) 3D views the head of a digital mouse. (b) Excitation sources on the plane of $z=18mm$. The four black circles denote the point excitation sources. (c) The views of simulation models with the single source. (d)-(f) The views of simulation models with the double source, where the edge-to-edge distance between the two source is different. (g) The mesh for inverse reconstruction.
Table 1. Optical parameters in numerical simulation and in-vivo experiments
In order to verify the feasibility of the $L_{1}$-$L_{2}$ via FBS algorithm, single and double source experiments were designed in feasible experiment. In the single source experiment, one sphere fluorescent source with a $1mm$ radius placed in the brain with center at $(9,13,18mm)$, as shown in Fig. 1(c). In double sources experiment, two sphere fluorescent sources with a $1mm$ radius are placed in the center positions placed at $(6,13,18mm)$ and $(12,13,18mm)$ respectively, shown in Fig. 1(d).
In addition, to verify the robustness of our proposed algorithm, two groups of numerical simulations were conducted by considering the number of excitation sources and the edge-to-edge distance between two sphere fluorescent sources in robustness experiment. Firstly, in the single source experiments, the number of the excitation sources was decrease to two (1 and 3 was used) and one (1 was used), respectively. Secondly, in the double source experiments, two sphere fluorescent sources with a $1mm$ radius were placed at $(6,13,18mm)$ and $(11,13,18mm)$, the edge-to-edge distance (EED) was $2mm$, as shown in Fig. 1(e). Furthermore, the center positions of two sphere fluorescent sources were placed at $(7.5,13,18mm)$ and $(10.5,13,18mm)$ with the edge-to-edge distance was $1mm$, as shown in Fig. 1(f).
In all numerical simulation experiment,the fluorescent yield $\eta {\mu _{af}}(r)$ of the fluorescent target was set to be 0.05 $m{m^{ - 1}}$.The digital mouse model was discretized into a tetrahedral mesh through Amria 5.2. A uniform tetrahedral mesh which include 10158 nodes and 52696 tetrahedral elements was used in inverse reconstruction. In addition, we used diffusion approximation model based on FEM to simulate the propagation of light through tissues to obtain specific surface energy distributions. The all-experiment codes were written in MATLAB2020 and were performed on a desktop computer with 2.90 GHz Intel Processor I5-9400F and 16G RAM
3.2 FMT imaging system
In order to acquire the CT data and optical image, an FMT imaging system was adopted in our study. Which consist of three parts: (1) Optical acquisition module mainly was composed of an electron multiplying charge coupled device (EMCCD) camera (iXon Ultra DU-897, Andor), a continuous wave optical laser (BS750-6WGA, BoSion Optoelectronic Technology, CHN) and optical filter (FB820-10, THORLABS, USA). (2) Micro-CT system included a X-tray detector (1512N-C90-HRCC, Dexela, USA), and an X-ray generator (L9181-02 MT2195, HAMAMATSUPHOTONICS, CHN); (3) Control module was consist of a 360-degree motorized turntable (RAK100, Zolix, CHN), and a controller (Zolix Instruments Co., Beijing, CHN). As show in Fig. 2, X-Ray Source, Turntable, and X-Ray Detector placed on the same line. In addition, The EMCCD was perpendicular with this line.
3.3 In-vivo imaging experiment
To further evaluate the performance of $L_{1}$-$L_{2}$ via FBS algorithm, in-vivo glioma mouse experiments were performed. Animal experiments were carried out under isoflurane gas anaesthesia, and every effort was made to relieve the pain of the male mice. Additionally, $5 \times 10^5$U-87 MG cells in $6\mu l$ phosphate buffer solution was injected into the brain of the mouse to construct the orthotopic glioma model. After seven days, the tumor-bearing mice were injected with Tf-IRDye800 (excitation spectrum: $750nm$, emission spectrum: $820nm$ ) via tail vein. Six hours later, fluorescence images and CT data were acquired. And the MRI data obtained subsequently, which was used to determine the location of the tumor.
To provide the excitation illumination, a $750nm$ continuous wave semiconductor laser with the output power of $450mW$ was used was used as the single excited source to reflex excite a fluorescent probe. And an electron multiplying charge coupled device (EMCCD) camera ((iXon Ultra DU-897, Andor) coupled with $820 \pm$ $20nm$ bandpass filter was cooled to $-80^\circ$C to collect surface fluorescence image with a $120^\circ$ field of view (FOV). The exposure time ($1.6s$), high shift speed ($12.9us$) and high-speed readout rate ($10MHz$ at $16$-$bit$) were set up to restrain the readout noise and the shifting noise. After that, X-ray source was set $70kVp$ and $39W$ to obtaining $360^\circ$ projection data.
Prior to the FMT reconstruction, some essential preprocessing operations were carried out. First, the raw CT data was converted into 3D volume data via the Feldkamp-Davis-Kress (FDK) algorithm for obtaining anatomical structural information. Second, the major organs including muscle, skull and brain were segmented by using Amria 5.2 (Amria, Visage Imaging, Australia). Lastly, the mouse model was discretized into 11504 nodes and 58970 tetrahedral elements for 2D-3D surface energy mapping and FMT reconstruction.
3.4 Algorithm comparison and evaluation index
For comparison, three effective algorithms based on sparse regularization (the $L_{1}$-$L_{2}$ via DCA, IVTCG, and the ISTA-$L_{1/2}$) were chosen as comparison methods. To promise all algorithms convergent, by experience the maximum iteration numbers and the minimum residuals were set as 1000 and 1e-6, respectively. In addition, to ensure the reliability of the results, the Generalized Cross-Validation method was adopted to determine the regulation parameters $\lambda$ [39]. The final regularization parameter values for different methods in different experiments have been listed in Table 2.
Table 2. Regularization parameters for different methods in different experiments.
To further quantitative evaluate the accuracy of FBS algorithm in both source location and shape recovery, location error (LE), contrast-to-noise ratio (CNR), means square error RMSE, and Dice index were used as the quantitative indexes. The LE measures the distance variation between the centers of the actual region and the reconstructed region. LE is defined as:
where ${L_r}$ is the center of the reconstructed area with the highest value in non-zero value of x. ${L_0}$ is the barycenter of the real fluorescent area. $||\begin {array}{*{20}{c}} . \end {array}|{|_2}$ is the operator of Euclidean distance. A lower LE index indicates that the reconstruction is better. Dice was introduced to evaluate the similarity of the reconstruction area and the real fluorescent area, which defined as: where ${s_r}$ and ${s_0}$ represent the reconstruction area and the actual area, respectively. A high dice index indicates that the two regions have better similarities in both location and morphology.More details, the reconstructed fluorescence area was determined by the non-zero tetrahedron region based on $x$. CNR was performed to demonstrate the contrast of the reconstructed signal and background. Which defined as follows:4. Experiment results
4.1 Numerical simulations results
4.1.1 Feasible experiment
Figure 3 showed the reconstruction results of the feasible experiment. Figure 3(a)-(d) and (e)-(h) show the reconstruction results in 3D and transverse view of the single and double source experiment using $L_{1}$-$L_{2}$ via FBS and three comparing algorithms. In 3D view, the actual and reconstructed fluorescent source were delineated with red meshes and green areas, respectively. Meanwhile, in transverse view the black circles indicated the actual positions of the fluorescent source in the slice of $z=18mm$. For better evaluated the reconstruction results, reconstructed source center and five quantitative indicators are listed in Table 3.
Fig. 3. Result of the feasible experiments (a)-(d) 3D-view and transverse view of the single source reconstructed results (e)-(h) 3D-view and transverse view of the double sources reconstructed results, with the edge-to-edge distance was $4mm$. In 3D view, green areas denote the reconstructed fluorescent sources and red grid represent actual fluorescent source. In transverse view, the black circles represent the actual sources in the slice over the centers of the sources at $z=18mm$.
Table 3. Quantitative results of single and double source reconstruction experiment.
From these results, it is obvious that the recovery source by our proposed algorithm is closer the true source compared with other approaches both in single and double source reconstruction. Furthermore, the quantified results in Table 3 also indicate that proposed $L_{1}$-$L_{2}$ via FBS algorithm achieves the smallest LE and RMSE, and largest Dice similarity and CNR. This indicates that $L_{1}$-$L_{2}$ via FBS algorithm is feasible and more suited for the application of early tumor detection. However, in slice views, $L_{1}$-$L_{2}$ via DCA have some location deviation in x-axis as shown in Fig. 3(b) and Fig. 3(f), which resulted in some artifacts and poor shape recovery performance with smallest Dice similarity and largest location errors. In particular, LE obtained by $L_{1}$-$L_{2}$ via DCA in double source reconstruction is the worst because of the discrete artifacts around the left source. Meanwhile, in double source reconstruction, compare with $L_{1}$-$L_{2}$ via FBS algorithm, two sources given by IVTCG and ISTA-$L_{1/2}$ algorithm have a little tendency to move towards the center with some location deviation from actual area, as shown in Fig. 3(e) and (h). In addition, in single source reconstruction, ISTA-$L_{1/2}$ algorithm has some deviation in z-axis from the 3D view in Fig. 3(d), which is consistent with the largest LE of $0.37mm$ in Table 3.
4.1.2 Robustness experiment
Figure 4 and Fig. 5 shows the robustness of the $L_{1}$-$L_{2}$ via FBS algorithm. In Fig. 4(a)-(d) and (e)-(f) show the reconstructed results in 3D view and transverse view with decreasing the number of excitation source. The corresponding quantitative indicators are presented in Table 4. With the reduction of the number of excitation source, although the reconstruction precision has some different degrees of deterioration, the location error and Dice similarity of all four algorithms can be kept within $0.6mm$ and above $58\%$ in two excitation source case. When one excitation source used, location error and Dice similarity of $L_{1}-L_{2}$ via FBS is still $0.28mm$ and $73\%$ without much deterioration compared with that in 4 excitation source. However, the offset of the green reconstructed area in Fig. 4(f), (g), and (h) are obvious compared with the left corresponding figures. When have one excited source, the location error and Dice similarity of these three contrast algorithm are floating about $0.2mm$ and $6\%$ than that in two excited source. The worst of it is that the Dice similarity of IVTCG is decrease to 35% with a significant shape information losing. In a word, compare with other algorithms, $L_{1}$-$L_{2}$ via FBS algorithm is more robust in terms of location error and Dice similarity when the number of excited sources is decreased.
Fig. 4. Result of the different number excite source experiment (a)-(d) 3D-view and transverse view of two excitation source experiment results (e)-(h) 3D-view and transverse view of the one excitation sources experiment results. In 3D view, green areas denote the reconstructed fluorescent sources and red grid represent actual fluorescent source. In transverse view, the black circles represent the actual sources in the slice over the centers of the sources at $z=18mm$.
Fig. 5. Result of the different edge-to-edge distance experiment (a)-(d) 3D-view and transverse view of the edge-to-edge distance was $2mm$ experiment results (e)-(h) 3D-view and transverse view of the edge-to-edge distance was $1mm$ experiment results. In 3D view, green areas denote the reconstructed fluorescent sources and red grid represent actual fluorescent source. In transverse view, the black circles represent the actual sources in the slice over the centers of the sources at $z=18mm$.
Table 4. Quantitative results of different number excitation source experiment.
Furthermore, the robustness in double source reconstruction were proved by Fig. 5 and Table 5. Figure 5(a)-(d) and (e)-(h) show the results in 3D view and transverse view of the different edge-to-edge distance experiment.The reconstruction challenges increased with decreasing edge-to-edge distance. And these robustness experiment also reflects the spatial resolving capability. From these results, we can see that $L_{1}$-$L_{2}$ via FBS algorithm can distinguish two source with edge-to-edge distance of $1mm$ clearly, as shown in Fig. 5(e). Meanwhile, the average location error and Dice of $L_{1}$-$L_{2}$ via FBS can reach $0.35mm$ and $63\%$ with a good performance in terms of location accuracy, shape recovery capability. In addition, from a series of comparative results, the performance of ISTA-$L_{1/2}$ would be the second-best in terms of LE and Dice. However, the reconstructed image of $L_{1}$-$L_{2}$ via DCA were blurred with some image artifacts even when edge-to-edge distance was $2mm$, as shown in Fig. 5(b). Thus, the corresponding Dice, CNR and LE were worst in Table 5. From Fig. 5(f), two sources given by $L_{1}$-$L_{2}$ via DCA cannot be distinguish because of many existing discrete artifacts. Similarly to Fig. 3(f), the two reconstructed source have a tendency to move towards the center when edge-to-edge distance is $2mm$. These tendencies are even more pronounced when the edge-to-edge distance is $1mm$ as shown in Fig. 5(g). In conclusion, all these results demonstrated that the $L_{1}$-$L_{2}$ via FBS algorithm have a best robustness in dual-sources resolution.
Table 5. Quantitative results of different edge-to-edge distance experiment
4.1.3 In-vivo imaging experiment
In in-vivo imaging experiments, to determine the actual source region, MRI image was adopted to draw the outline of the glioma. The reconstructed results of four comparing algorithm in 3D view, transverse view were displayed in Fig. 6(a)-(d), respectively. The reconstructed images in the CT coordinate system were merged with the MRI slice, as shown in Fig. 6(e)-(h). The registration process can divide into three steps: First, by using bwperim and imfill function in MALTAB, we obtained the contour of the CT and MRI images; Second, we used imshowpair and imregconfig function in MALTAB to obtain the optimization metrics and registration parameters; Third, CT and MRI image was put into imregister function and conducted MRI/CT registration based on the obtained optimization metrics. The white outline is the registration results of the contour of cross-sectional images in the plane $z=17mm$. And the color area and blue dotted line represent the reconstruction tumor region and the actual tumor margin. Meanwhile, to analyze the results quantitatively, reconstructed source center, LE and Dice were calculated and summarized in Table 6.
Fig. 6. Reconstruction results of the in-vivo experiments. (a)-(d) 3D-view and transverse view of in-vivo experiment result. (e)-(h) Images fused by the transverse view and MRI image. In 3D-view the red areas represent reconstruction result. In fused images the white curve is the contour of cross-sectional images, the red area is the reconstruction result, and the blue dotted line outlines the actual location of the tumor.
Table 6. Quantitative results of in-vivo experiment
According to the in-vivo reconstruction results, the tumors were successfully reconstructed using four different approaches. However, for $L_{1}$-$L_{2}$ via DCA and IVTCG algorithm, the reconstructed tumor were tend to brain interior with large location error and small Dice similarity. Fortunately, the reconstructed area obtained by the $L_{1}$-$L_{2}$ via FBS and ISTA-$L_{1/2}$ algorithm were consistent with the actual tumor area in MRI image, and their Dice value can reach 84% and 70% respectively. However, as can be from Table 6, the location deviation of the ISTA-$L_{1/2}$ algorithm mainly manifest in z-axial with LE of $1.24mm$. In general, $L_{1}$-$L_{2}$ via FBS algorithm showed the best accuracy with the least LE and largest Dice similarity. These results further revealed the superior performance of $L_{1}$-$L_{2}$ via FBS in obtaining the morphology and location tracking of the in-vivo fluorescence probe distribution in orthotopic glioma mouse models.
5. Discussion and conclusion
As a promising preclinical imaging technique, FMT has been paid much attention in imaging theories, acquisition equipment, and biomedical applications. However, because of the severe ill-posedness of inverse reconstruction, the accuracy has limited in many biomedical applications. In this work, a novel $L_{1}$-$L_{2}$ norm regularization via forward-backward splitting ($L_{1}$-$L_{2}$ via FBS) method was proposed for better recovering the 3D distributions of early tumor. Theoretically, $L_{1}$-$L_{2}$ via FBS algorithm employ proximity operators and gradient descent to obtain the search point and full consider the highly coherent of the system matrix of FMT. Meanwhile, it handles non-differentiable objectives and non-convex constraints caused by $L_{1}$-$L_{2}$ via FBS norm by splitting the objective to be minimized into simpler functions that are dealt with individually. Thus, $L_{1}$-$L_{2}$ via FBS method outperforms $L_{p}$-norm $(0<p<1)$ and $L_{1}$-norm in terms of promoting sparsity when sensitivity matrix A is highly coherent. Furthermore, FBS iterations are convergence and not trapped at stationary points, which can guarantee the robust of the solution.
To validate the performance of the proposed $L_{1}$-$L_{2}$ via FBS algorithm, numerical simulations and in-vivo experiments were conducted. $L_{1}$-$L_{2}$ via DCA, IVTCG and ISTA-$L_{1/2}$ were employed for qualitative and quantitative comparisons. In numerical simulations, feasibility and robustness experiments were conducted. The feasibility experiment results show that $L_{1}$-$L_{2}$ via FBS algorithm has a best performance in terms of location error and Dice similarity both in single and double source experiment. Which demonstrated the advantage of $L_{1}$-$L_{2}$ via FBS algorithm in location accuracy, and morphology recovery for early tumor detection. In robustness experiments, with the gradually reducing of the excitation source, the ill-posedness of FMT reconstruction is aggravated. Fortunately, the performance of the $L_{1}$-$L_{2}$ via FBS algorithm was still satisfactory. Even in the extreme case where only one excited node was used, the proposed algorithm still successfully recovered the source with a deviation of $0.21mm$. Which demonstrates that the accuracy of location and morphology of the proposed algorithm remained steady. And it also indicates the potential of the $L_{1}$-$L_{2}$ via FBS algorithm in single view FMT reconstruction. By considering the application requirements for multiple tumor detection, with the edge-to-edge distance decreasing, the spatial resolving capability of inverse reconstruction algorithms was verified in another robustness experiments. From these results, we can see that the $L_{1}$-$L_{2}$ via FBS algorithm can distinguish two sources with an edge-to-edge distance of 1 mm clearly. While, the tumor region reconstructed by the comparison algorithms exhibited problems of inaccuracy, spatial discontinuity and image blur with some image artifacts. Lastly, in in-vivo experiment, tumor-bearing mouse model was developed to further evaluate the practicability of $L_{1}$-$L_{2}$ via FBS algorithm. The in-vivo experimental results exhibited that $L_{1}$-$L_{2}$ via FBS algorithm was great practicality for early tumor detection of living animals. And compared to the other methods, the region reconstructed using $L_{1}$-$L_{2}$ via FBS algorithm was a better approximation of the actual tumor distribution. It further revealed the superior performance of $L_{1}$-$L_{2}$ via FBS algorithm in obtaining the morphology and location tracking of the in-vivo fluorescence probe distribution. In summary, this result demonstrates the usefulness of the $L_{1}$-$L_{2}$ via FBS algorithm in resolving multiple sources and improving reconstruction quality and spatial resolution.
It needs to be emphasized that in this work, the absorption and scattering coefficients assigned to each organ were estimated on the basis of a compilation of relevant tissue optical property data reported in the literature [40]. In the pre-clinical practical application, the educated guesses of the optical properties will influence the accuracy of the reconstructed results. In order to solve this problem, some more effective measuring method, like diffusion optical tomography (DOT) [41], or some other near infrared imaging method will be adopted to measure the optical properties of tissue and background. Based on that, in in-vivo experiment, the reconstruction accuracy of our proposed algorithm will be further improved. Besides that, due to the lack of absolute calibration technique for photon number and the standard for determining the distribution of tumors in living animals, in in-vivo experiment, we did not quantify the CNR and RMSE. In the future, based on the integrating sphere method [42], it is necessary to quantify the reconstruction results precisely by adopting some quantitative analysis method and multimodal imaging technology.
Although $L_{1}$-$L_{2}$ via FBS algorithm has achieved much better results, the work on this breed of algorithms is far from done, as there are various analysis questions that remain open. Firstly, in our research, $L_{1}$-$L_{2}$ was limited by the equal penalty for the $L_{1}$ and $L_{2}$ norms. Considering the versatility, it can be generalized by considering the ${L_1}$-$\alpha {L_2}$ metric for $\alpha \ge 0$ in further work. Secondly, an analytical solution for the proximal operator of the $L_{1}$-$L_{2}$ metric was derived, leading to classical fast $L_{1}$ norm solvers such as fast iterative soft-thresholding algorithm and alternating direction method of multipliers (ADMM) that is applicable for $L_{1}$-$L_{2}$ norm. Therefore, some comparison study will be conducted to confirm the efficient of these optimal algorithms. Lastly, our research mainly focus on the application of FMT in early tumor detection. For some other preclinical application, the types of probe distributions will have a significant difference, which may reduce the performance of our proposed sparse regularization algorithm in some degree. In this case, $L_{1}$-$L_{2}$ norm should be coupled with some other non-sparse regularization terms to recovery the continuous bigger area probe distributions.
In conclusion, a novel $L_{1}$-$L_{2}$ norm regularization via forward-backward splitting method was proposed to improve the reconstruction accuracy of FMT in this paper. The $L_{1}$-$L_{2}$ norm regularization was designed by full considering the highly coherent of the system matrix of FMT and maintaining the simplicity of gradient descent algorithms and the sparsity of results. Compared with several sparse reconstruction algorithms, $L_{1}$-$L_{2}$ via FBS algorithm performs better in terms of location accuracy, dual-source resolution capability, robustness and in vivo practicability. We believed that this study would promote the preclinical and clinical applications of FMT in early tumor detection.
Funding
National Outstanding Youth Science Fund Project of National Natural Science Foundation of China (11871321, 61901374, 61906154, 61971350); Natural Science Foundation of Shaanxi Province (2019JQ-724); Postdoctoral Innovative Talents Support Program (BX20180254); Scientific and Technological projects of Xi’an (201805060ZD11CG44); Key Research and Development Program of Shaanxi (2020SF-036); Xi’an Science and Technology Project (2019218214GXRC018CG019-GXYD18.3).
Acknowledgment
The authors would like to thank the CAS key Laboratory of Molecular Imaging for in-vivo experiment data.
Disclosures
The authors declare no potential conflict of interests.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request
References
1. V. Ntziachristos, C. Bremer, and R. Weissleder, “Fluorescence imaging with near-infrared light: new technological advances that enable in vivo molecular imaging,” European Radiology 13(1), 195–208 (2003). [CrossRef]
2. R. Weissleder, “Scaling down imaging: molecular mapping of cancer in mice,” Nat. Rev. Cancer 2(1), 11–18 (2002). [CrossRef]
3. 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]
4. M. Koch, P. Symvoulidis, and V. Ntziachristos, “Tackling standardization in fluorescence molecular imaging,” Nat. Photonics 12(9), 505–515 (2018). [CrossRef]
5. 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]
6. L. Kong, Y. An, Q. Liang, L. Yin, Y. Du, and J. Tian, “Reconstruction for fluorescence molecular tomography via adaptive group orthogonal matching pursuit,” IEEE Trans. Biomed. Eng. 67(9), 2518–2529 (2020). [CrossRef]
7. J. K. Willmann, N. Van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discovery 7(7), 591–607 (2008). [CrossRef]
8. A. Ale, V. Ermolayev, E. Herzog, C. Cohrs, M. H. De Angelis, and V. Ntziachristos, “Fmt-xct: in vivo animal studies with hybrid fluorescence molecular tomography–x-ray computed tomography,” Nat. Methods 9(6), 615–620 (2012). [CrossRef]
9. X. He, H. Meng, X. He, K. Wang, X. Song, and J. Tian, “Nonconvex laplacian manifold joint method for morphological reconstruction of fluorescence molecular tomography,” Mol Imaging Biol 23(3), 394–406 (2021). [CrossRef]
10. H. Meng, K. Wang, Y. Gao, Y. Jin, X. Ma, and J. Tian, “Adaptive gaussian weighted laplace prior regularization enables accurate morphological reconstruction in fluorescence molecular tomography,” IEEE Trans. Med. Imaging 38(12), 2726–2734 (2019). [CrossRef]
11. A. Jin, B. Yazici, A. Ale, and V. Ntziachristos, “Preconditioning of the fluorescence diffuse optical tomography sensing matrix based on compressive sensing,” Opt. Lett. 37(20), 4326–4328 (2012). [CrossRef]
12. H. Meng, Y. Gao, X. Yang, K. Wang, and J. Tian, “K-nearest neighbor based locally connected network for fast morphological reconstruction in fluorescence molecular tomography,” IEEE Trans. Med. Imaging 39(10), 3019–3028 (2020). [CrossRef]
13. X. He, X. Wang, H. Yi, Y. Chen, X. Zhang, J. Yu, and X. He, “Laplacian manifold regularization method for fluorescence molecular tomography,” J. Biomed. Opt. 22(4), 045009 (2017). [CrossRef]
14. H. Wang, C. Bian, L. Kong, Y. An, Y. Du, and J. Tian, “A novel adaptive parameter search elastic net method for fluorescent molecular tomography,” IEEE Trans. Med. Imaging 40(5), 1484–1498 (2021). [CrossRef]
15. S. Zhang, X. Ma, Y. Wang, M. Wu, H. Meng, W. Chai, X. Wang, S. Wei, and J. Tian, “Robust reconstruction of fluorescence molecular tomography based on sparsity adaptive correntropy matching pursuit method for stem cell distribution,” IEEE Trans. Med. Imaging 37(10), 2176–2184 (2018). [CrossRef]
16. H. Guo, X. He, M. Liu, Z. Zhang, Z. Hu, and J. Tian, “Weight multispectral reconstruction strategy for enhanced reconstruction accuracy and stability with cerenkov luminescence tomography,” IEEE Trans. Med. Imaging 36(6), 1337–1346 (2017). [CrossRef]
17. H. Guo, J. Yu, Z. Hu, H. Yi, Y. Hou, and X. He, “A hybrid clustering algorithm for multiple-source resolving in bioluminescence tomography,” J. Biophotonics 11(4), e201700056 (2018). [CrossRef]
18. Y. Zhou, M. Chen, H. Su, and J. Luo, “Self-prior strategy for organ reconstruction in fluorescence molecular tomography,” Biomed. Opt. Express 8(10), 4671–4686 (2017). [CrossRef]
19. N. Cao, A. Nehorai, and M. Jacob, “Image reconstruction for diffuse optical tomography using sparsity regularization and expectation-maximization algorithm,” Opt. Express 15(21), 13695–13708 (2007). [CrossRef]
20. S. Jiang, J. Liu, G. Zhang, Y. An, H. Meng, Y. Gao, K. Wang, and J. Tian, “Reconstruction of fluorescence molecular tomography via a fused lasso method based on group sparsity prior,” IEEE Trans. Biomed. Eng. 66(5), 1361–1371 (2019). [CrossRef]
21. 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]
22. 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]
23. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]
24. Y. Liu, J. Yu, X. Qin, and J. Guo, “The split bregman iteration algorithm for bioluminescence tomography,” Sci. Sin.-Inf. 44(2), 284–294 (2014). [CrossRef]
25. C. Leng, D. Yu, S. Zhang, Y. An, and Y. Hu, “Reconstruction method for optical tomography based on the linearized bregman iteration with sparse regularization,” Computational and mathematical methods in medicine 2015, 1–8 (2015). [CrossRef]
26. X. He, J. Yu, X. Wang, H. Yi, Y. Chen, X. Song, and X. He, “Half thresholding pursuit algorithm for fluorescence molecular tomography,” IEEE Trans. Biomed. Eng. 66(5), 1468–1476 (2019). [CrossRef]
27. H. Guo, Z. Hu, X. He, X. Zhang, M. Liu, Z. Zhang, X. Shi, S. Zheng, and J. Tian, “Non-convex sparse regularization approach framework for high multiple-source resolution in cerenkov luminescence tomography,” Opt. Express 25(23), 28068–28085 (2017). [CrossRef]
28. Y. Lou, P. Yin, Q. He, and J. Xin, “Computing sparse representation in a highly coherent dictionary based on difference of l_1 and l_2,” J Sci Comput 64(1), 178–196 (2015). [CrossRef]
29. H. Zhang, G. Geng, X. Wang, X. Qu, Y. Hou, and X. He, “Fast and robust reconstruction for fluorescence molecular tomography via regularization,” BioMed Res. Int. 2016, 1–9 (2016). [CrossRef]
30. Y. Lou and M. Yan, “Fast l1–l2 minimization via a proximal operator,” J Sci Comput 74(2), 767–785 (2018). [CrossRef]
31. X. Zhang, W. Zhou, X. Zhang, and H. Gao, “Forward–backward splitting method for quantitative photoacoustic tomography,” Inverse problems 30(12), 125012 (2014). [CrossRef]
32. Q. Ding, G. Chen, X. Zhang, Q. Huang, H. Ji, and H. Gao, “Low-dose ct with deep learning regularization via proximal forward–backward splitting,” Phys. Med. Biol. 65(12), 125009 (2020). [CrossRef]
33. X. He, J. Liang, X. Wang, J. Yu, X. Qu, X. Wang, Y. Hou, D. Chen, F. Liu, and J. Tian, “Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method,” Opt. Express 18(24), 24825–24841 (2010). [CrossRef]
34. A. X. Cong and G. Wang, “A finite-element-based reconstruction method for 3d fluorescence tomography,” Opt. Express 13(24), 9847–9857 (2005). [CrossRef]
35. M. Schweiger, S. Arridge, M. Hiraoka, and D. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11), 1779–1792 (1995). [CrossRef]
36. D. Wang, X. Liu, Y. Chen, and J. Bai, “A novel finite-element-based algorithm for fluorescence molecular tomography of heterogeneous media,” IEEE Trans. Inform. Technol. Biomed. 13(5), 766–773 (2009). [CrossRef]
37. 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(18), 8211–8223 (2006). [CrossRef]
38. T. K. Pong, P. Tseng, S. Ji, and J. Ye, “Trace norm regularization: reformulations, algorithms, and multi-task learning,” SIAM J. Optim. 20(6), 3465–3489 (2010). [CrossRef]
39. H. Liao and M. K. Ng, “Blind deconvolution using generalized cross-validation approach to regularization parameter estimation,” IEEE Trans. on Image Process. 20(3), 670–680 (2011). [CrossRef]
40. 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]
41. W. Lu, J. Duan, D. Orive-Miguel, L. Herve, and I. B. Styles, “Graph-and finite element-based total variation models for the inverse problem in diffuse optical tomography,” Biomed. Opt. Express 10(6), 2684–2707 (2019). [CrossRef]
42. O. A. Savchuk, J. Carvajal, J. Massons, M. Aguiló, and F. Díaz, “Determination of photothermal conversion efficiency of graphene and graphene oxide through an integrating sphere method,” Carbon 103, 134–141 (2016). [CrossRef]
