Spectral-resolved cone-beam X-ray luminescence computed tomography with principle component analysis

: Cone-beam X-ray luminescence computed tomography (CB-XLCT) has become a promising technique for its higher utilization of X-ray and shorter scanning time compared to the narrow-beam XLCT, but it suffers from the low-spatial resolution that results in the insufficiency to resolve the adjacent multiple probes. In multispectral CB-XLCT, multiple probes show different emission behaviors in the dimension of the spectrum. In this work, a spectral-resolved CB-XLCT method combining multispectral CB-XLCT with principle component analysis (PCA) was proposed to improve the imaging resolution. Results of digital simulation and the phantom experiment illustrated that the proposed method was capable of resolving adjacent multiple probes accurately and had better performance than the common multispectral CB-XLCT with spectrum information priori .


Introduction
The feasibility of X-ray luminescence computed tomography (XLCT) was first demonstrated by Xing's group in narrow-beam X-ray mode [1]. The principle of XLCT is that the nanoparticles can emit visible or near-infrared (NIR) luminescence when excited by X-ray beams and then the distribution of nanoparticles can be reconstructed from the projection data of luminescence [2,3]. Compared to the regular optical molecular tomography, the main advantages of XLCT are the increased excitation depth in tissue and the elimination of tissue autofluorescence. The narrow-beam XLCT showed good spatial resolution, but it suffered from the long scanning time which hindered its practical application to some extent. For fast imaging, different excitation modes of X-ray were studied. Cong et al. used fan-beam X-ray to excite nanoparticles in XLCT and can reduce scanning time, but it still need slice-by-slice scanning [4]. The cone-beam XLCT (CB-XLCT) was first proposed by Chen et al. [5], and applied to small animal experiment by Liu et al. [6]. Compared to the narrow-beam and fan-beam mode, CB-XLCT can illuminate the whole body and stimulate all the nanoparticles simultaneously, so it can improve the X-ray utilization and reduce the scanning time greatly which is critical for clinical application. However, CB-XLCT suffered from the ill-posed reconstruction which resulted in low-spatial resolution, and the ill-posed problem became worse for in vivo imaging due to the high scattering characteristic of photons in biological tissue. As a result, it is a challenge for CB-XLCT to resolve multiple probes especially when they are adjacent.
Efforts have been done to resolve the multiple probes in CB-XLCT recently. Multivoltage excitation scheme were used to resolve multiple targets of different kinds by Liu et al. [7] and of different concentrations in our previous work in CB-XLCT [8]. Compared to multivoltage excitation scheme, multispectral scheme could be another optional choice. Multispectral scheme has been applied in optical molecular tomography successfully [9][10][11][12][13] Compared with conventional CB-XLCT, multispectral CB-XLCT cannot only obtain more valid luminescence data but also get the different emission information of the multiple probes in the dimension of spectrum. An exposition of resolving multiple probes can come down to an analysis of the correlative data obtained over spectrums, so multivariate analysis can be applied to multispectral CB-XLCT for the purpose. . In optical molecular imaging, PCA was used to analyze optical data varying with times [22,23], spectrums [12] and energies [7,8]. Despite of these successful applications of PCA, no trials have been reported to use PCA in multispectral CB-XLCT. As PCA is a data-driven method and CB-XLCT has its own characteristics of poor spatial resolution, investigation on the application of PCA in multispectral CB-XLCT is very important.
In this work, a new method called spectral-resolved CB-XLCT by combining multispectral CB-XLCT with PCA was proposed to resolve adjacent multiple probes. Digital simulation and phantom experiments were performed to test the performance of the proposed method, and multispectral CB-XLCT with spectrum information priori was also studied for comparison. The qualitative and quantitative evaluations suggest that the method of performing PCA on multispectral CB-XLCT images is capable of resolving adjacent multiple probes accurately.
The outline of this paper is as follows. In Section 2, the basic theories of multispectral CB-XLCT with spectrum information priori and the proposed spectral-resolved CB-XLCT are detailed. In Section 3, the methods of experiments and evaluations are described. The results are illustrated and evaluated in Section 4. Finally, we discuss and conclude the major findings of this work.

Forward model
In CB-XLCT, the light propagation mainly contains two parts: propagations of the X-ray to the nanoparticles and the emission light from nanoparticles [5]. Based on Lambert-Beer law, during the propagation course of X-ray through the tissues, the X-ray intensity ( ) χ r along the primary path can be expressed as: where 0 ( ) χ r is the X-ray source intensity at the initial position 0 r , and x μ is the X-ray attenuation coefficient which can be obtained from XCT. The emission light source ( ) s S r excited by X-ray at position s r can be described as follows: where ρ is the nanoparticle concentration and χ is the X-ray intensity. η is the light yield which is the intrinsic property of the nanoparticle and can be experimentally determined. Light yield is defined as the total quantum yield at all emission wavelengths per unit nanoparticle concentration and X-ray intensity usually. The emission light transported in the tissues can be expressed by the diffusion equation (DE) approximately [5]: where ( ) D r is the diffusion coefficient, Φ is the photon density and a μ is the absorption coefficient. By discretizing the reconstruction region to N voxels, the above imaging model can be converted into the following linear relationship based on the finite-element method where x is the unknown N 1 × column vector for reconstruction which is proportional to nanoparticle concentration distribution ρ with the relationship of η = x ρ and the weighted matrix W describes the forward model from the unknown vector x to the known measurements m Φ .
The emission photons at all wavelengths are taken into account to define η in conventional CB-XLCT with no optical filters, and η can be taken as constant for a certain nanoparticle. But in multispectral CB-XLCT, η should be taken as variable. And Eq. (4) can be rewritten as follows:

Multispectral CB-XLCT with spectrum information
In the common multispectral CB-XLCT, the spectral information (SI) about ( ) But for multiple nanoparticles of different kinds, the Eq. (6) is not applicable. The position information priori of multiple nanoparticles are not available in common cases, so light yields cannot be assigned to each kind nanoparticle correspondingly at each wavelength. The reconstruction targets ( ) in the Eq. (5) are the combination of multiple nanoparticles, at wavelength k λ can be expressed as the following: where P is the number of nanoparticle kinds and P ρ is distribution of the P-th kind nanoparticle. Then ( ) k λ x for all wavelengths can be denoted as: To get P ρ from Eq. (8), it is critical that number of wavelengths M should satisfy the condition M P ≥ . In the work, the number of wavelengths taken into account was set to be equal to P for simplification. A can be obtained by spectral analysis experiments as spectrum information priori, and then the distribution of each kind nanoparticle can be recovered from Eq. (8).

Multispectral CB-XLCT with PCA
Another way to obtain the individual distributions of multiple nanoparticles from ( ) k λ x data is using an unmixing technique [27]. Multivariate analysis such as PCA can be used for the purpose. The reconstructed results in multispectral CB-XLCT were assembled as the input matrix for PCA, which can be denoted as . X is normalized to the maximal value in this work. After PCA, a N M × matrix P of the principle components (PCs) can be obtained as follows [12]: where 0 X is the centered data of X by subtracting its column means and the j-th principle component which is the j-th column of P can be denoted as: where j e is the j-th column vector of E . The absolute values of positive and negative elements in PCs are further utilized to create images respectively, from which the spatial distribution of multiple nanoparticles can be obtained. The function (princomp) in MATLAB software (MathWorks, Natick, MA, USA) was used to perform PCA in this work. The methods of conventional CB-XLCT with no optical filters, the multispectral CB-XLCT with spectrum information priori and the multispectral CB-XLCT with PCA are denoted hereafter as C-NF, M-SI and M-PCA methods, respectively.

Digital simulation
The model for digital simulation was composed of a cylinder containing two spheres, as shown in Fig. 1. The radius of the cylinder was 1.5 cm and the height was 3.0 cm. Probes A and B with the same radius of 0.2 cm were placed at the same height of Z = 1.5 cm and have an edge-to-edge distance (EED) of 0.2 cm. The cylinder was supposed to be filled with 1% intralipid solution with homogenous optical properties. The concentrations of probes A and B were set to 2 mg/ml and 2.5 mg/ml, respectively. The total quantum yields of two probes were set same for simplification. The relative fractions of probes A and B at four different spectral bands which have the center wavelengths of 580 nm, 620nm, 660 nm and 700 nm with about 40 nm full width at half maximum (FWHM) were set as Table 1.  The simulation model was rotated for 360°, with 12 projections acquired from the cylinder surface with an angular increment of 30°. Then, zero-mean white Gaussian noise was added to simulate the real noisy boundary measurements with signal-to-noise ratio (SNR) set as 35 dB.

Phantom experiment
A custom-made multispectral CB-XLCT system developed by our laboratory was adopted in the phantom experiments for validation. As shown in Fig. 2(a), the CB-XLCT system consists of a micro-focus X-ray source (Oxford Instrument, U.K.) with the maximal power of 80 W, a highly sensitive electron-multiplying CCD (EMCCD) camera (iXon DU-897, Andor, U.K.), which is coupled with a 50 mm f/1.8D lens (Nikon, Melville, N.Y.) for luminescence data acquisition, and a CMOS X-ray detector (2923, Dexela, U.K.) for the transmitted X-ray collection. EMCCD is placed in a 4 mm thickness lead box for protecting it from X-ray irradiation. The optical filter (Rayan Technology, China) in front of the EMCCD is used to choose the pass band in multispectral CB-XLCT. As shown in Fig. 2(b), a transparent glass cylinder (outer diameter 3.0 cm) embedded two glass tubes (outer diameter 0.4 cm) with an EED about 0.22 cm was used for the phantom experiments. The glass cylinder was filled with 1% intralipid solution diluted by water to mimic biological tissue. As shown in Fig. 2(c), the glass tubes filled with different nanoparticle powders Y 2 O 3 :Eu 3+ and Gd 2 O 2 S:Tb 3+ were used as multiple probes. Y 2 O 3 :Eu 3+ and Gd 2 O 2 S:Tb 3+ have several obvious emission peaks with the maximum at the wavelengths of 610 nm and 545 nm respectively [15,28]. Two band-pass optical filters which have the center wavelengths of 580 nm and 620 nm with about 40 nm FWHM were used for luminescence data collection in multispectral CB-XLCT. The region for reconstruction is about 3 cm height and the corresponding optical parameters were set at each wavelength [29].
To keep the total number of projections same as that in conventional CB-XLCT with no optical filters, the number of projections in multispectral CB-XLCT was reduced to 12 at each spectral band for two spectral bands are studied. The luminescence data were acquired at 24 projections with 15° step in the conventional CB-XLCT with no optical filters, and at 12 projections with 30° step in the multispectral CB-XLCT. During the luminescence data collection, the X-ray tube voltage and current were set to 40 kVp and 1 mA, respectively. The integration time, EM gain and binning of EMCCD were set to 1 s, 260 and 1 × 1, respectively.

Evaluation
To evaluate the spatial similarity between the actual regions and the computed regions, the Dice coefficient is calculated as [30]: where |·| denotes the number of voxels. X describes actual regions and Y describes the computed regions which are obtained by M-SI or M-PCA methods. The bigger value of Dice coefficient stands for the higher spatial similarity and the max value is 1. The location error (LE) is used to evaluate the localization accuracy of the proposed methods. It is given by calculating the Euclidean distance between centroids of the actual and computed regions: where a p and c p denote the centroids of the actual and computed regions which are obtained by M-SI or M-PCA methods, respectively.
In order to further quantitatively analyze the performance of the proposed methods in resolving the two targets, a relative coefficient is defined as [31]: values corresponding to two targets. According to the definition, the value of R denotes the spatial resolution, where R = 1 represents the highest spatial resolution and R = 0 represents the lowest spatial resolution. The reconstruction results at different spectral bands are illustrated in Fig. 3. As shown in Fig.  3, two probes A and B are mixed together, and it is difficult to resolve them from these reconstructed tomographic images at each single spectral band. The height of the displayed reconstruction slice is 1.5 cm. The actual boundaries of two probes are located by two small green circles.

The results of M-SI method
The results of traditional M-SI method are illustrated in Fig. 4. As shown in Fig. 4, the distributions of probes A and B can be recovered by M-SI method, respectively. The M-SI method is based on the first two spectral bands with center wavelengths of 580 nm and 620 nm.

Application of M-PCA on different combinations of multi-spectral data
Different combinations of spectral bands were used to test the performance of M-PCA method. As shown in Table 2, four spectral bands with center wavelengths of 580 nm, 620 nm, 660 nm and 700 nm were used as case 1, three spectral bands with center wavelengths of 580 nm, 620 nm and 660 nm were used as case 2, and two spectral bands with center wavelengths of 580 nm and 620 nm were used as case 3.  PCs obtained by M-PCA method of case 1, case 2 and case 3 are shown in Fig. 5. For each case, the images in the first row are the positive PCs, and the images in the second row are the negative PCs. In case 1, as shown in Fig. 5(a), the probes A and B are resolved by N2 image and P2 image, which were the negative and positive elements of PC2, respectively. Higher-order components, PC3 and PC4, contained less structure information of targets and their intensities were much smaller than that of PC2. In case 2, as shown in Fig. 5(b), the probes A and B are also well resolved by N2 image and P2 image, respectively. Similarly, the higher-order component, PC3, kept less information of target signals. As shown in Fig. 5(c), the probes A and B are resolved by P2 image and N2 image in case 3, respectively.
The transformation course from centered original data X 0 to the obtained P is described in Eq. (9). While transforming X 0 to P, the total variance keeps the same value which means that PCA does not change the total energy. In this work, the cumulated percentages of variances were calculated [32]. The percentages of variance (eigenvalue) in each PC and cumulated percentages of the first several PCs are displayed in Table 3. As shown in Table 3, the first two PCs account for most of the total variance for all cases.  As shown in Fig. 6, the images of two probes are merged in the same image in order to evaluate the spatial resolution. Figure 6(a) is the merged results of probes A and B recovered by M-SI method illustrated in Figs. 4(a) and 4(b). Figure 6(b) is the merged results of PC2 images by M-PCA method in case 3 illustrated in Fig. 5(c). All images are normalized to their maximal intensive value before mergence for comparison. The profiles across the two probes along the red dash lines (Y-axis) in Figs. 6(a) and 6(b) are shown in Fig. 6(c) in order to better clarify the spatial resolution. The FWHM of the profiles obtained with M-SI is larger than that of M-PCA and the valley value between the two peaks obtained with M-SI is much larger than that of M-PCA. To quantitatively analyze the spatial resolution, the index R defined in Eq. (13) for each profile in Fig. 6(c) is listed in Table 4. As shown in Table 4, R value is 0.9750 by M-PCA method in case 3 and 0.0639 by M-SI method, which illustrated that the M-PCA method can obtain much higher spatial resolution than M-SI method.
The performances of M-SI and M-PCA methods for all cases are further evaluated quantitatively with Dice, LE, and R indexes, as shown in Table 4. The merged PC2 images obtained by M-PCA and the merged image by M-SI were used to compute these indexes and LE was the mean value of LEs of both probes. Compared with M-SI, higher Dice values indicated that M-PCA reserved the target shape better, lower average LE indicated that the positions of the two probes determined by M-PCA method was closer to their actual positions, and higher R indicated that M-PCA obtained much higher spatial resolution. In addition, the proposed M-PCA method exhibited similar performances for all cases, especially those for case 1 and case 2. It indicates that the proposed M-PCA method could resolve two targets with only two spectral bands, which was validated qualitatively by Fig. 6 and quantitatively by Table 4.

Reconstructions of C-NF with different projections
The reconstruction images of C-NF method with different projections are shown in Fig. 7. It is difficult to resolve the adjacent multiple probes from these reconstructed images even with 24 projections in C-NF method. The qualities of the reconstruction images become worse with the number of projections decreasing.

The results of the M-SI method
In order to get the spectrum information priors for M-SI method, a spectral analysis were performed in advance. The emitted photons from nanoparticles Y 2 O 3 :Eu 3+ and Gd 2 O 2 S:Tb 3+ at the wavelengths from 400 nm to 800 nm with an interval of 1 nm were measured by spectrometer. Then ( ) k η λ of the pass bands of optical filters and thus the relative fractions ( ) k ω λ of the band pass optical filters with center wavelengths of 580 nm and 620 nm can be calculated and are shown in Table 5. The reconstruction results at spectral bands with the center wavelengths of 580 nm and 620 nm are shown in Figs. 8(a) and 8(b), respectively. It is hard to resolve the multiple probes from the reconstruction results. The reconstruction values are normalized by the maximum value of both spectral bands. As shown in Figs. 8(c) and 8(d), the distributions of nanoparticles Y 2 O 3 :Eu 3+ and Gd 2 O 2 S:Tb 3+ can be obtained by M-SI method, respectively.  The PCs obtained by M-PCA method are shown in Fig. 9. The images in the first row and in the second row are the positive and negative PCs. The nanoparticle powder Y 2 O 3 :Eu 3+ and Gd 2 O 2 S:Tb 3+ in the glass tubes on behalf of multiple probes are illustrated by the positive and negative elements in PC2 respectively, as shown in P2 image in Fig. 9(b) and N2 image in Fig. 9(d).

Evaluation
The performance of C-NF, M-SI and M-PCA methods are compared in Fig. 10 Table 6, M-PCA method has the higher Dice coefficient, smaller LE value and much larger R coefficient than M-SI method. These indexes quantitatively confirm M-PCA has better performance on resolving multiple probes than M-SI as shown in Fig. 10.

Discussion and conclusion
CB-XLCT has got development for its higher utilization of X-ray and shorter scanning time compared to the narrow-beam XLCT, but it suffers from the low-spatial resolution due to the ill-posed optical reconstruction. So it is hard to resolve the multiple reconstruction targets especially when they are adjacent in conventional CB-XLCT with no optical filters. In this work, a spectral-resolved CB-XLCT method which combined multispectral CB-XLCT with PCA was proposed for multiple probe imaging, and the performance was evaluated by digital simulation and phantom experiment further. Two probes with EED of 0.2 cm in digital simulation and two tubes with EED about 0.22 cm in phantom experiments were used as the adjacent multiple targets. As shown in the reconstruction images, the two targets cannot be separated through C-NF (Fig. 7) or CB-XLCT with single spectral band (Figs. 3, 8(a) and 8(b)). Meanwhile it is difficult to recognize the structure of the targets. But two targets can be resolved both by M-SI (Figs. 4, 8(c) and 8(d)) and M-PCA (Figs. 5, 9(b) and 9(d)) methods. Compared to M-SI method, the proposed M-PCA method has better performance on resolving adjacent multiple targets, which can be validated by the merged images (Figs. 6 and 10) qualitatively and the evaluation indexes (Tables 4 and 6) quantitatively. In addition, the performances of the proposed M-PCA method are tested with different combinations of spectral bands in digital simulations. The similar performance in case 3 to that in other cases indicated that two targets can be resolved by the proposed M-PCA method with only two spectral bands which was further validated by the phantom experiments in this work.
According to the principle of PCA, the recovered signals generally appeared in the first several components which has been testified in the previous studies [7,8,12,22,23]. How many principal components should be chosen can be determined to many rules [33]. In this work, the cumulated percentages of variances were calculated to choose the proper principal components which have been validated in optical imaging modalities [32]. For the resolving of two targets, most information of target signals is kept in the second PC, as shown in previous studies [7,8,12]. In this study, two targets with different behaviors were successfully resolved in PC2, too. But for the cases of low signal to noise ratio (SNR), two targets may appear in the higher-order PCs, such as in PC3 [23], so more PCs may needs to be considered for those special cases. In addition, more PCs should be considered if there are three or more targets, and combination of several components may be work in the cases for which single component is incapable.
Clearly with the multispectral detection, just as multi-voltage excitation, the proposed method would increase the imaging time and thus the radiation dose compared to the conventional CB-XLCT. In the phantom experiments, the projections at each spectral band were reduced to 12 in order to keep the total projection numbers same as that in the conventional CB-XLCT with no optical filters. Another purpose for this is to eliminate the effects of the increasing amount of the detection data on the performance of the proposed method. In addition, as described in our previous work in which two voltage excitations were used, radiation dose measured for 48 projections are still acceptable for in vivo imaging [8]. Several schemes could be taken to decrease the radiation dose further. Firstly, the total number of projections should be as small as possible. On one hand, fewer spectral bands should be used. In this work, experimental results indicated that two adjacent probes could be accurately resolved by the proposed method with only two spectral bands. On the other hand, the projections at each spectral band in multispectral CB-XLCT should be limited. But with the decreased projections, the reconstruction quality became worse just as shown in Fig. 7. The reconstruction method with limited projections can be considered [6]. Secondly, using the mirror design to collect projections of different angles simultaneously in the multispectral CB-XLCT as that in other optical imaging method may increase the time efficiency and decrease the X-ray dose [34].
Compared to other optical imaging modalities such as fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT), XLCT is newly developed and its advance depends greatly on the development of X-ray excitable luminescent nanoparticles (XLNP). Currently, the synthesis of water-soluble and bio-compatible XLNPs is still a challenge for most groups working on XLCT, resulting in few in vivo XLCT studies reported so far. Our group has put great effort on the synthesis of XLNPs and successfully prepared a waterdispersible β-NaGdF 4 :15% Eu 3+ nanoparticles in the end of last year [35]. However, in the present work, different types of soluble XLNP were needed to exhibit different excitation spectra, which have not been available yet. Therefore, different kinds of nanoparticle powders instead of solutions were used to test the feasibility of the proposed method. With the development of water-dispersible XLNPs, it is expected that XLCT imaging would advance greatly by phantom experiments, in vitro and in vivo small animal studies with XLNP solutions.
In conclusion, the proposed method which combined multispectral CB-XLCT with PCA provides an attractive approach for imaging multiple probes according to the results of digital simulation and phantom experiment. For multiple probes, the emission luminescence shows different variation trends with spectrum. Since PCA can catch emission behaviors of different varying trends, the M-PCA method can be applied to resolve targets filled with different nanoparticles. It would not work for targets containing the same nanoparticle with different concentrations because they have the same varying trend. If the linear assumption is not satisfied or the emission spectrum changes with concentrations, more research is needed. The independent component analysis could be considered further in multispectral CB-XLCT to get the concentration information.

Disclosures
The authors declare that there are no conflicts of interest related to this article.