Limited view cone-beam x-ray luminescence tomography based on depth compensation and group sparsity prior

Abstract. Significance: As a promising hybrid imaging technique with x-ray excitable nanophosphors, cone-beam x-ray luminescence computed tomography (CB-XLCT) has been proposed for in-depth biological imaging applications. In situations in which the full rotation of the imaging object (or x-ray source) is inapplicable, the x-ray excitation is limited by geometry, or a lower x-ray excitation dose is mandatory, limited view CB-XLCT reconstruction would be essential. However, this will result in severe ill-posedness and poor image quality. Aim: The aim is to develop a limited view CB-XLCT imaging strategy to reduce the scanning span and a corresponding reconstruction method to achieve robust imaging performance. Approach: In this study, a group sparsity-based reconstruction method is proposed with the consideration that nanophosphors usually cluster in certain regions, such as tumors or major organs such as the liver. In addition, depth compensation (DC) is adopted to avoid the depth inconsistency caused by a limited view strategy. Results: Experiments using numerical simulations and physical phantoms with different edge-to-edge distances were carried out to illustrate the validity of the proposed method. The reconstruction results showed that the proposed method outperforms conventional methods in terms of localization accuracy, target shape, image contrast, and spatial resolution with two perpendicular projections. Conclusions: A limited view CB-XLCT imaging strategy with two perpendicular projections and a reconstruction method based on DC and group sparsity, which is essential for fast CB-XLCT imaging and for some practical imaging applications, such as imaging-guided surgery, is proposed.

1 Introduction sensitivity due to the avoidance of background optical signals and autofluorescence. With the above advantages, great efforts have been devoted to XLCT imaging and several types of XLCT systems have been proposed according to the x-ray beam shapes. Narrow-beam 2,3 and pencilbeam XLCT 8,9 can achieve high spatial resolution, but the long imaging time hinders their application to fast biomedical applications. To reduce imaging time, cone-beam XLCT (CB-XLCT) systems [10][11][12][13] have been proposed for fast imaging to skip the translation step essential in narrowand pencil-beam XLCT systems.
However, due to high light scattering and low absorption properties in biological tissues, the reconstruction of CB-XLCT is an ill-posed problem. To improve the imaging performance, much prior information has been introduced for CB-XLCT reconstruction, such as permissible region, 14 CT images, 15 and sparsity. 16 Sparsity is the most frequently used prior information for CB-XLCT reconstruction since nanophosphors are usually sparsely distributed inside the body. With this information, in our previous studies, sparse view-based reconstruction methods have been proposed by Gao et al. 16 and Liu et al., 17 where two adjacent targets can be recovered accurately with only four projections evenly distributed in a 360-deg span. Considering this, in some applications, such as intraoperative breast cancer lumpectomy, 18 only limited view projections could be acquired due to geometry limitation and/or low-dose excitation. It is essential to develop a reconstruction method for the limited view strategy to achieve robust imaging performance in these situations.
Though Liu et al. proposed single-view CB-XLCT reconstruction methods based on sparsity 19,20 or group sparsity, 21 it either works for single target recovery or relies on CT images to obtain group prior, whereas the reconstruction is related to the group LASSO model and imaging performance may be affected by inaccurate group information. In optical imaging, nanophosphors often cluster in some areas, such as tumor regions where luminescent signals are simultaneously nonzero, and rarely distribute in other areas, such as normal tissues where luminescent signals are simultaneously near zero. Thus, beyond sparsity, nanophosphors can also be characterized by group sparsity. 22 With this prior information, instead of prior CT information, it is possible to improve the image quality of limited view CB-XLCT for multiple targets.
In addition, due to the nonlinear depth sensitivity of measurements detected, the optical intensity of targets far from the detector would be less than those near the detector, 23,24 resulting in low reconstructed values for far targets. Therefore, for optical imaging, especially for limited view acquisition, depth compensation (DC) is essential for avoiding depth inconsistency and obtaining high quality reconstructed images.
In this paper, a DC method is first adopted to level off the differences in detection sensitivities by incorporating two weight matrices into the optimization function. 23 To alleviate the ill-posedness, group sparsity is then introduced as the prior information. Unlike Liu et al. who used the group LASSO model with CT information, here, we add the fused LASSO (FL) penalty to the objective function, which was solved efficiently with a split Bregman iterative method. To evaluate the performance of the proposed method with the DC-FL penalty for CB-XLCT, both numerical simulation and physical phantom experiments were performed and the proposed algorithm was compared with the adaptive Tikhonov (Adaptik), 25 fast iterative shrinkage-thresholding algorithm (FISTA), 26 and our previously proposed sparse view-based T-FISTA method 16 in terms of location error, imaging contrast, shape, and spatial resolution.
The structure of the paper is organized as follows. In Sec. 2, the methodology is presented. First, the forward model for CB-XLCT is briefly described. Second, the DC matrix and FL penalty are formulated. Finally, the split Bregman FL method is described to solve the objective function. In Sec. 3, both numerical simulation and physical phantom experiments are performed to illustrate the performance of the proposed algorithm. In Sec. 4, the results from numerical simulations and physical phantom experiments are presented. In Sec. 5, some discussion and a conclusion are given.

Forward Problem
In the CB-XLCT system, x-rays penetrate the imaging object based on Lambert-Beers' law, and the x-ray intensity (W cm −3 ) distribution XðrÞ in the imaging object is expressed as follows: 10 where Xðr 0 Þ is the intensity of x-ray (W cm −3 ) at the initial position r 0 , and μ t ðτÞ is the x-ray attenuation coefficient (cm −1 ) obtained from x-ray transmission data. When irradiated by x-rays, nanophosphors distributed in the imaging object can emit visible or near-infrared (NIR) light, which is expressed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 6 5 1 SðrÞ ¼ εXðrÞnðrÞ; where SðrÞ is the source energy density (W cm −3 ), ε is the light yield defined as the quantum yield per unit nanophosphor concentration (mg mL −1 ), and nðrÞ is the nanophosphor concentration in position r. The light transportation in scattering media is modeled by the radiative transfer equation (RTE), but solving the RTE directly is extremely difficult. Considering the highly scattering and weakly absorbing properties of biological tissues in the visible or NIR spectral window, the RTE model is simplified to the following diffusion equation model: 7 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 5 4 8 −∇½DðrÞ∇ΦðrÞ þ μ a ðrÞΦðrÞ ¼ XðrÞ ðr ∈ ΩÞ ΦðrÞ þ 2αDðrÞ½ν∇ΦðrÞ ¼ 0 ðr ∈ ∂ΩÞ ; (3) where DðrÞ is the diffusion coefficient that is calculated by DðrÞ ¼ ½3½μ a ðrÞ þ μ 0 s ðrÞ −1 and μ a and μ 0 s are the absorption and reduced scattering coefficients of the tissue, respectively. ΦðrÞ is the photon fluence at position r, Ω is the domain of the imaged object, ∂Ω denotes the boundary of Ω, α is a factor describing the optical reflective index mismatch, and ν is the outward unit normal vector on ∂Ω.
With the finite-element method, Eq. (3) is discretized into a matrix equation, which builds a linear relationship between the nanophosphor concentration N and photon measurements on the object surface Φ meas : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 4 0 7 where A is a weight matrix used to map the unknown nanophosphor distribution to known measurements.

Limited View Reconstruction with DC-FL for CB-XLCT
It is known that the detection sensitivity in optical imaging decreases nonlinearly with increased depth. 24 This makes CB-XLCT measurements hypersensitive to targets near the detector, which typically happens in limited view imaging of multiple targets. The ill-posed weight matrix will lead to reconstructed biases toward the superficial targets. To level off this bias, DC is first implemented by incorporating two weight matrices into the objective function. Then, to take advantage of the priors in CB-XLCT, a FL penalty is added for reconstruction based on group sparsity priors. The final objective function is solved with the split Bregman method considering its computation efficiency. Figure 1 gives the flowchart of the proposed algorithm.

Depth compensation for the inverse problem
Generally, the weight matrix A is ill-posed, and it is impractical to solve N directly. To obtain a unique and stable solution, Eq. (4) is solved by the following minimization problem with the L p regularization term: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 1 4 6 min n fkAN − Φ meas k 2 þ λkNk P g; where λ is the regularization parameter and kNk p is the L p -norm of N.
For effective DC, two weight matrices, i.e., data weight matrix W d and model weight matrix W m are introduced, 23 and Eq. (5) are further converted to the following equations: In detail, W d provides constraints or a priori information for the solution, which is constructed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 6 ; 4 0 1 whereÑ is a rough approximation to the true solution N and q controls the compromise between a close fit to the data and the stability of the solution, which was set according to Ref. 24. In this paper, the rough approximation is obtained by the T-FISTA method. W m is used to level off the differences in detection sensitivities and is given by where β j is a normalization factor that is inversely proportional to the absolute largest difference between the elements within each column.

Fused LASSO penalty for the inverse problem
In this paper, the group sparsity of the luminescent signals, which means that the luminescent signals should be roughly piecewise constant, is introduced as a priori information. As a result, the FL penalty is incorporated into Eq. (6) to make use of the group sparsity: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 6 ; 1 0 6 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 6 ; 7 2 3 The regularization term with parameter λ 1 encourages the sparsity of the reconstructed signals, while the other term with parameter λ 2 shrinks the differences between neighboring signals.

Split Bregman method for the fused LASSO problem
Generally, solving Eq. (9) is computationally challenging due to the nondifferentiability of the objective function. Here, in this paper, we adopted a method based on the split Bregman iteration for solving the general FL problem. 27 The Bregman iteration technique is based on the Bregman distance that generalizes the concept of metric associating a distance to a convex function not necessarily differentiable, and it addresses the objective problem by analyzing it into several functions and minimizing them separately in an efficient simple way. We introduce two auxiliary variables, a and b, to represent Y and RY, and u and v are dual variables corresponding to the linear constraints Y ¼ a and RY ¼ b, respectively. The split Bregman iteration can be implemented as follows: Repeat:

Convergence
In this paper, the regularization parameters λ 1 and λ 2 were set to 0.1 and 1, respectively. The iteration numbers were set to 150 in simulations and 300 in phantom experiments according to the results, respectively. More detailed information on this method can be found in Ref. 28. 3 Experimental Setup

Numerical Simulation Setup
Numerical simulations were implemented with a cylinder phantom with two luminescent targets embedded to evaluate the performance of the proposed method. The geometry configuration of the phantom is shown in Fig. 2. It was a cylinder placed on a rotating stage, with the rotational axis defined as the Z axis and the bottom plane set as Z ¼ 0 cm. The diameter and height of the cylindrical phantom were 3.0 and 2.5 cm, respectively. To approximate high scattering media, the phantom was filled with 1% intralipid, where the absorption coefficient and reduced scattering coefficient were set to 0.02 and 10 cm −1 , respectively. Two luminescent targets with nanophosphors of Y 2 O 3 ∶Eu 3þ were placed in the phantom at a depth of 1.0 cm. The targets were 0.4 cm in diameter and 0.5 cm in height. The edge-to-edge distance (EED) of the two targets was set to 0.30 cm in case 1 and 0.10 cm in case 2. The nanophosphors were excited by CB x-rays with a tube voltage of 40 kV. The initial phantom position setup is shown in Fig. 3. Luminescent images were collected at every 90 deg for a full span of 360 deg. A total of four projections were collected, some of which were used for reconstruction. The detectors located on finite-element nodes of the boundary were inside 160 deg. To make the simulations more realistic, white Gaussian noise was added to generate noisy boundary measurements with a signal-to-noise ratio set to 20 dB. In this paper, the analytical simulation was conducted with COMSOL Multiphysics 3.3 (COMSOL, Inc., Burlington, Massachusetts).

Phantom Experiments Setup
A physical phantom similar to the simulations was imaged based on the custom-made CB-XLCT system developed by our lab. The system consists of a microfocus CB x-ray source that is used to excite the phantom, a highly sensitive electron-multiplying charge-coupled device (EMCCD) adopted for luminescent projections collection, and an x-ray flat panel detector to collect transmitted x-ray signals, as described in Fig. 3. Figure 4 illustrates the setup of the physical phantom. It was a transparent glass cylinder with a diameter of 3.0 cm and a height of 7.0 cm, filled with 1% intralipid and water (with an absorption coefficient of 0.02 cm −1 and a reduced scattering coefficient of 10 cm −1 ). Two glass tubes of 0.4 cm diameter that contained nanophosphors of Y 2 O 3 ∶Eu 3þ were implanted in the phantom. The concentration of Y 2 O 3 ∶Eu 3þ was 0.1 g mL −1 . Similar to simulation studies, two cases of phantom experiments with different EEDs (0.50 cm in

Quantitative Evaluation
To quantitatively evaluate the performance of the proposed reconstruction method, in this paper, three state-of-the-art methods were implemented for comparison, including the adaptive Tikhonov regularization (Adaptik, L-2 norm), FISTA (L-1 norm), and our previously proposed revised T-FISTA method. The regularization parameters were set according to Ref. 24, respectively. The iteration numbers were chosen empirically according to the results and to ensure the convergence of the calculation.
Images reconstructed by different methods were compared in terms of location error, imaging contrast, shape, and spatial resolution. For quantitative evaluation, the position error (PE), Dice similarity coefficient (DICE), contrast-to-noise ratio (CNR), and spatial resolution index (SPI) were calculated, respectively. 16 PE is defined as the Euclidean distance between the actual and reconstructed luminescent positions: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 6 ; 1 5 9 where p r and p t denote the centers of the reconstructed and true targets, respectively. DICE evaluates the similarity between the actual and reconstructed luminescent areas: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 6 ; 1 0 2 DICE ¼ 2jROI r ∩ ROI t j jROI r j þ jROI t j ; where ROI r and ROI t denote the reconstructed and true luminescent areas, respectively. The closer the reconstruction is to the true target, the closer the DICE is to 1. Otherwise, the DICE is closer to 0. CNR is used for quantitative evaluation of image contrast: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 6 ; 6 8 7 where the subscripts ROI and BK denote the target and background regions, respectively, mean and var denote the mean value and variance, respectively, and w is the weighting factor calculated by the relative volumes of the ROI. SPI is a spatial resolution quantitative index to analyze the performance of the algorithms in resolving two targets: where n l denotes the value of the profile along a given line that connects the two centers on the reconstructed cross-section. n l max , n l min , and n l valley are the maximal, minimal, and valley value between the two peak values, respectively. The more clearly the two targets are separated, the closer the SPI is to 1. Figure 5 shows the reconstruction results of case 1 with two projections in numerical simulations. The height of the cross-section is 1.3 cm. The red circle in each image depicts the phantom boundary, and the two yellow circles depict the true positions of the two targets. Two groups of reconstructions were conducted. In the first group, the projection pair of 0 deg to 180 deg was adopted for reconstruction, and the results are shown in the first row. It can be seen that both the traditional L2-norm based method Adaptik and the L1-norm based method FISTA could not resolve the two targets and introduce more noise in the images. Though the T-FISTA method could resolve the two targets clearly, some noise existed in the image when compared with that obtained by the proposed method. In the second group, a reduced imaging span was adopted and the projection pair of 0 deg to 90 deg was used for reconstruction. The results are shown in the second row of Fig. 5. In this group, the other three methods could not resolve two targets effectively, while the proposed algorithm could do it with good image quality.

Numerical Simulations
To demonstrate whether the initial position of the 90-deg span would affect the imaging performance, we performed another numerical simulation, where three other 90-deg span projection pairs, 30 deg to 120 deg, 60 deg to 150 deg, and 90 deg to 180 deg, were used, respectively. The results are shown in Fig. 6. It can be seen that, no matter where the projection started, the quality of reconstructed images was quite similar. This indicates that the initial position of the limited view span has little effect on the reconstruction results. Figure 7 shows the reconstruction results of simulation case 2, where two projection pairs of 0 deg to 90 deg were used. The upper row in Fig. 7 gives 2D slices while the lower row shows corresponding 3D visualization results. In 3D results, values of <10% of the maximum value were neglected. The cylinders depict the phantom, and the red objects represent the luminescent targets. We can see that, as the two targets get closer, the tomographic images obtained with the first three methods [Figs. 5(a)-5(c)] become much worse than those in case 1, and none of these methods could separate the two targets clearly. By contrast, the proposed DC-FL method could achieve high-quality imaging with accurate position and less noise, which can be further demonstrated by the 3D renderings.
To further evaluate the performance of the proposed method, quantitative analysis was performed, as shown in Table 1, which corresponds to the upper row in Fig. 7. The proposed method  achieves the least PE, which means the reconstructed positions of the two targets were most accurate compared with those obtained by the other three methods. In addition, the proposed method yields the highest DICE, CNR, and SPI, indicating that the targets were recovered with the best shape, least noise, and best spatial resolution. Figure 8 shows the cross-sectional and 3D results of the phantom experiment. The first and third rows show the tomographic fused XLCT/CT images corresponding to the CT slice indicated by  To demonstrate whether the initial position would affect the imaging performance of the phantom experiments, three other 90-deg span projection pairs, including 30 deg to 120 deg, 60 deg to 150 deg and 90 deg to 180 deg were tested, respectively. As shown in Fig. 9, images reconstructed from different projection pairs are quite similar, indicating that the initial position of the limited view span has little effect on the reconstruction results.

Phantom Experiments
To further validate the proposed algorithm, a quantitative evaluation of phantom experiment 2 (third row in Fig. 8) was carried out, as presented in Table 2. Similar to results in simulations, this demonstrates that, compared with the other methods, the proposed method yields the smallest LE and the highest CNR, DICE, and SPI, indicating that the targets were recovered with the least relative error.

Discussion and Conclusion
In this study, we established a group sparsity-based limited view CB-XLCT reconstruction method. To reduce the depth-dependent value inconsistency caused by the limited view strategy, DC was first adopted. Then, the group sparsity was introduced as prior information in consideration of nanophosphors usually being clustered in groups. Both numerical simulations and phantom experiments validate the performance of the proposed method.  The reconstruction results of numerical simulations and phantom experiments show that, with the proposed method, we can resolve the distribution of the nanophosphors even when the targets are close to each other, using only two projections at 0 deg and 90 deg (Figs. [6][7][8]. Quantitative analysis (Tables 1 and 2) together with the 3D renderings indicates that, compared with Adaptik, FISTA, and T-FISTA methods, the proposed DC-FL method achieved the highest location accuracy, contrast, and resolution and the most consistent shape. All of these results demonstrate the potential of the proposed method in improving imaging performance and reducing the imaging time and radiation dose of limited view CB-XLCT. Though DC-FL is designed for CB-XLCT in this paper, it can be extended to other XLCT systems, such as narrow-beam XLCT or fan-beam XLCT systems. Further, with appropriate forward models, DC-FL can be applied to other optical reconstructions, such as FMT and BLT.
Our results indicate that no matter which projection is used as the initial view (0 deg, 30 deg, 60 deg, or 90 deg), the imaging performance is similar (Fig. 6), which suggests that 90-deg span could be a robust imaging strategy for limited view imaging. We also find that a smaller span (single view in simulations or 60 deg in phantom experiments) may also achieve comparable performance. However, this depends on the position of the initial view, which may lead to unstable results. In addition, the imaging cost has been greatly reduced, i.e., from 94 s (with four projections evenly collected in a 360-deg span) to 17 s (with two perpendicular projections collected in a 90-deg span) in our experiments.
Although the proposed method achieved better performance, some important issues need to be further addressed. In this work, all of the hyperparameters including regularization parameters, iteration numbers, and initial values were empirically optimized based on value ranges suggested by references. The automatic selection of these parameters would substantially benefit CB-XLCT. In addition, the rough approximation of the solution used for constructing W d is critical for the imaging performance of the proposed method. The more accurate the approximation is, the better the result that can be achieved. Furthermore, effective optimization methods used to solve the FL problem, such as the primal-dual Newton conjugate gradient method, 29 may further improve the limited view imaging. In future work, more effort will be devoted to resolving more targets with smaller EED by the proposed method.

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