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

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


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 [1][2][3]. 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 [6][7][8].
In general, narrow-beam XLCT, pencil-beam XLCT (PB-XLCT) [9], and cone-beam XLCT (CB-XLCT) [1][2][3] 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 [15][16][17][18][19].
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,[20][21][22][23][24]. 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.

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.

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]: 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 µ 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: where S(r) is the target, which is often used to mimic the small tumor. ρ(r) is the nanophosphor concentration at position r, and ε 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 [27][28]: where D(r) = (3(µ a (r) + µ s (r))) −1 is the diffusion coefficient, where µ a (r) and µ s (r) is the absorption and reduced scattering coefficient, respectively. Ω is the domain of the image object and ∂Ω is the corresponding boundary. Φ(r) is the photon fluence at position r. ν is the outward unit normal vector to boundary, and κ 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 ρ and the NIR measurement J as following matrix-form equation: where A is the system matrix which is used to map the unknown ρ to known measurement J [8].

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 ρ 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 [12][13][14]. Herein, the compressed sensing (CS) theory can be applied to solve ρ by converting Eq. (4) into the following optimization problem [30]. min ||Aρ − J|| 2 2 + τ|| ρ|| p (5) where τ is the regularization coefficient. || ρ|| p is the L p -norm of ρ. 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].

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: As soon as P k is determined, A and ρ in Eq.
respectively. ⊗ 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: min 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 ρ k without manual intervention, we designed an adaptive kernel function to automatically speed up the shrinking process at a multi-scale rate as follows: where k is the iterations which is set from 1. α and β is multi-scale factor and nonnegative coefficient, respectively. α and β are all determined experientially to confine the shrinking rate. For the reconstruction, the range of α and β is [2.5, 3.5] and [3.5, 7.5], respectively. Figure 1 shows the value of ζ function with the increase of k (α = 3, β = 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 Num(R) · ζ k and is the integer operation. With the increase of K, the trend of solution ρ(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 ρ 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 ρ can be acquired. Thus, the function ζ k can adaptively fit this case with a multi-scale rate. Table 1 shows details of the proposed framework for CB-XLCT imaging. Algorithm: ADS_PSR framework for CB-XLCT imaging Step 1: Parameter initialization. Set the iteration number k = 1 and compute the multi-scale rate ζ as in Eq. (8) if k>1. Set the final nodes N f = 10 according to experience.
Step 2: Solving Eq. (5) for Global reconstruction to determine a current PR.
Step 3: Sort recovered density of current nodes in the decreasing order, choose Num(R) · ζ k mesh nodes with higher recovered density to form a new PR.
Step 4: Get the vector P k by Eq. (6) and Computing the optimized system matrix A k .
Step 5: Local reconstruction based on Eq. (7) to refine a solution.
Step 6: if |A k+1 ρ k+1 − A k ρ k |/ |A k ρ k | ≤ ξ or the number of nodes of refined PR N k+1 ≤ N f , output the solution ρ based on the final PSR and quit; otherwise, go to step 2 for next round iteration.

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 [19][20][21][22][23].
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 [1][2][3], 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 [22][23]. In addition, the regularization coefficient τ of IVTCG and Tikhonov is automatically selected by the L-curve method [35][36].
The multi-scale factor α and nonnegative coefficient β of ζ 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.

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 µg/mm 3 , and the corresponding light yield were 0.15 cm 3 /mg [12]. 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.

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

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. As can be seen in Figs. 6-7, 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.

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 µg/mm 3 . Herein, it is obvious that the proposed method is robust in different view angles.  It is noted that noise was unavoidable in CB-XLCT imaging [1][2][3], 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 µg/mm 3 . It is clear that the proposed method is robust to measurement noise. 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

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 × 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. 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 L 2 -norm is not sparse, and the ART fails to obtain feasible reconstruction results.

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