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

Peng Gao; Junyan Rong; Tianshuai Liu; Wenli Zhang; Bin Lan; Xiaoping Ouyang; Hongbing Lu

doi:10.1117/1.JBO.25.1.016004

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

Peng Gao, Junyan Rong, Tianshuai Liu, Wenli Zhang, Bin Lan, Xiaoping Ouyang, Hongbing Lu

Author Affiliations +

Journal of Biomedical Optics, Vol. 25, Issue 1, 016004 (January 2020). https://doi.org/10.1117/1.JBO.25.1.016004

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

X-ray luminescence computed tomography (XLCT) has attracted much attention as a CT/optical dual-mode imaging technique.¹^–³ It utilizes x-ray excitable nanophosphors as biomarkers to produce visible or near-infrared (NIR) light upon x-ray irritation and naturally combines x-ray structure imaging of high resolution with optical molecular imaging of high sensitivity and specificity. Compared with traditional optical molecular imaging modalities, such as fluorescence molecular tomography (FMT)⁴^,⁵ and bioluminescence tomography (BLT),⁶^,⁷ XLCT can detect targets deeply seated in the imaging object with the use of x-rays and provides higher 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²^,³ and pencil-beam XLCT⁸^,⁹ 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¹⁰^–¹³ have been proposed for fast imaging to skip the translation step essential in narrow- and 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,¹⁴ CT images,¹⁵ and sparsity.¹⁶ 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.¹⁶ and Liu et al.,¹⁷ 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,¹⁸ 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¹⁹^,²⁰ or group sparsity,²¹ 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.²² 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,²³^,²⁴ 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.²³ 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),²⁵ fast iterative shrinkage-thresholding algorithm (FISTA),²⁶ and our previously proposed sparse view-based T-FISTA method¹⁶ 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.

2. Methodology

2.1.

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:¹⁰

Eq. (1)

X (r) = X (r_{0}) \exp {- \int_{r_{0}}^{r} μ_{t} (τ) d τ},

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:

Eq. (2)

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

Eq. (3)

{\begin{cases} - \nabla [D (r) \nabla Φ (r)] + μ_{a} (r) Φ (r) = X (r) & (r \in Ω) \\ Φ (r) + 2 α D (r) [ν \nabla Φ (r)] = 0 & (r \in \partial Ω) \end{cases},

where

D (r)

is the diffusion coefficient that is calculated by

D (r) = {[3 [μ_{a} (r) + μ_{s}^{'} (r)]]}^{- 1}

and

μ_{a}

and

μ_{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,

\partial Ω

denotes the boundary of

Ω

,

α

is a factor describing the optical reflective index mismatch, and

ν

is the outward unit normal vector on

\partial Ω

.

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}$ :

Eq. (4)

AN = Φ_{meas},

where

A

is a weight matrix used to map the unknown nanophosphor distribution to known measurements.

2.2.

Limited View Reconstruction with DC-FL for CB-XLCT

It is known that the detection sensitivity in optical imaging decreases nonlinearly with increased depth.²⁴ 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.

Fig. 1

The flowchart of the proposed DC-FL algorithm.

2.2.1.

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:

Eq. (5)

\min_{n} {{‖ AN - Φ_{meas} ‖}^{2} + λ {‖ N ‖}_{P}},

where

λ

is the regularization parameter and

{‖ N ‖}_{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,²³ and Eq. (5) are further converted to the following equations:

Eq. (6)

{\begin{cases} \min_{y} {{‖ {AW}_{d} W_{m}^{- 1} Y - Φ_{meas} ‖}_{2}^{2} + λ {‖ Y ‖}_{P}^{2}} \\ Y = W_{m} N \end{cases} .

In detail, $W_{d}$ provides constraints or a priori information for the solution, which is constructed as

Eq. (7)

W_{d} = diag (\tilde{N}) = (\begin{matrix} \tilde{n_{1}^{q}} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & \tilde{n_{n}^{q}} \end{matrix}),

where

\tilde{N}

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

Eq. (8)

{\begin{cases} β_{j} = \frac{1}{| \max (A_{j}) - \min (A_{j}) |} \\ W_{m} = diag {β_{j} * {‖ A_{j} ‖}_{2}^{- 1}} = (\begin{matrix} β_{j} * {‖ A_{j} ‖}_{2}^{- 1} & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & β_{n} * {‖ A_{n} ‖}_{2}^{- 1} \end{matrix}), j = 1,2, \dots, N \end{cases},

where

β_{j}

is a normalization factor that is inversely proportional to the absolute largest difference between the elements within each column.

2.2.2.

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:

Eq. (9)

\min_{y} {{‖ {AW}_{d} W_{m}^{- 1} Y - Φ_{meas} ‖}_{2}^{2} + λ_{1} {‖ Y ‖}_{1}^{2} + λ_{2} {‖ R Y ‖}_{1}^{2}},

where

Eq. (10)

R_{i j} = {\begin{cases} - 1 & j = i, i = 1,2, \dots \\ 1 & j = i + 1, i = 1,2, \dots \\ 0 & otherwise \end{cases} .

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.

2.2.3.

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

Initialization: $y^{0}$ , $a^{0}$ , $b^{0}$ , $u^{0}$ , and $v^{0}$ .
Let $V (y) = {‖ {AW}_{d} W_{m}^{- 1} Y - Φ_{meas} ‖}_{2}^{2}$
Repeat:
1) $y^{m + 1} = \arg \min_{y} {V (y) + ⟨ u^{m}, y - a^{m} ⟩ + ⟨ v^{m}, R y - b^{m} ⟩ + \frac{μ_{1}}{2} {‖ y - a^{m} ‖}_{2}^{2} + \frac{μ_{2}}{2} {‖ R y - b^{m} ‖}_{2}^{2}}$
2) $a^{m + 1} = ζ_{u_{1}^{- 1} λ_{1}} (y^{m + 1} + μ_{1}^{- 1} u^{m})$
3) $b^{m + 1} = ζ_{u 2^{- 1} λ_{2}} (R y^{m + 1} + μ_{2}^{- 1} v^{m})$
4) $u^{m + 1} = u^{m} + δ_{1} (y^{m + 1} - a^{m + 1})$
5) $v^{m + 1} = v^{m} + δ_{2} (R y^{m + 1} - b^{m + 1})$ where $ζ$ a soft thresholding operator.
Until
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

3.1.

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

Fig. 2

Setup for the cylinder phantom. (a), (b) The view of the $X Y$ plane of the two cases. (c), (d) The view of the $X Z$ plane of the two cases. Two luminescent targets of $Y_{2} O_{3} : {Eu}^{3 +}$ were placed inside a cylinder phantom. The EED along the $X$ axis is 0.30 cm in case 1 (first column) and 0.10 cm in case 2 (second column).

Fig. 3

Schematic diagram of the hybrid CB-XLCT/XCT imaging system and the initial phantom position setup used in simulation and phantom experiments.

3.2.

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 case 1 and 0.23 cm in case 2) were performed. The detectors located on finite-element nodes of the boundary were inside 160 deg. The imaging details can be found in Ref. 16.

Fig. 4

Setup for experimental studies. (a), (b) Representative x-ray projections of the phantoms with $EED = 0.50 cm$ and $EED = 0.23 cm$ , respectively. Regions between the blue and green lines are used for analysis. (c), (d) XCT slice of the phantom, corresponding to the slice indicated by the red line in (a) and (b), respectively. The left and right tubes are named as tube 1 and tube 2, respectively. White scale bars in (a) and (c): 1 cm.

Figures 4(a) and 4(b) show the representative x-ray projections of the physical phantoms with different EEDs. Representative CT slices of the phantoms with different EEDs are shown in Figs. 4(c) and 4(d). Regions between the blue and green lines were selected for the study.

3.3.

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

PE is defined as the Euclidean distance between the actual and reconstructed luminescent positions:

Eq. (11)

PE = {‖ p_{r} - p_{t} ‖}_{2},

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:

Eq. (12)

DICE = \frac{2 | {ROI}_{r} \cap {ROI}_{t} |}{| {ROI}_{r} | + | {ROI}_{t} |},

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:

Eq. (13)

CNR = \frac{| {Mean}_{ROI} - {Mean}_{BK} |}{{(w {Var}_{ROI}^{2} + (1 - w) {Var}_{BK}^{2})}^{1 / 2}},

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:

Eq. (14)

SPI = \frac{n_{\max}^{l} - n_{valley}^{l}}{n_{\max}^{l} - n_{\min}^{l}},

where

n_{l}

denotes the value of the profile along a given line that connects the two centers on the reconstructed cross-section.

n_{\max}^{l}

,

n_{\min}^{l}

, and

n_{valley}^{l}

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.

4. Results

4.1.

Numerical Simulations

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.

Fig. 5

CB-XLCT reconstruction results of simulation case 1 ( $EED = 0.30 cm$ ). (a)–(h) The first to the fourth column are results obtained with the Adaptik, FISTA, T-FISTA, and DC-FL. The first and second row are the results reconstructed with two projections of 180-deg and 90-deg view spans, respectively. The yellow circles depict the true positions of the two targets, and the red circles depict the boundary of the phantom. All images are normalized to the maximal value.

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.

Fig. 6

The reconstruction results for illustrating the effect of the initial position of the 90-deg span. (a)–(c) Three different 90-deg span projection pairs: (a) 30 deg to 120 deg, (b) 60 deg to 150 deg, (c) and 90 deg to 180 deg.

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.

Fig. 7

(a)–(h) CB-XLCT reconstruction results of simulation case 2 ( $EED = 0.10 cm$ ). From left to right: results obtained by the Adaptik, FISTA, T-FISTA, and DC-FL, respectively. From top to bottom: 2D and 3D results reconstructed with two projection pairs of 0-deg to 90-deg view span, respectively. The cylinders depict the phantom, and the red objects represent the luminescent targets.

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.

Table 1

Quantitative analysis of simulation case 2 with two projection pairs of 0 deg to 90 deg.

	PE (mm)		DICE		CNR	SPI
	Tube 1	Tube 2	Tube 1	Tube 2	CNR	SPI
Adaptik	3.79	2.97	0.07	0.06	0.70	0.37
FISTA	3.81	3.05	0.05	0.05	1.40	0.30
T-FISTA	5.94	0.82	0.06	0.06	1.84	0.23
DC-FL	0.22	0.29	0.82	0.73	8.38	0.94

4.2.

Phantom Experiments

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 the red line in Fig. 4(a), while the second and fourth rows depict the XLCT 3D visualization results of case 1 ( $EED = 0.50 cm$ ) and case 2 ( $EED = 0.23 cm$ ), respectively. Reconstructions obtained by the Adaptik, FISTA, T-FISTA, and DC-FL are shown from the first to the fourth column. In both cases, two projection pairs of 0 deg to 90 deg were adopted for reconstruction.

Fig. 8

Results of the phantom experiment. The first and third row shows the tomographic fused XLCT/CT images, while the second and fourth row depicts the XLCT 3D visualization results of (a)–(h) case 1 ( $EED = 0.50 cm$ ) and (i)–(p) case 2 ( $EED = 0.23 cm$ ), respectively. Reconstructions obtained by the Adaptik, FISTA, T-FISTA, and DC-FL are shown from the first to the fourth column.

The phantom experimental results indicate that the limited view used in the reconstruction could only resolve one target [Figs. 8(a)–8(c) and 8(i)–8(k)] with the Adaptik, FISTA, and T-FISTA methods, which might be due to one target being closer to the detector than the other in the two projections. With the DC-based FL reconstruction strategy, both targets could be resolved accurately as expected [Fig. 8(a)], even when the two tubes were close with an EED of 0.23 cm [Fig. 8(l)]. The performance was further confirmed by the 3D results [Figs. 8(h) and 8(p)].

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.

Fig. 9

Results of the phantom experiment case 2 with different 90-deg span projection pairs. (a)–(c) Three different 90-deg span projection pairs: (a) 30 deg to 120 deg, (b) 60 deg to 150 deg, and (c) 90 deg to 180 deg.

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.

Table 2

Quantitative analysis of phantom results of case 2.

	PE (mm)		DICE		CNR	SPI
	Tube 1	Tube 2	Tube 1	Tube 2	CNR	SPI
Adaptik	5.65	1.33	0.01	0.23	1.68	0.24
FISTA	4.85	3.44	0.11	0.17	1.91	0.10
T-FISTA	3.52	3.45	0.28	0.26	2.78	0.02
DC-FL	0.93	0.64	0.48	0.68	5.00	0.98

5. 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–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,²⁹ 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.

Acknowledgments

This project was supported by National Key Research and Development Program of China (Grant Nos. 2017YFC0107400, 2017YFC0107401, and 2017YFC0107403), National Natural Science Foundation of China (Grant No. 81230035), Natural Science Foundation of Shaanxi Province (Grant No. 2016JQ1012), and Key Project supported by Military Science and Technology Foundation (Grant No. BWS14C030).

References

1.

C. M. Carpenter et al., “Hybrid x-ray/optical luminescence imaging: characterization of experimental conditions,” Med. Phys., 37 (8), 4011 –4018 (2010). https://doi.org/10.1118/1.3457332 MPHYA6 0094-2405 Google Scholar

2.

G. Pratx et al., “Tomographic molecular imaging of x-ray-excitable nanoparticles,” Opt. Lett., 35 (20), 3345 –3347 (2010). https://doi.org/10.1364/OL.35.003345 OPLEDP 0146-9592 Google Scholar

3.

G. Pratx et al., “X-ray luminescence computed tomography via selective excitation: a feasibility study,” IEEE Trans. Med. Imaging, 29 (12), 1992 –1999 (2010). https://doi.org/10.1109/TMI.2010.2055883 ITMID4 0278-0062 Google Scholar

4.

H. Pu et al., “Separating structures of different fluorophore concentrations by principal component analysis on multispectral excitation-resolved fluorescence tomography images,” Biomed. Opt. Express, 4 (10), 1829 –1845 (2013). https://doi.org/10.1364/BOE.4.001829 BOEICL 2156-7085 Google Scholar

5.

X. Cao et al., “Reconstruction for limited-projection fluorescence molecular tomography based on projected restarted conjugate gradient normal residual,” Opt. Lett., 36 (23), 4515 –4517 (2011). https://doi.org/10.1364/OL.36.004515 OPLEDP 0146-9592 Google Scholar

6.

Y. Gao et al., “Bioluminescence tomography based on Gaussian weighted Laplace prior regularization for in vivo morphological imaging of glioma,” IEEE Trans. Med. Imaging, 36 (11), 2343 –2354 (2017). https://doi.org/10.1109/TMI.2017.2737661 ITMID4 0278-0062 Google Scholar

7.

W. Cong et al., “Practical reconstruction method for bioluminescence tomography,” Opt. Express, 13 (18), 6756 –6771 (2005). https://doi.org/10.1364/OPEX.13.006756 OPEXFF 1094-4087 Google Scholar

8.

C. Li, A. Martinez-Davalos and S. R. Cherry, “Numerical simulation of x-ray luminescence optical tomography for small-animal imaging,” J Biomed. Opt., 19 (4), 046002 (2014). https://doi.org/10.1117/1.JBO.19.4.046002 Google Scholar

9.

M. C. Lun, W. Zhang and C. Li, “Sensitivity study of x-ray luminescence computed tomography,” Appl. Opt., 56 (11), 3010 –3019 (2017). https://doi.org/10.1364/AO.56.003010 APOPAI 0003-6935 Google Scholar

10.

D. Chen et al., “Cone beam x-ray luminescence computed tomography: a feasibility study,” Med. Phys., 40 (3), 031111 (2013). https://doi.org/10.1118/1.4790694 MPHYA6 0094-2405 Google Scholar

11.

X. Liu, Q. Liao and H. Wang, “Fast x-ray luminescence computed tomography imaging,” IEEE Trans. Biomed. Eng., 61 (6), 1621 –1627 (2014). https://doi.org/10.1109/TBME.2013.2294633 IEBEAX 0018-9294 Google Scholar

12.

P. Gao et al., “Resolving adjacent nanophosphors of different concentrations by excitation-based cone-beam x-ray luminescence tomography,” Biomed. Opt. Express, 8 (9), 3952 –3965 (2017). https://doi.org/10.1364/BOE.8.003952 BOEICL 2156-7085 Google Scholar

13.

H. Pu et al., “Spectral-resolved cone-beam x-ray luminescence computed tomography with principle component analysis,” Biomed. Opt. Express, 9 (6), 2844 –2858 (2018). https://doi.org/10.1364/BOE.9.002844 BOEICL 2156-7085 Google Scholar

14.

H. Yi et al., “A permissible region extraction based on a knowledge priori for x-ray luminescence computed tomography,” Multimedia Syst., 25 (2), 147 –154 (2019). https://doi.org/10.1007/s00530-017-0576-3 MUSYEW 1432-1882 Google Scholar

15.

P. Gao et al., “Direct prior regularization from anatomical images for cone beam x-ray luminescence computed tomography reconstruction,” Proc. SPIE, 10578 105781I (2018). https://doi.org/10.1117/12.2293205 PSISDG 0277-786X Google Scholar

16.

P. Gao et al., “Sparse view cone beam x-ray luminescence tomography based on truncated singular value decomposition,” Opt. Express, 26 (18), 23233 –23250 (2018). https://doi.org/10.1364/OE.26.023233 OPEXFF 1094-4087 Google Scholar

17.

T. Liu et al., “Regularized reconstruction based on joint

L_{1}

and total variation for sparse-view cone-beam x-ray luminescence computed tomography,” Biomed. Opt. Express, 10 (1), 1 –17 (2019). https://doi.org/10.1364/BOE.10.000001 BOEICL 2156-7085 Google Scholar

18.

C. M. Carpenter et al., “Limited-angle x-ray luminescence tomography: methodology and feasibility study,” Phys. Med. Biol., 56 (12), 3487 –3502 (2011). https://doi.org/10.1088/0031-9155/56/12/003 PHMBA7 0031-9155 Google Scholar

19.

X. Liu, Q. Liao and H. Wang, “In vivo x-ray luminescence tomographic imaging with single-view data,” Opt. Lett., 38 (22), 4530 –4533 (2013). https://doi.org/10.1364/OL.38.004530 OPLEDP 0146-9592 Google Scholar

20.

X. Liu et al., “A wavelet-based single-view reconstruction approach for cone beam x-ray luminescence tomography imaging,” Biomed. Opt. Express, 5 (11), 3848 –3858 (2014). https://doi.org/10.1364/BOE.5.003848 BOEICL 2156-7085 Google Scholar

21.

X. Liu et al., “Single-view cone-beam x-ray luminescence optical tomography based on Group_YALL1 method,” Phys. Med. Biol., 64 (10), 105004 (2019). https://doi.org/10.1088/1361-6560/ab1819 PHMBA7 0031-9155 Google Scholar

22.

S. Jiang et al., “Reconstruction of fluorescence molecular tomography via a fused LASSO method based on group sparsity prior,” IEEE Trans. Biomed. Eng, 66 (5), 1361 –1371 (2019). https://doi.org/10.1109/TBME.10 IEBEAX 0018-9294 Google Scholar

23.

F. Liu et al., “Weighted depth compensation algorithm for fluorescence molecular tomography reconstruction,” Appl. Opt., 51 (36), 8883 –8892 (2012). https://doi.org/10.1364/AO.51.008883 APOPAI 0003-6935 Google Scholar

24.

F. Tian and H. Liu, “Depth-compensated diffuse optical tomography enhanced by general linear model analysis and an anatomical atlas of human head,” Neuroimage, 85 166 –180 (2014). https://doi.org/10.1016/j.neuroimage.2013.07.016 NEIMEF 1053-8119 Google Scholar

25.

X. Cao et al., “An adaptive Tikhonov regularization method for fluorescence molecular tomography,” Med. Biol. Eng. Comput., 51 (8), 849 –858 (2013). https://doi.org/10.1007/s11517-013-1054-5 MBECDY 0140-0118 Google Scholar

26.

A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci., 2 (1), 183 –202 (2009). https://doi.org/10.1137/080716542 Google Scholar

27.

J. F. P.-J. Abascal et al., “Fluorescence diffuse optical tomography using the split Bregman method,” Med. Phys., 38 (11), 6275 –6284 (2011). https://doi.org/10.1118/1.3656063 MPHYA6 0094-2405 Google Scholar

28.

G.-B. Ye and X. Xie, “Split Bregman method for large scale fused LASSO,” Comput. Stat. Data. Anal., 55 (4), 1552 –1569 (2011). https://doi.org/10.1016/j.csda.2010.10.021 CSDADW 0167-9473 Google Scholar

29.

I. Dassios, K. Fountoulakis and J. Gondzio, “A preconditioner for a primal-dual Newton conjugate gradient method for compressed sensing problems,” SIAM J. Sci. Comput., 37 (6), A2783 –A2812 (2015). https://doi.org/10.1137/141002062 SJOCE3 1064-8275 Google Scholar

Biography

Peng Gao received his MS degree in biomedical engineering from the Fourth Military Medical University of China in 2016, and he is currently pursuing his PhD. His research interests include low-dose CT imaging and multimodality imaging.

Junyan Rong received her PhD in physics from the Institute of High Energy Physics, Chinese Academy of Sciences, Beijing, China, in 2009. Currently, she is an assistant professor in the Department of Biomedical Engineering at the Fourth Military Medical University of China. Her research interests include CT techniques and multimodality biomedical imaging.

Tianshuai Liu received his PhD degree from Fourth Military Medical University of China in 2018. His current research interests are in the field of x-ray imaging model and image reconstruction techniques.

Wenli Zhang received his PhD degree from Louisiana Tech University in 2015. His current research interests are in the field of x-ray-excitable nanophosphors and x-ray induced photodynamic therapy.

Bin Lan is a MD candidate in the Department of Biomedical Engineering at Fourth Military Medical University of China. His current research interests are the synthesis of x-ray-excitable nanophosphors.

Xiaoping Ouyang is an academician of the Chinese Academy of Engineering. He received his PhD in particle physics and nuclear physics from Fudan University in 2002. His current research interests are radiophysical measurement and pulsed neutron diagnosis.

Hongbing Lu is a professor in the Department of Biomedical Engineering at the Fourth Military Medical University of China. She received her PhD in biomedical engineering from Tsinghua University of China, in 1998, followed by postdoctoral training in radiology at the State University of New York at Stony Brook, until 2002. Her research interests cover a spectrum from medical image reconstruction to image analysis for computer-aided detection (CADe) and diagnosis (CADx).

CC BY: © The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.

Citation Download Citation

Peng Gao, Junyan Rong, Tianshuai Liu, Wenli Zhang, Bin Lan, Xiaoping Ouyang, and Hongbing Lu "Limited view cone-beam x-ray luminescence tomography based on depth compensation and group sparsity prior," Journal of Biomedical Optics 25(1), 016004 (22 January 2020). https://doi.org/10.1117/1.JBO.25.1.016004

Received: 8 November 2019; Accepted: 6 January 2020; Published: 22 January 2020

Access the abstract

JOURNAL ARTICLE
14 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 7 scholarly publications.

Explore citations on Lens.org

KEYWORDS

X-rays

Luminescence

Tomography

X-ray imaging

Numerical simulations

Imaging systems

Spatial resolution

1.

Introduction

2.

Methodology

2.1.

Forward Problem

Eq. (1)

Eq. (2)

Eq. (3)

Eq. (4)

2.2.

Limited View Reconstruction with DC-FL for CB-XLCT

Fig. 1

2.2.1.

Depth compensation for the inverse problem

Eq. (5)

Eq. (6)

Eq. (7)

Eq. (8)

2.2.2.

Fused LASSO penalty for the inverse problem

Eq. (9)

Eq. (10)

2.2.3.

Split Bregman method for the fused LASSO problem

3.

Experimental Setup

3.1.

Numerical Simulation Setup

Fig. 2

Fig. 3

3.2.

Phantom Experiments Setup

Fig. 4

3.3.

Quantitative Evaluation

Eq. (11)

Eq. (12)

Eq. (13)

Eq. (14)

4.

Results

4.1.

Numerical Simulations

Fig. 5

Fig. 6

Fig. 7

Table 1

4.2.

Phantom Experiments

Fig. 8

Fig. 9

Table 2

5.

Discussion and Conclusion

Disclosures

Acknowledgments

References

Biography

Show All Keywords

Keywords/Phrases

Search In:

Publication Years