Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Adaptive shrinking reconstruction framework for cone-beam X-ray luminescence computed tomography

Open Access Open Access

Abstract

Cone-beam X-ray luminescence computed tomography (CB-XLCT) emerged as a novel hybrid technique for early detection of small tumors in vivo. However, severe ill-posedness is still a challenge for CB-XLCT imaging. In this study, an adaptive shrinking reconstruction framework without a prior information is proposed for CB-XLCT. In reconstruction processing, the mesh nodes are automatically selected with higher probability to contribute to the distribution of target for imaging. Specially, an adaptive shrinking function is designed to automatically control the permissible source region at a multi-scale rate. Both 3D digital mouse and in vivo experiments were carried out to test the performance of our method. The results indicate that the proposed framework can dramatically improve the imaging quality of CB-XLCT.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

As a novel hybrid X-ray CT/molecular imaging modality, X-ray luminescence computed tomography (XLCT) can be used for drug research, metabolic research and clinical diagnosis [13]. Compared with conventional optical molecular modalities, i.e. fluorescence molecular tomography (FMT) [4] and bioluminescence tomography (BLT) [5], XLCT has been shown to be able to avoid obvious background noise while imaging deep tissue. Furthermore, XLCT has been used in small animals in vivo [68].

In general, narrow-beam XLCT, pencil-beam XLCT (PB-XLCT) [9], and cone-beam XLCT (CB-XLCT) [13] are three main types of common XLCT. Xing et al., firstly proposed a narrow-beam XLCT to realize the 3D reconstruction of nanoparticles [2]. Further, using a collimated pencil-beam X-ray, Li et al., designed a pencil-beam XLCT to recover the deep targets [9]. Besides, narrow beam XLCT was proposed based on a limited-view imaging technique [10]. However, although the above two XLCT techniques have high spatial resolution, long data acquisition time significantly limits their application in drug research and early tumor detection [11]. To resolve the problem, Chen et al., proposed a cone-beam XLCT (CB-XLCT) imaging system, which can sharply reduce the imaging time and improve the X-ray dose utilization efficiency [12].

Directly caused by the insufficiency of external measurements, high ill-posedness is still a technical bottleneck for XLCT [3]. Extensive researches have been conducted to improve reconstruction results. Chen et al., used the multi-spectral data to realize quantitative reconstruction [13]. Liu et al., proposed a single-view reconstruction by a wavelet transform method [14]. Gao et al., presented a truncated singular value decomposition (TSVD) approach for CB-XLCT imaging [7]. Generally, a prior information, such as structural prior, optical properties, sparsity of target distribution and permissible source region (PSR) can effectively improve reconstruction results [1519].

The PSR method can reduce the scale of inverse problem by reducing the number of unknown variables, and has been widely employed to BLT, FMT and XLCT [13,2024]. However, the permissible region (PR) is traditionally determined by a rough and subjective estimation from the surface photon distribution [13]. When the light source is deep in imaging body, the accuracy of the method is affected [25]. Although the iterative estimation can be adopted to improve the above deficiency, the threshold for the selection of PR is usually a fixed manually set rate, leading to excessive iterations and high computational cost [26]. To overcome this obstacle, we developed an adaptive shrinking permissible source region (ADS_PSR) framework for the imaging of small targets without a prior information in this study. Meanwhile, a multi-scale kernel function is designed to accelerate the shrinking rate automatically. To the best of our knowledge, this is the first time that a special PSR framework was proposed to alleviate the ill-posedness of inverse problem for CB-XLCT imaging. Several numerical simulation experiments and an in vivo experiment have been performed to validate the effectiveness and robustness of the proposed framework.

The contributions of the paper can be summarized as follows: 1) a novel PSR framework was proposed to alleviate the ill-posedness of inverse problem of CB-XLCT; 2) instead of the traditionally manual experience, a new shrinking function was designed to automatically obtain permissible region; and 3) the proposed framework is suitable for combining with mainstream reconstruction algorithms for CB-XLCT.

In Section 2, the imaging model of CB-XLCT and the proposed ADS_PSR framework are introduced. Then, numerical simulation experiments and the in vivo experiment are demonstrated in Section 3. Finally, we discuss the results and draw a conclusion in Section 4.

2. Method

2.1 Imaging system of the CB-XLCT

As illustrated in Fig. 1, the traditional scheme diagram of CB-XLCT imaging system mainly includes the following parts: a cone beam X-ray source to achieve X-ray excitation, a CMOS X-ray detector panel to collect X-ray projection data, an electron-multiplying CCD (EMCCD) camera to collect luminescent data for optical imaging, and a rotation stage to realize multi-angles X-ray excitation. Hence, X-ray CT imaging and 3D optical imaging can be achieved together by CB-XLCT system.

 figure: Fig. 1.

Fig. 1. Schematic diagram of CB-XLCT system [1112].

Download Full Size | PDF

2.2 Imaging model of the XLCT

In detail, X-rays are emitted by the cone-beam X-ray source and traveled through biological tissues. Based on Beer-Lambert’s law, the propagation process of X-ray traveled in tissues can be modeled as follows [11]:

$$X(r) = {X_0}exp\{ - \int_{{r_0}}^r {{\mu _x}(\tau )d\tau } \}$$
where $X(r)$ is the X-ray intensity at position r in tissue, ${X_0}$ is the initial X-ray intensity at ${r_0}$, and ${\mu _x}$ is the X-ray attenuation coefficient which can be calculated via an attenuation-based CT technique [1].

Sequentially, once nanoparticles distributed in the imaging object are irradiated by X-rays, the visible or NIR light will be emitted. This emission light can be expressed as the following linear relationship:

$$\textrm{S(}r) = \varepsilon X(r)\rho (r)$$
where $\textrm{S(}r)$ is the target, which is often used to mimic the small tumor. $\rho (r)$ is the nanophosphor concentration at position r, and $\varepsilon$ is the light yield of nanoparticles.

Due to the highly scattering and weakly absorbing properties of light in biological tissues, the visible or near infrared (NIR) light transmitted in the tissues can be modeled by the following diffusion equation (DE) model with the Robin boundary condition as follows [2728]:

$$\left\{ \begin{array}{ll} - \nabla [D(r)\nabla \Phi (r)] + {\mu_a}(r)\Phi (r) = S(r) &r \in \Omega \\ \Phi (r) + 2\kappa D(r)[\nu \nabla \Phi (r)] = 0 &r \in \partial \Omega \end{array} \right.$$
where $D(r) = {(3({\mu _a}(r) + \mu_s^{\prime}(r)))^{ - 1}}$ is the diffusion coefficient, where ${\mu _a}(r)$ and $\mu_s^{\prime}(r)$ is the absorption and reduced scattering coefficient, respectively. $\Omega$ is the domain of the image object and $\partial \Omega$ is the corresponding boundary. $\Phi (r)$ is the photon fluence at position r. $\nu$ is the outward unit normal vector to boundary, and $\kappa$ is a constant describing the optical reflective index mismatch.

Based on the finite-element (FEM) method [29], Eqs. (1) and (3) can be converted into a linear relationship between the unknown nanophosphor distribution $\rho$ and the NIR measurement J as following matrix-form equation:

$$J = A\rho$$
where A is the system matrix which is used to map the unknown $\rho$ to known measurement J [8].

2.3 Reconstruction problem

The goal of CB-XLCT imaging is to recover the unknown distribution of the nanophosphor based on the NIR measurement J captured by the EMCCD camera. However, the high scattering character of light in biological tissue leads to a poorly conditioned system matrix $A$. Thus, that it is impractically to solve $\rho$ directly from Eq. (4). In addition, the usually limited optical measurements and other experiment conditions further increase the difficulty of solution to the above inverse problem. The distribution of nanophosphors in biological tissue are usually small and sparse, compared to the entire imaging body [1214]. Herein, the compressed sensing (CS) theory can be applied to solve $\rho$ by converting Eq. (4) into the following optimization problem [30].

$$\min ||{A\rho - J} ||_2^2 + \tau {||\rho ||_p}$$
where $\tau$ is the regularization coefficient. ${||\rho ||_p}$ is the ${L_p}$-norm of $\rho$. When $p = 0$, Eq. (5) is ${L_0}$ norm, which is NP-hard. When $p = 1$, Eq. (5) is the well-known ${L_1}$-norm as the convex relaxation of ${L_0}$-norm, which has been a widely used sparsity-inducing norm for CB-XLCT imaging [13]. When $p = 2$, Eq. (5) is the ${L_2}$-norm, which is the commonly used Tikhonov regularization [31].

2.4 Adaptive shrinking permissible source region framework for CB-XLCT imaging

Generally, the size of system matrix A in Eq. (5) is usually large and this obviously aggravates the ill-posedness of the inverse problem [11,14]. To achieve high imaging quality, ADS_PSR framework has been developed in this study. The implementation of the ADS_PSR framework is started from the entire domain to a small area without a prior information. From every shrinking iteration $k$, the permissible region (PR) ${R_k}$ is updated by removing nodes with lower recovered density. We define a following vector ${P_k}$ to reduce the scale of inverse problem:

$${P_k} = \left\{ \begin{array}{l} 1\;\; \hbox{if node}\; i \;\hbox{is within the permissible region} \\ 0\;\; \textrm{otherwise} \end{array} \right.$$
As soon as ${P_k}$ is determined, A and $\rho$ in Eq. (5) become ${A_k} = {A_{k - 1}} \otimes {P_k}$ and ${\rho _k} = {\rho _{k - 1}} \otimes {P_k}$, respectively. ${\otimes}$ represents the operation that removes columns in previous matrix ${A_{k - 1}}$ with zero-element in column vector ${P_k}$. Then, an updated inverse problem has been converted as follows:
$$\mathop {\min }\limits_{{R_k}} ||{{A_k}{\rho_k} - J} ||_2^2 + {\tau _k}{||{{\rho_k}} ||_p}$$
After solving Eq. (7), a new PR is determined by ranking the magnitude of ${X_k}$ in the decreasing order for next iterative reconstruction.

Aiming to automatically select the group of mesh nodes to obtain a shrinking ${\rho _k}$ without manual intervention, we designed an adaptive kernel function to automatically speed up the shrinking process at a multi-scale rate as follows:

$${\zeta _k} = 1 - \frac{{\beta \cdot \exp [{{ - (k - 1)} \mathord{\left/ {\vphantom {{ - (k - 1)} \alpha }} \right.} \alpha }]}}{{1 + \beta \cdot \exp [{{ - (k - 1)} \mathord{\left/ {\vphantom {{ - (k - 1)} \alpha }} \right.} \alpha }]}}$$
where k is the iterations which is set from 1. $\alpha$ and $\beta$ is multi-scale factor and nonnegative coefficient, respectively. $\alpha$ and $\beta$ are all determined experientially to confine the shrinking rate. For the reconstruction, the range of $\alpha$ and $\beta$ is $[2.5,3.5]$ and $[3.5,7.5]$, respectively. Figure 1 shows the value of $\zeta$ function with the increase of k ($\alpha = 3,\; \; \;\beta = 6$).

In the implementation, $length(R)$ is the number of nodes related to the PR determined in the current iteration. Next, the number of nodes for the next round reconstruction with higher recovered density is automatically obtained by $\lfloor{Num(R)\cdot {\zeta_k}} \rfloor$ and $\lfloor{} \rfloor$ is the integer operation. With the increase of $K$, the trend of solution $\rho ({R_i})$ gradually to be stable, which means that the difference of number of nodes between the adjacent rounds is gradually reduced. Especially, as can be seen in Fig. 2, about 90% of the nodes of the previous solution ${\rho _{k - 1}}$ will be maintained as the value of k increases to 10, compared with 10% of that at the beginning. Meanwhile, the number of discarded nodes with lower recovered density decrease correspondingly. Consequently, the PR will correspondingly shrink to a stable area and a stable solution $\rho$ can be acquired. Thus, the function ${\zeta _k}$ can adaptively fit this case with a multi-scale rate. Table 1 shows details of the proposed framework for CB-XLCT imaging.

 figure: Fig. 2.

Fig. 2. The multi-scale of $\zeta$ function with the increase of $k$.

Download Full Size | PDF

Tables Icon

Table 1. Procedure of the adaptive shrinking permissible source region framework for CB-XLCT imaging

3. Experiments and results

In this section, several numerical simulation experiments and an in vivo experiment were investigated to evaluate the performance of the proposed framework. For detailed comparison, the traditional iterative shrinking permissible region (ISPR) method was used in our study [1923]. Besides, some widely-used algorithms in CB-XLCT, including the algebraic reconstruction technique (ART) [32], the incomplete variables truncated conjugate gradient method for sparse ${L_1}$-norm minimization (IVTCG) [33] and the Tikhonov regularization algorithm for ${L_2}$-norm minimization [34] were adopted. Considering that ${L_1}$-norm has been proved efficient for the sparsity-inducing norm for detecting the sparsely distributed tumors [13], only our proposed framework and ISPR method incorporate IVTCG algorithm. The threshold value of ISPR was set as 50% and 70% according to previous works [2223]. In addition, the regularization coefficient $\tau$ of IVTCG and Tikhonov is automatically selected by the L-curve method [3536]. The multi-scale factor $\alpha$ and nonnegative coefficient $\beta$ of ${\zeta _k}$ function were set empirically as 3 and 5, respectively.

Several quantitative indexes, such as absolute location error (LE) [11], recovered density (RD) [13], relative quantity error (RQE) [13] and the percentage of non-zero coefficient (PNZ) [37] were adopted for assessment. The experiment codes were written in MATLAB and was performed on a desktop computer with 2.20 GHz Intel Xeon Processor I7-4702MQ and 16G RAM.

3.1 Numerical simulation experiments

A mouse atlas of CT data, as demonstrated in Fig. 3(a), with a single nanophosphor target and double nanophosphor targets respectively was adopted to test our method. The single nanophosphor target was located with center of (17.8, 6.6, 50.8) mm. The double nanophosphor targets were located with center of (11.2, 14.6, 60) mm and (22.8, 14.6, 60) mm. The 2D views of CT slices of the digital mouse model with single-target and double-targets were shown in Fig. 3(d) and Fig. 3(e), respectively. All targets were designed as a small ball with a radius of 1 mm to mimic small tumor. The voltage and current of X-ray source were set as 50kVp voltage and 1 mA, respectively. The attenuation coefficient of X-ray was set as 0.0535 mm−1 [38]. The homogeneous absorption coefficient and reduced scattering coefficient of the two digital mouse models were set as 0.3 mm−1 and 10 mm−1, respectively [39]. The concentration of all targets was 0.3183 $\mu g/m{m^3}$, and the corresponding light yield were 0.15 cm3/mg [12].

 figure: Fig. 3.

Fig. 3. Reconstruction model in simulation experiment. (a) 3D Mouse Model (b) Forward simulation of single-target case (c) Forward simulation of double-targets case (d) 2D view of CT slice for the single-target case (e) 2D view of CT slice for the double-targets case.

Download Full Size | PDF

The digital mouse model for single-target case was discretized into 47733 nodes and 251509 tetrahedral elements. The digital mouse model for double-targets case was discretized into 47656 nodes and 251161 tetrahedral elements. A cone beam X-ray source was used to excite the two digital mouse models with every 20 deg intervals during a 360-deg scan. By using the Monte Carlo (MC) method, which is implemented by the Chinese Academy of Sciences Molecular Optical Simulation Environment (MOSE) software [40], we could get the forward simulation results for the single-target case and the double-targets case as demonstrated in Fig. 3(b) and Fig. 3(c), respectively.

3.1.1 Single-target experiment

For the single-target case, a mesh with 6367 nodes and 32158 tetrahedral elements (average mesh size of about 0.41 mm) was used for the inverse reconstruction. The 2D views of the recovered results at $Z = 50.8$ mm plane by ADS_PSR + IVTCG, IPSR with artificial threshold 50% + IVTCG, IPSR with artificial threshold 70% + IVTCG, IVTCG with no processing, Tikhonhov with no processing and ART with no processing are shown in Fig. 4, where the real position of nanophosphor target is marked as the red circle. Table 2 presents the corresponding quantitative results.

 figure: Fig. 4.

Fig. 4. The 2D views ($Z = 50.8$ mm plane) of the reconstructed results for the single-target reconstruction.

Download Full Size | PDF

Tables Icon

Table 2. Quantitative results of the single-target case

From Fig. 4 and Table 2, it is significant that ADS_PSR and ISPR outperform the traditional algorithms. The ADS_PSR method with IVTCG achieves the least LE and RQE, the highest RD. However, with no processing, the LE of IVTCG is larger the 1 mm, which means the improvement of reconstruction by the proposed framework is remarkable. Furthermore, the PNZ of the ADS_PSR method with IVTCG is the lowest, denoting that the solution of our proposed method is the sparest. As shown in Fig. 5, the iteration number of ADS_PSR is significantly less than that of ISPR with artificial threshold.

 figure: Fig. 5.

Fig. 5. The iteration number of ADS_PSR and ISPR with artificial threshold in single-target case.

Download Full Size | PDF

3.1.2 Double-targets case

In the double-targets simulation experiment, a mesh with 6746 nodes and 35189 tetrahedral elements (average mesh size of about 0. 38 mm) was adopted. Figure 6 gives the 2D views of the recovered results at $Z = 60$ mm plane by ADS_PSR + IVTCG, IPSR with artificial threshold 50% + IVTCG, IPSR with artificial threshold 70% + IVTCG, IVTCG with no processing, Tikhonhov with no processing and ART with no processing. Table 3 presents the corresponding results and Fig. 7 is the iteration number of ADS_PSR and ISPR with artificial threshold in double-targets case.

 figure: Fig. 6.

Fig. 6. 2D views ($Z = 60$ mm plane) of the reconstructed results for the double-targets reconstruction.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. The iteration number of ADS_PSR and ISPR with artificial threshold in double-targets case.

Download Full Size | PDF

Tables Icon

Table 3. Quantitative results of the double-targets case

As can be seen in Figs. 67, and Table 3, all methods can separate the two targets. The proposed ADS_PSR with IVTCG yields the best results. However, the LE of one target are larger than 1 mm by ISPR method and IVTCG algorithm. Besides, the results of Tikhonhov with no processing and ART with no processing are similar to that of single target reconstruction.

3.1.3 Stability analysis with single-target reconstruction

To further assess the stability and robustness of the proposed method for CB-XLCT imaging, we conduct a series of experiments to investigate the reconstructed results in this section. Specifically, the influence of nodes and tetrahedral elements, the influence of measurements, and the influence of the measurement with different noise level were all evaluated with single-target reconstruction.

In the previous experiments, the surface data of the phantom were obtained with 20° intervals for reconstruction. The reduction in measurement data due to decreased view angles will directly aggravate the ill-posed inverse problem and thus affect the reconstruction. Thus, the influence of different view numbers on the proposed method was investigated firstly. We gradually reduced the number of view numbers to 12, 9, 6, and 3. The corresponding quantitative results were reported in Table 4, and Fig. 8 gives the graph results on LE and RD. As shown Fig. 8, although the quality of the proposed is influenced with the reducing of view angles, the LE are all below 0.7 mm. Meanwhile, the value of RD decreased correspondingly. But the value of RD is almost above 0.2 $\mu g/m{m^3}$. Herein, it is obvious that the proposed method is robust in different view angles.

 figure: Fig. 8.

Fig. 8. Illustration of LE and RD at different number of view angles

Download Full Size | PDF

Tables Icon

Table 4. Comparison results for reconstruction single-target by ADS_PSR + IVTCG with different measurements

It is noted that noise was unavoidable in CB-XLCT imaging [13], a group of experiments with noise were also performed. Concretely, Gaussian noise with five different noise levels (5%, 10%, 15%, 20% and 25%) was added to the measurement data of single-target reconstruction. Table 5 and Fig. 9 show the corresponding reconstruction results. Obviously, the LE are all around (even below) 0.5 mm. The RD are all above 0.2 $\mu g/m{m^3}$. It is clear that the proposed method is robust to measurement noise.

 figure: Fig. 9.

Fig. 9. Illustration of LE and RD at different noise levels.

Download Full Size | PDF

Tables Icon

Table 5. Comparison results for the reconstruction single-target by ADS_PSR + IVTCG with different noise levels

To further evaluate the stability of the proposed method, the influence of nodes and tetrahedral elements are investigated. Commonly used different level numbers of nodes and tetrahedral elements (from around 10000 nodes to 2000 nodes) were adopted for reconstruction. The detailed quantitative results are shown in Table 6 and Fig. 10. We found that although the numbers of nodes decreases from 10635 to 2593, the LE changes basically between 0.4 and 0.8. Specially, the number of nodes from 5000 to 6000 give the best results on LE. Similar to the results of previous experiment with noise, the RD with different numbers of nodes and tetrahedral elements are all above 0.2 $\mu g/m{m^3}$.

 figure: Fig. 10.

Fig. 10. Illustration of LE and RD at different number of nodes.

Download Full Size | PDF

Tables Icon

Table 6. Comparison results for the reconstruction single-target by ADS_PSR + IVTCG with different numbers of nodes and tetrahedral elements

3.2 In vivo experiment

To ensure the feasibility of our proposed method in in vivo applications of CB-XLCT, a female eight-week-old mouse was used for the in vivo experimental [14]. The homogeneous absorption coefficient and reduced scattering coefficient of in vivo data were set to be 0.3 mm−1 and 10 mm−1, respectively [39]. The voltage and current of a cone-beam X-ray source (Apogee, Oxford Instruments, USA) with a micro 55-mm f/2.8 lens (Nikkor, Nikon, Japan) were set to be 50 KVp and 1 mA, respectively. A highly sensitive CCD camera (PIXIS 2048, Princeton Instruments, USA) with a micro 55-mm f/2.8 lens (Nikkor, Nikon, Japan) was mounted vertical to the X-ray axis to capture the optical data (620 nm). The integrating time and binning of the CCD were set to be 30 s and $2 \times 2$, respectively. For the optical imaging, the 360 deg luminescent photons data acquired by the CCD with 45 deg intervals was adopted for the inverse reconstruction. Meanwhile, a flat-panel detector (C7921CA-02, Hamamatsu, Japan) was used to fulfill the 360 deg micro-CT imaging with 1 deg intervals. The micro-CT imaging was performed via the Feldkamp-Davis-Kress (FDK) method [41]. The micro-CT result of the female mouse is demonstrated in Fig. 11.

 figure: Fig. 11.

Fig. 11. The micro-CT result of the female mouse.

Download Full Size | PDF

According to the result of stability analysis, the mesh with 5987 nodes and 28816 tetrahedral elements (average mesh size of about 0.46 mm) was utilized in in vivo experiment for inverse reconstruction. The reconstruction results overlaid with CT data are shown in Fig. 12. The detailed information of results are listed in Table 7. Figure 13 is the iteration number of ADS_PSR and ISPR with artificial threshold. As we can see from Fig. 12 and Table 7, the quantitative result of the proposed method is satisfactory for in vivo experiment, where the LE is less than 1 mm. Among all comparison methods, the proposed method yields the sparest solution. Moreover, the solution of L2-norm is not sparse, and the ART fails to obtain feasible reconstruction results.

 figure: Fig. 12.

Fig. 12. Reconstruction results overlaid with CT data of in vivo experiment. (a) coronal-view results at $X = 12$ mm. (b) sagittal-view results at $Y = 13.2$ mm; (c) transversal-view results at $Z = 6.8$ mm.

Download Full Size | PDF

 figure: Fig. 13.

Fig. 13. The iteration number of ADS_PSR and ISPR with artificial threshold of in vivo experiment.

Download Full Size | PDF

Tables Icon

Table 7. Quantitative results of the in vivo experiment

4. Discussion and conclusions

In this study, we proposed an adaptive shrinking permissible source region framework, named ADS_PSR for the recovered result of CB-XLCT imaging. In this framework, the shrinking for reconstruction is beginning from the whole body without a prior information. Meanwhile, an adaptive multi-scale kernel function is designed to accelerate the shrinking processing. Specially, in each iteration, the kernel function can automatically select the nodes for next round reconstruction, at a multi-scale rate according to the shrinking process. Thus, a feasible solution by the ADS_PSR can be obtained only with a few mesh nodes. Simulation experiment and in vivo experiment were applied to test the performance of the proposed ADS_PSR framework. The widely-used ISPR method, IVTVG, Tikhonov and ART algorithms were employed for comparison. The results demonstrate that all the LE of the simulation case and in vivo case by the proposed framework were all less than 1 mm, and corresponding RQE were almost below 20%. Compared with traditional ISPR method with artificial threshold, our proposed method gets more sparse solution with the fewer iterations. These illustrate the potential of the proposed framework in improving CB-XLCT imaging for feasible and effective reconstruction. Meanwhile, the results of stability case further demonstrated the robustness of the proposed framework for CB-XLCT imaging. Furthermore, the proposed framework can also be extended to other optical molecular modalities, such as FMT and BLT.

However, it should be noted that the size and shape of nanophosphor target cannot be accurately recovered in this paper. Besides, the selection of different number of nodes for reconstruction are often based on experience. Thus, we will focus on these challenges and conduct further research. In conclusion, the proposed framework has a better performance in terms of accuracy, stability and practicability. We hope the proposed framework can facilitate the development of CB-XLCT, as well as other optical molecular tomography technologies.

Funding

National Natural Science Foundation of China (61902317); Science and Technology Plan Program in Shaanxi Province of China (2019JQ-166).

Acknowledgements

The authors would like to thank the School of Life Science and Technology of Xidian University for providing experimental data acquisition system.

Disclosures

The authors declare no conflicts of interest.

References

1. G. Pratx, C. M. Carpenter, C. Sun, and L. Xing, “X-ray luminescence computed tomography via selective excitation: a feasibility study,” IEEE T. Med. Imaging 29(12), 1992–1999 (2010). [CrossRef]  

2. G. Pratx, C. M. Carpenter, C. Sun, R. P. Rao, and L. Xing, “Tomographic molecular imaging of x-ray-excitable nanoparticles,” Opt. Lett. 35(20), 3345–3347 (2010). [CrossRef]  

3. M. Ahmad, G. Pratx, M. Bazalova, and L. Xing, “X-ray luminescence and x-ray fluorescence computed tomography: new molecular imaging modalities,” IEEE Access 2(2), 1051–1061 (2014). [CrossRef]  

4. 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 T. Med. Imaging 37(10), 2176–2184 (2018). [CrossRef]  

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

6. X. Liu, Q. Liao, and H. Wang, “In vivo x-ray luminescence tomographic imaging with single-view data,” Opt. Lett. 38(22), 4530–4533 (2013). [CrossRef]  

7. P. Gao, J. Rong, H. Pu, T. Liu, W. Zhang, X. Zhang, and H. Lu, “Sparse view cone beam x-ray luminescence tomography based on truncated singular value decomposition,” Opt. Express 26(18), 23233–23250 (2018). [CrossRef]  

8. G. Zhang, F. Liu, J. Liu, J. Luo, Y. Xie, J. Bai, and L. Xing, “Cone beam x-ray luminescence computed tomography based on bayesian method,” IEEE T. Med. Imaging 36(1), 225–235 (2017). [CrossRef]  

9. C. Li, K. Di, J. Bec, and S. R. Cherry, “X-ray luminescence optical tomography imaging: experimental studies,” Opt. Lett. 38(13), 2339–2341 (2013). [CrossRef]  

10. C. M. Carpenter, G. Pratx, C. Sun, and L. Xing, “Limited-angle x-ray luminescence tomography: methodology and feasibility study,” Phys. Med. Biol. 56(12), 3487–3502 (2011). [CrossRef]  

11. X. Liu, Q. Liao, and H. Wang, “Fast x-ray luminescence computed tomography imaging,” IEEE Trans. Biomed. Eng. 61(6), 1621–1627 (2014). [CrossRef]  

12. D. Chen, S. Zhu, H. Yi, X. Zhang, D. Chen, J. Liang, and J. Tian, “Cone beam x-ray luminescence computed tomography: a feasibility study,” Med. Phys. 40(3), 031111 (2013). [CrossRef]  

13. D. Chen, S. Zhu, X. Chen, T. Chao, X. Cao, F. Zhao, L. Huang, and J. Liang, “Quantitative cone beam x-ray luminescence tomography/x-ray computed tomography imaging,” Appl. Phys. Lett. 105(19), 191104 (2014). [CrossRef]  

14. X. Liu, H. Wang, M. Xu, S. Nie, and H. Lu, “A wavelet-based single-view reconstruction approach for cone beam x-ray luminescence tomography imaging,” Biomed. Opt. Express 5(11), 3848 (2014). [CrossRef]  

15. Y. Tan and H. Jiang, “Dot guided fluorescence molecular tomography of arbitrarily shaped objects,” Med. Phys. 35(12), 5703–5707 (2008). [CrossRef]  

16. F. Stuker, C. Baltes, K. Dikaiou, D. Vats, L. Carrara, E. Charbon, J. Ripoll, and M. Rudin, “Hybrid small animal imaging system combining magnetic resonance imaging with fluorescence tomography using single photon avalanche diode detectors,” IEEE T. Med. Imaging 30(6), 1265–1273 (2011). [CrossRef]  

17. A. Ale, V. Ermolayev, E. Herzog, C. Cohrs, M. H. 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]  

18. J. Zhang, J. Shi, X. Cao, F. Liu, J. Bai, and J. Luo, “Fast reconstruction of fluorescence molecular tomography via a permissible region extraction strategy,” J. Opt. Soc. Am. A 31(8), 1886–1894 (2014). [CrossRef]  

19. M. A. Naser and M. S. Patterson, “Bioluminescence tomography using eigenvectors expansion and iterative solution for the optimized permissible source region,” Biomed. Opt. Express 2(11), 3179–3193 (2011). [CrossRef]  

20. J. Feng, K. Jia, G. Yan, S. Zhu, C. Qin, Y. Lv, and J. Tian, “An optimal permissible source region strategy for multispectral bioluminescence tomography,” Opt. Express 16(20), 15640–15654 (2008). [CrossRef]  

21. X. He, F. Dong, J. Yu, H. Guo, and Y. Hou, “Reconstruction algorithm for fluorescence molecular tomography using sorted L-one penalized estimation,” J. Opt. Soc. Am. A 32(11), 1928–1935 (2015). [CrossRef]  

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

23. H. Yi, P. Jiao, X. Li, J. Peng, and X. He, “Three-way decision based reconstruction frame for fluorescence molecular tomography,” J. Opt. Soc. Am. A 35(11), 1814–1822 (2018). [CrossRef]  

24. P. Svenmarker, C. T. Xu, H. Liu, X. Wu, and S. Andersson-Engels, “Multispectral guided fluorescence diffuse optical tomography using upconverting nanoparticles,” Appl. Phys. Lett. 104(7), 073703 (2014). [CrossRef]  

25. C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59(1), R1–R64 (2014). [CrossRef]  

26. J. Feng, C. Qin, K. Jia, S. Zhu, X. Yang, and J. Tian, “Bioluminescence tomography imaging in vivo: recent advances,” IEEE J. Select. Topics Quantum Electron. 18(4), 1394–1402 (2012). [CrossRef]  

27. A. D. Klose, B. J. Beattie, H. Dehghani, L. Vider, C. Le, V. Ponomarev, and R. Blasberg, “In vivo bioluminescence tomography with a blocking-off finite-difference sp3 method and mri/ct coregistration,” Med. Phys. 37(1), 329–338 (2010). [CrossRef]  

28. A. D. Klose and E. W. Larsen, “Light transport in biological tissue based on the simplified spherical harmonics equations,” J. Comput. Phys. 220(1), 441–470 (2006). [CrossRef]  

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

30. H. Zhang, G. Geng, S. Zhang, K. Li, C. Liu, Y. Hou, and X. He, “Sparse non-convex lp regularization for cone-beam x-ray luminescence computed tomography,” J. Mod. Opt. 65(19), 2220–2230 (2018). [CrossRef]  

31. M. Chen, H. Su, Y. Zhou, C. Cai, D. Zhang, and J. Luo, “Automatic selection of regularization parameters for dynamic fluorescence molecular tomography: a comparison of l-curve and u-curve methods,” Biomed. Opt. Express 7(12), 5021 (2016). [CrossRef]  

32. G. T. Herman and L. B. Meyer, “Algebraic reconstruction techniques can be made computationally efficient (positron emission tomography application),” IEEE T. Med. Imaging 12(3), 600–609 (1993). [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,” IEEE T. Med. Imaging 18(24), 24825–24841 (2010). [CrossRef]  

34. X. Cao, B. Zhang, X. Wang, F. Liu, K. Luo, and J. Bai, “An adaptive tikhonov regularization method for fluorescence molecular tomography,” Med. Biol. Eng. Comput. 51(8), 849–858 (2013). [CrossRef]  

35. P. C. Hansen and D. P. O’Leary, “The use of the l-curve in the regularization of discrete ill-posed problems,” SIAM J. Sci. Comput. 14(6), 1487–1503 (1993). [CrossRef]  

36. P. C. Hansen, “Analysis of discrete ill-posed problems by means of the l-curve,” SIAM Rev. 34(4), 561–580 (1992). [CrossRef]  

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

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

39. D. Hyde, R. Schulz, D. Brooks, E. Miller, and V. Ntziachristos, “Performance dependence of hybrid x-ray computed tomography/fluorescence molecular tomography on the optical forward problem,” J. Opt. Soc. Am. A 26(4), 919–923 (2009). [CrossRef]  

40. S. Ren, X. Chen, H. Wang, X. Qu, G. Wang, J. Liang, and J. Tian, “Molecular Optical Simulation Environment (MOSE): A Platform for the Simulation of Light Propagation in Turbid Media,” PLoS One 8(4), e61304 (2013). [CrossRef]  

41. L. A. Feldkamp, L. C. Davis, and J. W. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1(6), 612–619 (1984). [CrossRef]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)

Fig. 1.
Fig. 1. Schematic diagram of CB-XLCT system [1112].
Fig. 2.
Fig. 2. The multi-scale of $\zeta$ function with the increase of $k$ .
Fig. 3.
Fig. 3. Reconstruction model in simulation experiment. (a) 3D Mouse Model (b) Forward simulation of single-target case (c) Forward simulation of double-targets case (d) 2D view of CT slice for the single-target case (e) 2D view of CT slice for the double-targets case.
Fig. 4.
Fig. 4. The 2D views ( $Z = 50.8$ mm plane) of the reconstructed results for the single-target reconstruction.
Fig. 5.
Fig. 5. The iteration number of ADS_PSR and ISPR with artificial threshold in single-target case.
Fig. 6.
Fig. 6. 2D views ( $Z = 60$ mm plane) of the reconstructed results for the double-targets reconstruction.
Fig. 7.
Fig. 7. The iteration number of ADS_PSR and ISPR with artificial threshold in double-targets case.
Fig. 8.
Fig. 8. Illustration of LE and RD at different number of view angles
Fig. 9.
Fig. 9. Illustration of LE and RD at different noise levels.
Fig. 10.
Fig. 10. Illustration of LE and RD at different number of nodes.
Fig. 11.
Fig. 11. The micro-CT result of the female mouse.
Fig. 12.
Fig. 12. Reconstruction results overlaid with CT data of in vivo experiment. (a) coronal-view results at $X = 12$ mm. (b) sagittal-view results at $Y = 13.2$ mm; (c) transversal-view results at $Z = 6.8$ mm.
Fig. 13.
Fig. 13. The iteration number of ADS_PSR and ISPR with artificial threshold of in vivo experiment.

Tables (7)

Tables Icon

Table 1. Procedure of the adaptive shrinking permissible source region framework for CB-XLCT imaging

Tables Icon

Table 2. Quantitative results of the single-target case

Tables Icon

Table 3. Quantitative results of the double-targets case

Tables Icon

Table 4. Comparison results for reconstruction single-target by ADS_PSR + IVTCG with different measurements

Tables Icon

Table 5. Comparison results for the reconstruction single-target by ADS_PSR + IVTCG with different noise levels

Tables Icon

Table 6. Comparison results for the reconstruction single-target by ADS_PSR + IVTCG with different numbers of nodes and tetrahedral elements

Tables Icon

Table 7. Quantitative results of the in vivo experiment

Equations (8)

Equations on this page are rendered with MathJax. Learn more.

X ( r ) = X 0 e x p { r 0 r μ x ( τ ) d τ }
S( r ) = ε X ( r ) ρ ( r )
{ [ D ( r ) Φ ( r ) ] + μ a ( r ) Φ ( r ) = S ( r ) r Ω Φ ( r ) + 2 κ D ( r ) [ ν Φ ( r ) ] = 0 r Ω
J = A ρ
min | | A ρ J | | 2 2 + τ | | ρ | | p
P k = { 1 if node i is within the permissible region 0 otherwise
min R k | | A k ρ k J | | 2 2 + τ k | | ρ k | | p
ζ k = 1 β exp [ ( k 1 ) / ( k 1 ) α α ] 1 + β exp [ ( k 1 ) / ( k 1 ) α α ]
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.