X-ray luminescence computed tomography imaging based on X-ray distribution model and adaptively split Bregman method.

X-ray luminescence computed tomography (XLCT) has become a promising imaging technology for biological application based on phosphor nanoparticles. There are mainly three kinds of XLCT imaging systems: pencil beam XLCT, narrow beam XLCT and cone beam XLCT. Narrow beam XLCT can be regarded as a balance between the pencil beam mode and the cone-beam mode in terms of imaging efficiency and image quality. The collimated X-ray beams are assumed to be parallel ones in the traditional narrow beam XLCT. However, we observe that the cone beam X-rays are collimated into X-ray beams with fan-shaped broadening instead of parallel ones in our prototype narrow beam XLCT. Hence we incorporate the distribution of the X-ray beams in the physical model and collected the optical data from only two perpendicular directions to further speed up the scanning time. Meanwhile we propose a depth related adaptive regularized split Bregman (DARSB) method in reconstruction. The simulation experiments show that the proposed physical model and method can achieve better results in the location error, dice coefficient, mean square error and the intensity error than the traditional split Bregman method and validate the feasibility of method. The phantom experiment can obtain the location error less than 1.1 mm and validate that the incorporation of fan-shaped X-ray beams in our model can achieve better results than the parallel X-rays.


Introduction
Since the 1970s, the rare earth phosphors, especially the oxysulfide phosphors, due to its advantages such as high luminescence efficiency and innocuity, etc., have caused great attention as the host materials of X-ray phosphors [1]. While radioluminescence of phosphor material has long been used in radiation detectors, the use of nanophosphors in biological contexts is just beginning to be explored [2]. Among these applications, X-ray luminescence computed tomography (XLCT) is one of the newly-developed and most attractive imaging modality, which has been demonstrated by in vitro imaging of nanoparticles in tissue phantoms for oxysulfides (Gd 2 O 2 S: Eu) [3]. The combination of insignificant scattering of X-rays in tissues, and the high tissue penetration of near-infrared (NIR) optical photons emitted from appropriate phosphors, opens up the possibility of achieving deep tissue optical imaging in vivo with unprecedented spatial resolution [4]. Meanwhile researchers have synthesized potentially less toxic particles recently [5,6]. This can overcome the limitation of phosphor nanoparticle synthesis of XLCT molecular probes and further promote its future application.
Based on the promising nanophosphors and its attractive imaging principles, many researchers have studied and developed X-ray luminescence computed tomography (XLCT) these years. To take advantage of linear propagation of the X-rays, Xing et. al. proposed the narrow beam X-ray excitation technology known as XLCT in their papers [7]. Then, Li et. al. used a collimated pencil beam X-ray with a beam size of 1 mm diameter to more selectively excite deep targets [8]. However, since the scanning area of the collimated pencil beam X-ray is smaller than narrow beam XLCT, this technology demands relatively longer scanning time under X-ray exposure which limits its application in biology processes. In order to overcome this issue, a cone beam XLCT imaging system was designed by our group [9] to improve scanning speed, and was applied to small animal experiments by Liu et. al. [10]. Nevertheless, the introduction of cone beam XLCT may cause limited spatial resolution for deep targets due to light scatter, which is similar to fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT) [8]. Wang et. al. also analyzed that the cone beam scanning mode does not sufficiently utilize the primary benefit of XLCT in terms of a reduced permissible source region and presented a fan-beam stimulation scheme for XLCT, which can be regarded as one transformation of narrow beam XLCT [11]. Hence, among these XLCT imaging strategies, narrow beam XLCT can be regarded as an optimal balance between the pencil beam mode and the cone-beam mode in terms of imaging efficiency and image quality.
In the traditional narrow beam XLCT, the collimated X-ray beams are modeled as parallel ones [7]. However, in our prototype narrow beam XLCT system, we observed that the X-rays were collimated into fan-shaped X-ray beams instead of parallel X-rays, which was different from traditional narrow beam XLCT. The description of precise X-rays distribution can contribute to the XLCT imaging quality. Hence we first established the X-ray distribution model through our experiment and applied it in our propagation model. Meanwhile to further speed up the scanning time, we collected the optical data from only two perpendicular directions. In order to recover the distribution of the nanophosphors within such limited information, we also proposed a depth related adaptive regularized split Bregman (DARSB) method to overcome the ill-posedness of our reconstruction problem.
The remaining sections are organized as follows. Section 2 describes the photon propagation model and reconstruction method in our system. In Section 3, the numerical simulation with different phantoms is conducted as well as the phantom experiment. Finally, the conclusion and discussion are given in Section 4. A prototype XLCT system is built in our laboratory with the configuration shown in Fig. 1. The system can perform conventional computed tomography (CT) imaging with a microfocus X-ray source (Apogee, Oxford Instruments, USA) and an X-ray flat panel detector (C7921CA-02, Hamamatsu, Japan). Meanwhile, the X-ray beam is collimated with the X-ray collimator made up of two L-shaped lead shields to excite the nanophosphors. A liquid-cooled back illuminated CCD camera (PIXIS 2048B, Princeton Instruments, USA) with a focus lens (Micro-Nikkor 55 mm f/2.8, Nikon, Japan) is positioned to capture the emitted light for the X-ray luminescence imaging. The phantom is placed on a motorized rotation stage which is mounted on a motorized linear stage to realize the scanning of the object with selectable step interval at selectable directions. A lead made X-ray shield is placed between the camera and the phantom to reduce x-ray radiation on the CCD.

Imaging System
It is necessary to perform a fundamental assessment after building our experiment system. We collect the X-ray data from the X-ray detector with the collimator positioned in front of the X-ray source. First we adjusted the collimator with the motorized linear stage to align the X-ray beam with the X-ray source focus spot, which was realized by comparing the maximum of the X-ray data on the X-ray detector for every movement. Then we adjusted the X-ray detector with different distances from the X-ray source (715 cm, 690 cm, 665 cm, and 640 cm) and obtained corresponding X-ray data. The width of X-ray beam was measured by selecting the full width at half maxima (FWHM) of the average value of the columns in the X-ray data and multiplying the pixel of the detector as shown in Fig. 2(a). The widths of the X-ray were 1.275, 1.2, 1.125, and 1.05 mm at the corresponding position. The experiment shows that the width of the X-ray beam is linearly changing with the position of the detector as shown in Fig. 2(b). It can be concluded that the profile of the X-ray after the collimator is fan-shaped broadening at the plane of the stage. Meanwhile, based on the relationship between the X-ray beam width and the distances from the X-ray source, the X-ray beam was 0.6 mm at the end of the collimator. Based on the X-ray distribution model and the widths of X-ray beams obtained on the X-ray detector, we can calculate the width of the X-ray beam at any distance between the collimator and the X-ray detector, which can be applied to establish the accurate fan-shaped X-ray distribution model in the following photon propagation model.

Photon propagation model
In XLCT imaging process, the imaging model involves the following three steps: First, the Xrays are emitted from the X-ray source and suffer from attenuation through the tissues or the phantoms as follows: where X 0 is the X-ray source intensity with the initial position r 0 , and µ t (τ) is the X-ray attenuation coefficient at position τ; Second, once the nanophosphors are excited by the X-rays, they will emit NIR light, which follows linear relationship with the nanophosphor density and the X-ray intensity as follows: where S(r) is the light source, X(r) is the X-ray intensity, ρ(r) is the nanophosphor density at position r and ε is the light yield; Last, the NIR light will transport to the surface of the biological tissues and be captured by our CCD camera which can be accurately modeled by the radiative transfer equation (RTE) or diffusion approximation [12] and expressed as: where r is the position vector, Ω is the domain under consideration, D(r) is the diffusion coefficient, Φ(r) is the photon flux density and S(r) is the source. Based on the finite element method (FEM), the above model can be converted into the following matrix between the nanophosphor concentration ρ and the NIR measurement Φ, whose details can be found in [9]: where A = A Ψ, ρ = ρ Ψ. A is the system matrix constructed with FEM, represents the operation that removes rows in matrix A corresponding with the zero-elements in columns vector Ψ while Ψ is a columns vector composed of zero or unity elements and defined by the specific X-rays distribution as follows: if node n is within the X-ray scanning area 0 otherwise , Our above fundamental assessment convealed that the cone beam X-rays became X-ray beams with fan-shaped broadening after collimation. Based on the X-ray distribution model and the widths of X-ray beams obtained on the X-ray detector, we can establish the accurate fan-shaped distribution model using Eq. (5).

Reconstruction Method
The solution to the above linear inverse problem in Eq. (4) can be tackled as the recovery of the nanophosphor concentration ρ from the measured luminescent boundary data Φ. The inverse reconstruction problem suffers from ill-posedness because of the high scattering of light in biological tissues, the limited optical measurements and the other experiment conditions. This ill-posed problem can be addressed by means of implicit regularization and skillful iterative methods like the algebraic reconstruction technique [13], unconstrained optimization methods based on l 2 -regularization [14] or l 1 -regularization [15] and even l p -regularization [16]. Meanwhile, the priori information such as the permissible region strategy [17] and multispectral method [18] are also commonly applied in optical imaging and have also been applied in XLCT [19]. Compared with these advanced methods, the methods based on Bregman iteration transform the optimization problem to an equivalent constrained problem with the concept of Bregman distance [20]. Among these methods, the split Bregman (SB) method addresses the optimization problem by analyzing it into several functions and minimizing them separately on an efficient and simple way [21,22]. For instance, for l 1 -constraints such as the total variation (TV), the SB method decouples l 1 -functionals and l 2 -functionals, and minimizes them separately, the l 2 -part by using conventional methods, and the l 1 -part by straightforward shrinkage formulas [23]. The SB method has been widely applied in image denoising and magnetic resonance imaging [20,21,24]. Recently, it has also been extended to the application of fluorescence diffuse optical tomography and bioluminescence tomography [23,25]. The problem is converted into the following unconstrained optimization problem by the Lagrangian method [26]: In reconstruction, there are some problems that the reconstructed sources tend to become superficial while the localization errors are increasing particularly for deep sources [27]. Here we fully utilized the X-ray attenuation information as a related depth weighting factor and introduced two auxiliary variables when referring to the variable splitting technique. Hence, we apply the X-ray attenuation information as a related depth weighting factor to modify the l 1 -norm part of Eq. (6) as follows: whereρ = Π · * ρ in which Π is the vector that reflects the corresponding X-ray relative attenuation value at each vertex and · * represents the product of corresponding elements of the vectors Π and ρ. The corresponding X-ray relative attenuation value at each vertex is estimated by accumulating the attenuation value on the voxel data of CT reconstruction. Then we decouple it into two subproblems to minimize over ρ with d k fixed and next minimize over d with ρ k+1 fixed alternately [28], in which λ and θ are the regularization coefficient and "splitting" regularization coefficient respectively.
Based on the present SB scheme, we proposed a depth related adaptive regularized split Bregman (DARSB) method to overcome the ill-posedness of the inverse problem in Eq. (4).
Meanwhile, our method involves on the adaptive regularization parameter estimation with adjusting it in a closed form in each iteration, where Morozovs discrepancy principle is applied to provide a modular solver to update the nanophosphor density ρ [29].

Deduction of the proposed method
An auxiliary variable x and another auxiliary variable y are applied to represent Aρ and ρ respectively. Then, the Eq. (7) can be transformed into: With the Bregman function J(ρ, x, y) and the corresponding Bregman distance the following iterative framework is deduced based on the SB method and the alternating direction method (ADM) [21,29]: x k+1 = arg min Hence, Eq. (7) can be transformed into the problem of tackling the above subproblem (11)-(14) efficiently. The subproblem with respect to ρ is solved by preconditioned conjugate gradient method with the backtracking operation to guarantees the global convergence [30]. Meanwhile we apply the nonnegative constraint with respect to ρ in subproblem (11). Eq. (12) can be solved by shrinkage as follows: where For the minimization problem concerning x k+1 , we have to analyze the problem that whether x − Aρ 2 2 ≤ c to satisfy the discrepancy principle, which is realized by the following two situations in each iteration [29]:

Parameter selection
During the above algorithm, the choice of the noise dependent upper bound c is still undetermined. This parameter can reflect and constrain the error between the actual solution and the reconstructed one to a reasonable level, which is normally proportional to the noise variance and still an open problem [29]. Here we estimate it by referring to the minimum of the Aρ −Φ 2 2 iteratively. Meanwhile, we introduce two parameters β 1 and β 2 when calculating the regularization parameter adaptively. It is obvious that β 2 is the weight of the difference between y and ρ while β 1 is the weight of the difference between x and Aρ. We set β 1 = 10 10 * SNR × β 2 , β = 10 −3 , where the signal-noise ratio (SNR) of the optical signal is obtained by setting the maximum of the local variance as the variance of the signal while the minimum as the variance of the noise [31]. Even though the convergence of the proposed method can be guaranteed by the condition that β 1 and β 2 > 0 [32], a relatively higher SNR represents the nearer distance between Aρ and φ , which results in a relatively larger β 1 to control the variable x. All the settings of these parameters are obtained empirically.

Evaluation standard
In order to value our reconstruction results more clearly and further explain the validity of our proposed method, we propose two benchmarks that are regularly applied in optical imaging. The first one is the location error of the Euclidean distance between the actual source position and the reconstructed source position as follows [33]: where (x, y, z) is the center coordinate of the reconstruction source with the maximum density and (x 0 , y 0 , z 0 ) is that of the actual source center. We also introduce DC (dice coefficient) [34] to compare the similarity of the reconstructed results and the actual source distribution, which is obtained as follows: where |RS| and |AS| represent the number of tetrahedral-elements in reconstructed results and the actual source correspondingly, while |RS AS| is the number of tetrahedral-elements shared by reconstructed results and the actual source. The larger value of DC can reflect that the reconstructed results are more similar to the actual source.
To further evaluate the accuracy of the reconstruction, the MSE (mean square errors) is applied in our experiments as follows: where ρ t k is the true density while ρ r k is the reconstructed density on the kth mesh element. N is the number of elements where ρ t k > 0. Then the relative intensity error between the actual source and the reconstructed one is applied and defined as [11]: where ρ t k is the true density while ρ r k is the reconstructed density on the kth mesh element respectively, and Num(ρ t k > 0) is the number of elements where ρ t k > 0. The relative intensity error can accurately reveal the intensity error between the actual source and the reconstructed one in the actual region and can be seen as a combination of the relative error of the power density and the measurement of overlapped region, both of which are used to reflect the real reconstruction performance in [35]. The smaller intensity error can indicate the better quantitative reconstruction results.

Simulation study
Based on the above system test, we conducted simulation study and compared our method with the traditional SB method to further validate the feasibility and effectiveness of our proposed method. We used four different cylinder phantoms (phantom A, B, C and D) in the simulation study. The phantom with the size of a 30 mm diameter and 30 mm height contained a small cylinder with a 4 mm diameter and 4 mm height to represent the nanophosphors, whose centers were set to (15,15,20) mm, (15, 11.25, 20) mm, (15, 7.5, 20) and (15, 3.75, 20) mm respectively, as shown in Figs 3 (a), (b), (c) and (d). The cylinder phantom was applied to mimic the muscle tissue, whose absorption and reduced scattering coefficients of the phantom were set to 0.013 mm −1 and 0.93 mm −1 respectively [36]. The X-ray attenuation coefficient was set to 0.0475 mm −1 [37]. The light yield ε was set to 0.15 cm 3 /mg [38]. We applied the distribution of the X-ray beam in the above X-ray beam assessment as well as the distances of the X-ray source, X-ray detector and the stage to ensure the validity of the simulation. The phantom was supposed to move with step size of the width of the X-ray beam. The simulation data was acquired by employing the light yield ε, the optical parameters, X-ray attenuation coefficient, the discretized information of the phantoms and the system matrix A calculated by the Eq. (4). It is reported that measurements at two orthogonal projections are enough to reconstruct an X-ray luminescence optical tomography (XLOT) image of a single target [39] or even more targets with certain strategy. Hence the simulation was conducted with the assumption of scanning the phantom at two perpendicular directions presented by the black arrows in Fig. 3. During the simulation, we obtained 9 projections by moving the phantom 9 times with the step size of the width of the X-ray beam for each direction. The SB method was applied in the l 1 -regularized problem of Eq. (6) and the parameters applied in the method are regularization coefficient and "splitting" regularization, which were empirically set as 2 × 10 −2 and 10 −6 correspondingly [21].
The reconstruction results were obtained as shown in Fig. 3. Figures 3(a) Table 1. For each phantom, the proposed method has a better performance on location error, DC, MSE and intensity error than SB method. Moreover, the proposed method can better reveal the contour of the nanophorsphors than the traditional split method as shown in Fig. 3, which can be also seen in DC results in Table 1. This is mainly caused by the introduction of the weighting scheme and the adaptive regularization parameter estimation, which overcome the ill-posed problem in the reconstruction. Meanwhile, as shown in Table 1, the reconstruction results including location error, DC, MSE and intensity error become worse when nanophosphors locate deeper in the phantom. This can be explained by the facts that more energy of the X-rays is lost when the depth of the nanophosphors is increasing, which means that less energy of the X-rays is used to excite the nanophosphors.
Furthermore, we applied the reconstruction under different noise level to test the robustness of the proposed method. The noise of different levels was added which followed the normal  distribution with a mean of 0 and the standard deviation varied from 10% to 40% of the average of the surface measurements with a 10% interval. The reconstruction location errors remain the same under all of the noise levels in each phantom. Meanwhile, there is a little difference in the DC, MSE and intensity error as shown in Table 2. The DC, MSE and intensity errors stay the same level in each phantom, which can indicate that the measurement noise affect the imaging quality with little fluctuation. Hence it can be induced that the proposed method is robust to measurement noise.

Phantom experiment
Based on the above simulation study, we conducted a physical phantom experiment in our imaging system to further evaluate our proposed imaging strategy. A polyoxymethylene made cylinder phantom with a 20 mm diameter and 20 mm height was applied to mimic biological tissues. The absorption and reduced scattering coefficients of the phantom were set to 0.025 mm −1 and 1.12 mm −1 respectively, which were measured by diffuse optical tomography. The nanophosphors were put into a plastic capillary with a 1 mm radius and 3 mm height, which was embedded in the phantom. During the experiment, we used the X-ray source to excite the nanophosphors from two perpendicular directions and the CCD camera to collect the luminescent photons emitted from the phantom. The image acquisition system was enclosed in a light-tight environment to avoid the outside light effect. Before the luminescence signal collection, we measured the X-ray beams with the same methods applied in the above imaging system test to obtain the width of the Xray beam formed by the collimator. The width of the X-ray beam was 0.5 mm at the end of the collimator. During the luminescence signal collection, the phantom was moved with the step size of the X-ray beam. The exposure time was set to 2 s. It only took about 56 s to complete the data collection without considering the time of phantom movement. The fusion images of the collected data from one direction are shown in Fig. 4. Figures 4(a)-4(n) show the fusion images of the phantom when it moves at the direction parallel to the X-ray source.It can be observed that with the movement of the phantom, the collected data have changed obviously. If the center of the nanophosphors was near the X-rays beams, the collected data was becoming larger because more nanophosphors were excited by the X-rays. Then, the computed tomography (CT) projections were first performed (50 kVp, 1.0 mA, 360 views with 1 • intervals) and were reconstructed by the filtered back-projection method to get the physical structure and the corresponding X-ray attenuation coefficient of the phantom. The center of the capillary in the phantom was (42.67, 45.81, 21.79) mm, which was obtained from the CT reconstruction results. In the inverse reconstruction, the cylinder phantom was discretized into 56621 tetrahedral-elements and 11009 nodes from the CT results by AMIRA. From the measured data, the distribution of the luminescent sample could be reconstructed by the proposed method.
In the experiment, in order to further validate the accuracy of the proposed XLCT imaging, we compared the results in which the X-ray beams were approximate as parallel ones. We also compared the results of the proposed method with the traditional SB method. Figure 5 shows the results in the slices of z = 21.79 mm and x = 42. 67 mm respectively represented by red dotted lines shown in Fig. 5(a). Figure 5(a) shows the micro-CT result on the x-z plane.  Table 3. It can be observed that the results of the fan-shaped X-rays model were better than the results of the parallel X-rays model, which means that the model of the accurate attenuation of the X-rays can help improve the reconstruction effectively. Meanwhile, the results from two different methods were consistent with the simulation experiments. The reconstruction results from the proposed method have better performance than SB method in the location error. Moreover, the proposed method can better reveal the contour of the nanophorsphors than the traditional split method as shown in Fig. 5. Moreover, the reconstruction results with the parallel X-rays model are sparser than that with fan-shaped X-rays model as shown in Fig. 5. The parallel X-rays model is less accurate than the fan-shaped X-rays one, and the scanning area is disagree with actual experiment system. This may cause the sparser reconstruction results in the parallel X-rays model.  The results show that the distribution of the nanophosphors in the phantom can be revealed by our proposed method in the XLCT imaging system. However, due to the ill-posed problem in the phantom reconstruction and the measurement errors in the process of the data collection, the reconstruction error was slightly bigger than in the simulation experiments, which was normal in the phantom experiment. The light yield ε of the nanophosphors need related information of the X-ray energy, which couldn't be directly measured in our present system conditions. Hence,  the unit of Fig. 5 was ignored and the comparison of quantity results as well as MSE results was not discussed and will be studied in future.

Conclusion
In this paper, a narrow beam XLCT imaging modality is proposed, in which the fan-shaped Xrays distribution is modeled and a depth related adaptive regularized split Bregman (DARSB) method is proposed to overcome the ill-posedness of reconstruction. The simulation and phantom experiments were conducted in our imaging system, which could perform conventional CT imaging and X-ray luminescence imaging. The simulation experiments show that the proposed method can achieve better results both in the location error, DC, MSE and the intensity error than the traditional SB method. The simulation experiments showed that our proposed method was stable and robust against different measurement noise levels. The phantom experiment can obtain the location error of 1.02 mm and further validate that the incorporation of fan-shaped X-ray beams in our model can achieve better results than the parallel X-rays.
The cone beam XLCT can be viewed as one of the optical imaging modality such as bioluminescence tomography (BLT) and fluorescence molecular tomography (FMT), whose reconstruction quality are relatively lower than that of narrow beam XLCT. However, the imaging quality and scanning time of the narrow beam XLCT is closely related to the beam width and sampling [40]. For example, C. M. Carpenter et. al. reported that objects as small as 1 mm in size could be resolved using a 1 mm beam width [40] and that the phantom was translated 26 times in increments of 1 mm and was rotated 24 times to cover 360 • [7]. However, the scanning time in our imaging system could be shortened with the relatively larger scanning area of X-rays. Meanwhile, the scanning only performed from two perpendicular directions to accomplish the reconstruction. Even though the scanning time in our phantom experiment was shorter than the cone beam XLCT imaging [9], it cannot indicate that the scanning time in our system is less than the cone beam XLCT because different nanophosphors were applied and the scanning time is related to the size of the object. Nonetheless, compared with the traditional narrow beam XLCT imaging, the scanning time in our proposed imaging modality is greatly decreased. However, the effects of the width of the X-ray beam on the number of projections as well as the results still need further study and research.