Sparse view cone beam X-ray luminescence tomography based on truncated singular value decomposition

Peng Gao; Junyan Rong; Huangsheng Pu; Tianshuai Liu; Wenli Zhang; Xiaofeng Zhang; Hongbing Lu

doi:10.1364/OE.26.023233

1. Introduction

With the development of X-ray excitable nanophosphors that can produce visible or near-infrared (NIR) light upon X-ray irritation, X-ray luminescence computed tomography (XLCT) has attracted attentions as a new imaging modality [1–3]. It naturally combines X-ray structure imaging of high resolution and optical molecular imaging of high sensitivity and specificity. In addition, XLCT has the advantages of deep penetrating ability with the use of X-rays and higher sensitivity compared with other optical molecular imaging modalities such as fluorescence molecular tomography (FMT) [4, 5] and bioluminescence tomography (BLT) [6, 7] due to the avoidance of background optical signals and autofluorescence.

Since the feasibility of XLCT was first demonstrated by Xing’s group, much efforts have been devoted to XLCT systems including narrow/pencil beam XLCT [3, 8, 9], cone beam XLCT [10, 11] and the corresponding reconstruction methods. The main advantage of narrow/pencil beam XLCT is the high spatial resolution and high imaging depth with X-ray beam’s position as priors. However, the imaging time of narrow/pencil beam is long due to the selective excitation scheme. Though researchers proposed various methods to improve scanning speed with multiple pinhole collimator [12] or focused X-ray beam [13], the imaging time is comparatively longer than cone beam XLCT (CB-XLCT), which is more suitable for practical applications such as small animal whole-body imaging [14–17].

However, the image reconstruction of CB-XLCT is a severely ill-conditioned inverse problem due to the diffuse propagation of luminescence photons in biological tissues, which results in poor image quality and low spatial resolution. To improve the performance of CB-XLCT, wavelet-based [16] reconstruction and multi-wavelength based [15] reconstruction method have been proposed for tomographic imaging of a single luminescent target. To resolve two adjacent targets, an excitation-based method [18] and a reconstruction method using a forward model based on X-ray absorption dosage [19] have been proposed. Though the spatial resolution was improved effectively, data acquisition time was increased greatly in the former method due to multiple acquisition of projections at different X-ray tube voltages, and both data acquisition and preprocessing time was increased in the later method due to the estimation of X-ray absorption dosages. Both the two methods together with the Bayesian framework based CB-XLCT reconstruction methods [20, 21] proposed by Zhang et al require at least 24 luminescent projections to reconstruct. The in vivo applications of CB-XLCT, especially the visualization of some fast biological processes demands fewer imaging projections.

Considering that in most biological applications, the X-ray excitable nanophosphors are sparsely distributed and the luminescent targets are usually small compared with the entire reconstruction region especially when the nanophosphors accumulate in tumor area of early stages, the compressive sensing (CS) technique [22, 23] can be adopted for high quality reconstruction from fewer views. Since the success of CS depends on the restricted isometry property (RIP) of the measurement matrix, which means that the matrix should be incoherent or orthogonal, the transformation of the measurement matrix with a sparse transform such as a truncated singular value decomposition (TSVD), could improve the RIP of CB-XLCT effectively. For fast imaging of multiple targets that is required by in vivo applications, in this paper, (TSVD) based sparse view reconstruction method was proposed for CB-XLCT. To evaluate the performance of the proposed method, both numerical simulation and physical experiments with phantoms of two targets were designed and its performance was compared with three traditional reconstruction methods. It indicates that the proposed method could achieve better results in terms of localization accuracy, target shape, image contrast and spatial resolution with sparse views, making it a practical candidate for fast dynamic CB-XLCT imaging.

This paper is organized as follows. The imaging model and reconstruction method are described in Section 2. In Section 3, both numerical simulations and physical phantom experiments are detailed. Finally, the results are discussed and concluded in Section 4.

2. Method

2.1 Forward problem

In CB-XLCT system, X-rays travel through the imaging object based on Lambert-Beers’ law. The X-ray intensity distribution A(r) in the imaging object can be expressed as follows [10]:

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

where A(r₀) is the intensity of X-ray at the initial position r₀, and μ_t(τ) is the X-ray attenuation coefficient at position τ, which can be obtained from X-ray transmission data. After being irradiated, X-ray excitable nanophosphors distributed in the imaging object can emit visible or NIR light. The process can be expressed as follows,

S (r) = ε A (r) ρ (r)

where S(r) is the source energy density, ε is the light yield of the nanophosphors which is defined as the quantum yield per unit nanophosphor concentration, and ρ(r) represents the nanophosphor concentration in position r, respectively. The photon migration in biological tissues can be modeled using the coupled diffusion equation (DE) [6] considering the highly scattering and weakly absorbing properties of biological tissues in the visible or NIR spectral window:

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

where D(r) is the diffusion coefficient that can be calculated by D(r) = (3(𝜇_a(r) + 𝜇’_s(r)))⁻¹, 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, ∂Ω denotes the boundary of Ω, α is the boundary mismatch factor, and ν is the outward unit normal vector on ∂Ω.

With the finite-element method (FEM), Eq. (3) can be discretized into the following matrix equation

F ε A ρ = K Φ

with

{\begin{cases} F_{i j} = \int_{Ω} Φ (r) ψ_{i} ψ_{j} d r \\ A_{i j} = A_{i j} (r) \\ K_{i j} = \int_{Ω} (D \nabla ψ_{i} \nabla ψ_{j} + μ_{a} ψ_{i} ψ_{j}) d r + \frac{1}{2 α} \int_{\partial Ω} ψ_{i} ψ_{j} d r \end{cases}

where F_ij, K_ij, A_ij are the elements of matrix F, K and A, respectively. ψ_i and ψ_j are corresponding elements of test function of the FEM.

Since the matrix K is positive definite, the above model can be converted into a linear relationship between the nanophosphor concentration ρ and photon measurements on the object surface Φ_meas:

W ρ = Φ_{meas}

where W = K⁻¹FεA, and W is a weight matrix used to map the unknown nanophosphor distribution to known measurements.

2.2 Reconstruction for CB-XLCT

2.2.1 The inverse problem of XLCT

Solving the inverse problem of CB-XLCT is to recover the distribution of nanophosphors ρ from the measured values of the observable measurements Φ_meas. Generally, the weight matrix W is poorly conditioned, and it is impractical to solve ρ directly. To obtain a unique and stable solution, it is necessary to involve the regularization to the object function, and Eq. (6) can be solved by the following minimization problem with the L_p regularization term:

Q = {‖ W ρ - Φ_{meas} ‖}^{2} + λ_{p} {‖ ρ ‖}_{P}^{}

where λ_p is the regularization parameter, and ||ρ||_p is the L_p-norm of ρ. When p = 0, Eq. (7) becomes the L₀-norm regularization problem, which could ensure the sparsest image reconstruction. However, the L₀-minimization is a non-deterministic polynomial hard (NP-hard) problem that is almost impossible to solve. To avoid the disadvantage of unstable reconstruction of L₀-norm, researchers resort to a L₁-norm based reconstruction method, which approximates the L₀-norm relatively well. When p = 2, Eq. (7) becomes the L₂-norm regularization problem, which is commonly called the Tikhonov regularization and usually results in over-smoothed reconstruction.

2.2.2 Compressed sensing for CB-XLCT

In some cases, nanophosphors accumulate in specific areas such as certain organs or tumors which are usually sparse, especially when tumors are in the early stage. As a result, the CS theory is applied to the CB-XLCT reconstruction to obtain a sparse result. The performance of CS theory depends on the RIP of the weight matrix, which means that W is orthogonal or the columns of W are incoherent. Thus, we convert matrix W to a new incoherent matrix V by truncated singular value decomposition (TSVD) [24]. Firstly, the singular value decomposition (SVD) conversion of W is obtained by:

W_{m \times n} = U_{m \times m} Σ_{m \times n} V_{n \times n}^{T} ​

Σ_{m \times n} = [\begin{matrix} d i a g {δ_{i}} & 0 \\ 0 & 0 \end{matrix}]

where U and V are orthogonal matrices, composed of the singular vectors of W. Σ represents a matrix of singular values and δ_i is the i-th singular value of Σ. Substituting Eq. (8) into Eq. (6), Eq. (6) can be rewritten as

U_{m \times m} Σ_{m \times n} V_{n \times n}^{T} ​ ρ_{n \times 1} = Φ_{m \times 1}

To avoid the contamination caused by data errors and rounding errors, small δ_i could be filtered out from matrix Σ and only the first r elements of Σ is reserved. As a result, Eq. (10) can be changed as:

U_{m \times r} Σ_{r \times r} V_{r \times n}^{T} ​ ρ_{n \times 1} = Φ_{m \times 1}

Since the unitary matrix U is orthogonal and the diagonal matrix Σ is invertible, Eq. (11) can be converted into the following form:

V_{r \times n}^{T} ​ ρ_{n \times 1} = Σ_{r \times r}^{- 1} U_{m \times r}^{T} Φ_{m \times 1}

Thus, the primary matrix W_m _× _n has been converted to an incoherent matrix $V_{r \times n}^{T}$ , and at the same time, $Φ_{m \times 1}$ is converted to $Σ_{r \times r}^{- 1} U_{r \times n}^{T} Φ_{m \times 1}$ . As a result, Eq. (12) can be transformed into the following equation:

Q = {‖ V_{r \times n}^{T} ρ_{n \times 1} - Σ_{r \times r}^{- 1} U_{m \times r}^{T} Φ_{m \times 1} ‖}^{2} + λ_{p} {‖ ρ_{n \times 1} ‖}_{P}^{}

The proper value of r is crucial to obtain good performance. In our simulations and phantom experiments, r was set manually by testing a set of parameters. The selection criterion was that δ_r was retained empirically in matrix Σ when δ_r > 2%*δ₁, where all δ were arranged in numerical order and δ₁ is the largest singular value.

2.2.3 Sparse view reconstruction with L1-norm regularization

After the above transformation, Eq. (13) can be solved with CS based reconstruction method such as L₀-norm or L₁-norm regularization algorithm. Considering the solution of L₀-norm method is unstable, in this paper, a widely used L₁-norm regularization method (Fast Iterative Shrinkage-Thresholding Algorithm, FISTA) [25] was used for CS reconstruction.

Let f(x) = || $V_{r \times n}^{T} ρ_{n \times 1} - Σ_{r \times r}^{- 1} U_{r \times n}^{T} Φ_{m \times 1}$ ||² and g(x) = λ_p || $ρ_{n \times 1}$ ||₁, the FISTA method can be implemented as follows:

Initialization: L = L(f), where L is a Lipschitz constant of ∇f; y₁ = x₀ $\in R^{n}$ ; t₁ = 1;

Step k (k $\geq 1$ ):

x_{k} = p_{L} (y_{k})

t_{k + 1} = \frac{1 + \sqrt{1 + 4 t_{k}^{2}}}{2}

y_{k + 1} = x_{k} + (\frac{t_{k} - 1}{t_{k + 1}}) (x_{k} - x_{k - 1})

where

p_{L} = \underset{x}{\arg \min} {g (x) + \frac{L}{2} | | x - (y - \frac{1}{L} \nabla f (y)) | |^{2}}

2.3 Numerical simulation setup

To evaluate the performance of the proposed method, numerical simulations of a phantom with two luminescent targets were first implemented. The geometry configuration of the phantom is shown in Fig. 1. It was a cylinder placed on a rotating stage, with the rotational axis defined as Z axis and the bottom plane set as Z = 0 cm. The diameter and height of the cylindrical phantom were 3.0 cm and 2.5 cm, respectively. In order to approximate high scattering media, the absorption coefficient and reduced scattering coefficient were set to 0.02 cm⁻¹ and 10.0 cm⁻¹, respectively. Two luminescent Y₂O₃:Eu³⁺ targets of 0.4 cm in diameter and 0.5 cm in height were placed in the phantom at a depth of 1.0 cm. The edge-to-edge distance (EED) of the two targets were set to 0.50 cm (case 1), 0.30 cm (case 2) and 0.10 cm (case 3), respectively. The nanophosphors were excited by cone beam X-ray source located at the height of Z = 1.25 cm with tube voltage of 40 kV. Luminescent images were collected at every 15 degree for a full span of 360 degree. To make the simulations more realistic, white Gaussian noise was added to generate noisy boundary measurements with signal-to-noise ratio (SNR) set as 20 dB [5].

Fig. 1 Setup for simulation studies. (a)-(c) The view of XY plane of the three cases. (d)-(f) The view of XZ plane of the three cases. Two luminescent targets of Y₂O₃:Eu³⁺ are placed inside a cylinder phantom. The EED along X axis is 0.50 cm in case 1 (first column), 0.30 cm in case 2 (second column) and 0.10 cm in case 3 (third column).

Download Full Size | PDF

2.4 Phantom experiments setup

Physical phantom experiments were conducted based on a custom-made CB-XLCT system developed in our laboratory. In the system, a micro-focus cone beam X-ray source was adopted to excite the phantom, and luminescent projections were collected by a highly sensitive electron-multiplying charge coupled device (EMCCD) camera positioned at 90° toward the X-ray axis. The transmitted X-ray signals were collected by an X-ray flat panel. Please refer to [18] for a detailed configuration of the system. The phantom was a transparent glass cylinder with a diameter of 3.0 cm and height of 7.0 cm, filled with 1% intralipid and water. The absorption coefficient and the reduced scattering coefficient of the medium was 0.02 cm⁻¹ and 10.0 cm⁻¹, respectively. Two glass tubes of 0.40 cm diameter that contained nanophosphors of Y₂O₃:Eu³⁺ were implanted in the phantom. The concentration of Y₂O₃:Eu³⁺ was 0.1 g/ml. The phantom was fixed on the rotation stage of the system. Three cases of phantom experiments with different EEDs (0.50 cm in case 1, 0.23 cm in case 2 and 0.17 cm in case 3) were performed, as shown in Fig. 2.

Fig. 2 Setup for experiment studies. (a)-(c) Representative X-ray projections of the phantoms with EED = 0.50 cm, 0.23 cm and 0.17 cm, respectively. Regions between the blue and green lines are used for the study. (d)-(f) XCT slice of the phantom, corresponding to the slice indicated by the red line in (a)-(c), respectively. The left and right tubes are named as tube 1 and tube 2. Scale bars in (a) and (d): 1 cm.

Download Full Size | PDF

The X-ray tube voltage and current were set to 40 kV and 1 mA for both luminescence and X-ray imaging, respectively. In luminescence imaging, the EMCCD binning was set to 1 × 1, the EM gain was set to 260, and the integrating time was set to 0.5 s. Luminescent images of 360 degree full view were collected at every 15 degree. In X-ray scanning, the detector binning was set to 2 × 2, and the integrating time of each projection was 0.3 s. 360 X-ray projections were obtained with an interval of 1 degree. The Feldkamp–Davis–Kress (FDK) algorithm was used for CT reconstruction [26]. Figures 2(a)-2(c) 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. 2(d)-2(f). Regions between the blue and green lines were selected for the study.

2.5 Quantitative evaluation

To evaluate the performance of the proposed reconstruction method, several quantitative indexes, including location error (LE) [18], dice similarity coefficient (DICE) [19], contrast-to-noise ratio (CNR) [21] and spatial resolution index (SPI) [27], were calculated.

LE is used to evaluate the localization accuracy of the reconstruction, which is defined as the Euclidean distance between the true and reconstructed target positions:

LE = | | 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 reconstructed and true target regions:

DICE = \frac{2 | R O I_{r} \cap R O I_{t} |}{| R O I_{r} | + | R O I_{t} |}

where ROI_r and ROI_t denote the reconstructed and true target regions, respectively. |·| denotes the number of voxels. Generally, the more the DICE is close to 1, the better the reconstruction is close to the true target.

CNR is used for quantitative evaluation of image contrast:

CNR = \frac{| μ_{R O I} - μ_{B K} |}{{(w_{R O I} σ_{R O I}^{2} + w_{B K} σ_{B K}^{2})}^{1 / 2}}

where the subscript ROI and BK denote the target and background regions of the imaged object, respectively, μ and σ² denote the mean value and variance,

w_{R O I}

and

w_{B K}

are weighting factors calculated by the relative volumes of the ROI and background, respectively.

SPI is a quantitative index to analyze the performance of the algorithms in resolving two targets:

SPI = \frac{ρ_{\max}^{l} - ρ_{valley}^{l}}{ρ_{\max}^{l} - ρ_{min}^{l}}

where ρ^l denotes the value of the profile along a given line on the reconstructed cross-section.

ρ_{max}^{l}

,

ρ_{min}^{l}

and

ρ_{valley}^{l}

are the maximal, minimal and valley value between the two peak values, respectively. The closer SPI is to 1, the better.

In this paper, some widely-used traditional methods were implemented for comparison. The algebraic reconstruction technique (ART) was used for solving Eq. (6) with no regularization. Its relaxation factor was set as 0.1. The FISTA and Adaptive Tikhonov regularization (Adaptik) [28] were adopted for solving Eq. (7) with L₁-norm and L₂-norm regularization, respectively. The regularization parameter were empirically set to 10⁻⁵ and 10⁻⁴, respectively. To further demonstrate the performance of the proposed algorithm, reconstructions with different regularizations such as a Haar wavelet [29] and a TV regularizer [30,31], while using the SVD transformation, were performed on both simulations and phantom experiments. The Haar wavelet regularization imposes a sparsity constraint by penalizing the L₁-norm of the Haar wavelet coefficients ω as follows:

\min_{ρ} {‖ V_{r \times n}^{T} ρ_{n \times 1} - Σ_{r \times r}^{- 1} U_{m \times r}^{T} Φ_{m \times 1} ‖}^{2} + λ_{wt} {‖ ω_{n \times 1} ‖}_{1}

Given the piece-wise constant nature of the reconstructed image, the TV norm of the estimated distribution is added as below:

\min_{ρ} {‖ V_{r \times n}^{T} ρ_{n \times 1} - Σ_{r \times r}^{- 1} U_{m \times r}^{T} Φ_{m \times 1} ‖}^{2} + λ_{TV} {‖ ρ_{n \times 1} ‖}_{TV}

where

{‖ ρ_{n \times 1} ‖}_{TV} = {\int || \nabla ρ (x) ||}_{1} d x

. The Haar wavelet and TV regularization parameter were set to 10⁻³. The iteration numbers were chosen empirically according to the results and to ensure the convergence of the calculation.

3. Results

3.1 Numerical simulations

Figures 3-5 are the reconstruction results of the three cases in numerical simulations. In each case, the two luminescent targets are reconstructed with 24, 4 and 2 projections respectively to demonstrate the performance of the proposed method. All the projections used are evenly sampled over 360 span. The true locations of the two targets are depicted by the yellow circles, and the phantom boundary is depicted by the outer red circle. The height of the cross-sections is 1.25 cm. As shown in Fig. 3, when EED is 0.50 cm, all the methods can resolve the two targets with at least 4 imaging projections. However, the image reconstructed with wavelet and TV regularization from 4 views looks poorer than those reconstructed by other methods in terms of noise and spatial resolution, indicating that the direct use of the two regularizers may not be suitable for the sparse view optical reconstruction. Therefore, in following quantitative analysis and evaluation subsections, it was not included for comparison. It can be seen from Fig. 3(g) that the ART yields distorted target shape for it is sensitive to measurement noise. With the help of L₂-norm regularization, the Adaptik method improves the reconstructions (Fig. 3(h)) in noise suppression but the CNR is lower compared with L₁-norm regularization methods (Fig. 3(i) and Fig. 3(l)). When the number of projections decrease to two, only the proposed method can successfully recover the nanophosphors distribution, and the reconstructed signals obtained by other methods are mainly contaminated by noise.

Fig. 3 CB-XLCT reconstruction results of simulation with EED = 0.5 cm. The first to sixth column are results obtained with the ART, Adaptik, FISTA, wavelet regularization, TV regularization and proposed method. The first to third row are the results reconstructed with 24, 4 and 2 projections, respectively. The yellow circles depict the true locations of the two targets and the red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

Fig. 4 CB-XLCT reconstruction results of simulation with EED = 0.30 cm. The first to fourth column are results obtained with the ART, Adaptik, FISTA and proposed method. The first to third rows are the results reconstructed with 24, 4 and 2 projections, respectively. The yellow circles depict the true locations of the two targets and the red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

Fig. 5 CB-XLCT reconstruction results of simulation with EED = 0.10 cm. The first to fourth column are results obtained with the ART, Adaptik, FISTA and proposed method. The first to third rows are the results reconstructed with 24, 4 and 2 projections, respectively. The yellow circles depict the true locations of the two targets and the red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

Table 1 gives the quantitative analysis results of simulation case 1 with 2 imaging projections (corresponding to the third row in Fig. 3). The proposed method achieves the least LE which means the reconstructed positions of the two targets determined were the closest to the true positions. In addition, the proposed method yields highest DICE and CNR, indicating that the targets are best recovered with the least noise. Besides, the SPI obtained with our method is much higher than those obtained by the other three methods, indicating that the proposed method achieves much better separation of two targets. The quantitative results further confirm the observations in Fig. 3.

Table 1. Quantitative analysis of simulation case 1 with 2 projections

View Table | View all tables in this article

Figure 4 shows the reconstruction results of numerical simulation case 2 where the EED is 0.30 cm. It can be seen from Fig. 4(e) that the ART method is unable to resolve the two targets when 4 projections are used. With the regularization strategy, the Adaptik method obtains better resolution (Fig. 4(f)) compared with the ART method, but the results suffer from noise. Thanks to the CS theory, both FISTA and the proposed method can obtain good performance as shown in Fig. 4(g)-4(h). However, when only 2 projections were used, images reconstructed by the ART, Adaptik and FISTA methods are filled with noise (Fig. 4(i)-4(k)) and can hardly locate the targets. By contrast, our method can achieve high quality which can be further demonstrated by quantitative analysis results in Table 2. Similar to case 1, the proposed method yielded the smallest LE, highest CNR, Dice and SPI, indicating that the targets were recovered more accurately.

Table 2. Quantitative analysis of simulation case 2 with 2 projections

View Table | View all tables in this article

Figure 5 shows the reconstruction results of numerical simulation case 3 where the EED is 0.10 cm. All the three traditional methods can hardly separate the two targets when 4 projections were used, and the results with 2 projections are similar as those in Fig. 3 and Fig. 4, where images are full of noises. Quantitative analysis of simulation case 3 with 2 projections is shown in Table 3, which demonstrates that our method outperforms the other three.

Table 3. Quantitative analysis of simulation case 3 with 2 projections

View Table | View all tables in this article

To demonstrate the performance of our method for the reconstruction of all different slices, 3-D rendering of the CB-XLCT results obtained by four methods are shown in Fig. 5. Images in the first to the third row in Fig. 6 are results reconstructed with 2 projections of case 1 to case 3, respectively. Values less than 10% of the maximum value are not considered. The red objects represent the luminescent targets. It demonstrates that compared with other three methods, two targets can be better recovered with accurate position and less noise by the proposed method (Fig. 6(d), 6(h) and 6(l)).

Fig. 6 3-D rendering of the simulation reconstruction results from 2 projections. Images in the first to the third row are results reconstructed with 2 projections of case 1 to case 3, respectively. The red objects represent the recovered targets.

Download Full Size | PDF

3.2 Phantom experiments

Figure 7 shows the cross-section results of phantom experiment case 1, which are fused with the same corresponding CT slice as the red line indicated in Fig. 2(a). Reconstructions obtained by the ART, Adaptik, FISTA, wavelet regularization, TV regularization and the proposed method are shown from the first to the sixth column. The first to the third row represent the reconstruction results obtained with 24 projections, 8 projections, and 4 projections, respectively.

Fig. 7 CB-XLCT fusion results of phantom experiments with EED = 0.50 cm. The first to sixth column are results obtained with the ART, Adaptik, FISTA, wavelet regularization, TV regularization and proposed method. The first to third rows are the results reconstructed with 24, 8 and 4 projections, respectively. The slice is indicated by the red line in Fig. 2(a). The red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

When 24 projections were used for reconstruction, most tested methods could distinguish two targets effectively (Figs. 7(a)-7(d), and Fig. 7(f)), except for the TV regularization (Fig. 7(e)). The result with wavelet regularization has larger location error and the intensity values of the two targets differ greatly. Compared with the other four methods, both the TV regularization and wavelet regularization lead to poor image quality, therefore they were not included for comparison in the following reconstructions and quantitative analysis. The ART reconstruction result has a low contrast due to noise (Fig. 7(a)), this is mainly because no regularization is adopted in the ART method. When the projection number is 8, only the L₁-norm regularization method can resolve the two targets clearly with high contrast and less location error (Fig. 7(i) and Fig. 7 (l)). When the projection number decreases to 4, the ART, Adaptik and FISTA methods are all unable to resolve the targets, and the results contain much noise (Figs. 7(m)-7(o)). Image in Fig. 7(r) shows that our method outperforms the other methods in recovering spatial resolution and image contrast.

Table 4 shows the quantitative analysis of phantom experiment case 1 with 4 projections (third row in Fig. 7) to better clarify the proposed algorithm. Similar to results in simulations, the proposed method yielded the smallest LE, highest CNR, Dice and SPI, indicating that the targets have been recovered with the least relative error.

Table 4. Quantitative analysis of phantom experiment case 1 with 4 projections

View Table | View all tables in this article

Figure 8 shows the cross-section results of phantom experiment case 2, which are fused with the same corresponding CT slice as the red line indicated in Fig. 2(b). Reconstructions obtained by the ART, Adaptik, FISTA and proposed method are shown from the first to the fourth column. The first to the third row represent the reconstruction results obtained with 24 projections, 8 projections, and 4 projections, respectively.

Fig. 8 CB-XLCT fusion results of phantom experiments with EED = 0.23 cm. The first to fourth column are results obtained with the ART, Adaptik, FISTA and proposed method. The first to third row are the results reconstructed with 24, 8 and 4 projections, respectively. The slice is the same as that indicated by the red line in Fig. 2(b). The red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

Unlike case 1, in this case, the Adaptik and FISTA method cannot resolve the two targets even if 24 projections are used, for the two tubes are close to each other. Though two targets have been reconstructed by the ART method as shown in Fig. 8(a), the localization error is large and the shape is distorted. However, Fig. 8(d), 8(h) and 8(l) show that our method can recover the two targets no matter with 24, 8 or 4 imaging projections. Quantitative analysis in Table 5 further illustrates that our method performs better than the other three methods with least LE, highest DICE, CNR and SPI.

Table 5. Quantitative analysis of phantom experiment case 2 with 4 projections

View Table | View all tables in this article

Figure 9 shows the tomographic results of phantom experiment case 3, which are fused with the same corresponding CT slice as the red line indicated in Fig. 2(c). The results are similar with those in Fig. 8, which means traditional method are unable to resolve close targets, and our method can achieve the goal even with 4 projections. Quantitative analysis of phantom experiments case 3 with 2 projections is shown Table 6, which further indicates the effectiveness of the proposed method.

Fig. 9 CB-XLCT fusion results of phantom experiments with EED = 0.17 cm. The first to fourth column are results obtained with the ART, Adaptik, FISTA and proposed method. The first to third row are the results reconstructed with 24, 8 and 4 projections, respectively. The slice is the same as that indicated by the red line in Fig. 2(c). The red circles depict the boundary of the phantom. All images are normalized to the maximal value.

Download Full Size | PDF

Table 6. Quantitative analysis of phantom experiment case 3 with 4 projections

View Table | View all tables in this article

To demonstrate the performance of our method for phantom experiments at all slices, 3-D rendering of the CB-XLCT results obtained by four methods are shown in Fig. 10. Images in the first to the third in Fig. 10 are results reconstructed with 4 projections of case 1 to case 3, respectively. Values less than 10% of the maximum value are not considered. The red objects represent the luminescent targets. Similar to the results of simulations, the two targets can be recovered best by the proposed method (Figs. 10(d)-10(h)).

Fig. 10 3-D rendering of the simulation reconstruction results with 4 projections. Images in the first to the third row are results reconstructed with 4 projections of case 1 to case 3, respectively. The red objects represent the recovered targets.

Download Full Size | PDF

4. Discussion and conclusion

As a novel optical molecular technique, CB-XLCT could provide fast and dual-modality imaging and plays an important role in 3-D non-invasive imaging of biological and biochemical processes in living subjects. However, due to severe ill-posedness of CB-XLCT inverse problem, the image quality of the reconstruction is relatively low. It is worse when sparse projections are used. This limits the application of CB-XLCT especially the fast dynamic imaging in vivo. To address the problem, in this study, we have established a compressive sensing-based sparse view CB-XLCT reconstruction method. To meet the demand of CS theory, the weight matrix was transferred to an incoherent matrix by TSVD, and then an L₁-norm regularization based method (FISTA) was adopted to solve the new linear equation.

The performance of the proposed method is demonstrated by numerical simulations and phantom experimental results (Figs. 3-10). The reconstruction results show that the proposed method could resolve targets with an EED of 0.10 cm using only 2 imaging projections in simulation studies and an EED of 0.17 cm using 4 projections in phantom experiments. Quantitative analysis for XLCT reconstructions indicates that compared with ART, Adaptik and FISTA methods, the proposed method could achieve better results in location accuracy, shape recovery, image contrast and spatial resolution of the luminescent targets (Tables 1-6), especially when sparse views were used. All these results demonstrate the potential of the proposed method in improving CB-XLCT for fast and dynamic imaging. In addition, the proposed algorithm could be further extended to solving other similar inverse problems such as FMT and BLT.

Our results indicate that in phantom experiments, more projection data were required to resolve adjacent targets when compared with simulation studies. That may be due to several reasons. Firstly, the optical diffusion coefficients used for phantom experiments are set according to references, which may not be optimal for our system. With accurate coefficients obtained by diffuse optical tomography (DOT) [32], the reconstruction performance can be further improved. Secondly, the noise in the measured experimental data has not been modeled in this paper based on the fact that the measured luminescent signal is much higher than noise signal in our study. When the luminescent signal is weak, noise model such as Gaussian distribution model must be constructed. Besides, the effect of X-ray beam hardening and scattering is not considered for the size and structure of the phantom is small and simple. In addition, the tubes are longer in the z dimension than the proposed method is able to reconstruct, especially when the EED is small. We have shown that different sparse basis (Haar wavelet and TV regularization) have a great influence on the image quality. We believe that the accuracy in the z dimension can be improved with some modification to the proposed method such as by using a different sparse basis. These factors will be included to construct a more precise imaging model in our future work.

It should be noted that it is difficult to compare the performance of different reconstruction algorithms when considering the effect of regularization parameters, iteration numbers, and initial values on image reconstruction. In this paper, the relaxation factor for the ART method, and the regularization parameters for the Tikhonov, FISTA, and the proposed method were manually optimized according to value ranges suggested by references [28, 33, 34]. Automatic selection of the optimal or proper regularization parameters and iteration numbers by means of Bayesian framework will be our future work.

For optical reconstruction, the image quality benefits from a larger number of projections as the number increase in a certain range, while the increase would be limited beyond that range [35]. In our phantom experiments, the range is from 2 to 4. When the number of projections becomes lager, there is no major quality improvements. In can be deduced that 4 projections contain most of the required reconstruction information in that situation, and more views means redundant measurements are taken when a CS based method is used. In this paper, all the projections selected were evenly sampled over the 360° span. We found that in phantom experiments, reconstruction performance with 3 views was unstable and influenced by the selection of the first view. This was also happened when the 4 views used were not evenly distributed over 360 degrees. In some practical applications, the imaging object is unable to rotate 360 degrees or luminescence can only transmit from one side of the imaging object, making the optical signals be detected only in limited views [36]. In these cases, the proposed method should be further improved to make it more suitable for practical use.

In conclusion, we have presented TSVD based sparse view reconstruction method for CB-XLCT. Both simulations and phantom experiments of multi-targets demonstrates the advantage of our method compared with conventional methods. Future work will focus on applying the proposed method to in vivo experiments.

Funding

National Key Research and Development Program of China (2017YFC0107400, 2017YFC0107401, 2017YFC0107402, 2017YFC0107403, 2017YFC0107404 and 2017YFC0107405); National Natural Science Foundation of China (NSFC) (31700865, 11805274, 61871383); Natural Nature Science Foundation of Shaanxi Province (2016JQ1012); Key Science and Technology Program in Social Development of Shaanxi Province (2016SF-044).

Disclosures

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

References

1. C. M. Carpenter, C. Sun, G. Pratx, R. Rao, and L. Xing, “Hybrid x-ray/optical luminescence imaging: characterization of experimental conditions,” Med. Phys. 37(8), 4011–4018 (2010). [CrossRef] [PubMed]

2. G. Pratx, C. M. Carpenter, C. Sun, R. P. Rao, and L. Xing, “Tomographic molecular imaging of x-ray-excitable nanoparticles,” Opt. Lett. 35(20), 3345–3347 (2010). [CrossRef] [PubMed]

3. G. Pratx, C. M. Carpenter, C. Sun, and L. Xing, “X-ray luminescence computed tomography via selective excitation: a feasibility study,” IEEE Trans. Med. Imaging 29(12), 1992–1999 (2010). [CrossRef] [PubMed]

4. H. Yi, D. Chen, W. Li, S. Zhu, X. Wang, J. Liang, and J. Tian, “Reconstruction algorithms based on l1-norm and l2-norm for two imaging models of fluorescence molecular tomography: a comparative study,” J. Biomed. Opt. 18(5), 056013 (2013). [CrossRef] [PubMed]

5. H. Pu, G. Zhang, W. He, F. Liu, H. Guang, Y. Zhang, J. Bai, and J. Luo, “Resolving fluorophores by unmixing multispectral fluorescence tomography with independent component analysis,” Phys. Med. Biol. 59(17), 5025–5042 (2014). [CrossRef] [PubMed]

6. W. Cong, G. Wang, D. Kumar, Y. Liu, M. Jiang, L. Wang, E. Hoffman, G. McLennan, P. McCray, J. Zabner, and A. Cong, “Practical reconstruction method for bioluminescence tomography,” Opt. Express 13(18), 6756–6771 (2005). [CrossRef] [PubMed]

7. G. Wang, W. Cong, K. Durairaj, X. Qian, H. Shen, P. Sinn, E. Hoffman, G. McLennan, and M. Henry, “In vivo mouse studies with bioluminescence tomography,” Opt. Express 14(17), 7801–7809 (2006). [CrossRef] [PubMed]

8. C. Li, A. Martínez-Dávalos, and S. R. Cherry, “Numerical simulation of x-ray luminescence optical tomography for small-animal imaging,” J. Biomed. Opt. 19(4), 046002 (2014). [CrossRef] [PubMed]

9. M. C. Lun, W. Zhang, and C. Li, “Sensitivity study of x-ray luminescence computed tomography,” Appl. Opt. 56(11), 3010–3019 (2017). [CrossRef] [PubMed]

10. D. Chen, S. Zhu, H. Yi, X. Zhang, D. Chen, J. Liang, and J. Tian, “Cone beam x-ray luminescence computed tomography: A feasibility study,” Med. Phys. 40(3), 031111 (2013). [CrossRef] [PubMed]

11. X. Liu, Q. Liao, and H. Wang, “Fast X-ray luminescence computed tomography imaging,” IEEE Trans. Biomed. Eng. 61(6), 1621–1627 (2014). [CrossRef] [PubMed]

12. W. Zhang, D. Zhu, M. Lun, and C. Li, “Multiple pinhole collimator based X-ray luminescence computed tomography,” Biomed. Opt. Express 7(7), 2506–2523 (2016). [CrossRef] [PubMed]

13. W. Zhang, M. C. Lun, A. A.-T. Nguyen, and C. Li, “X-ray luminescence computed tomography using a focused x-ray beam,” J. Biomed. Opt. 22(11), 1–11 (2017). [CrossRef] [PubMed]

14. 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). [CrossRef] [PubMed]

15. D. Chen, S. Zhu, X. Chen, T. Chao, X. Cao, F. Zhao, L. Huang, and J. Liang, “Quantitative cone beam X-ray luminescence tomography/X-ray computed tomography imaging,” Appl. Phys. Lett. 105(19), 191104 (2014). [CrossRef]

16. X. Liu, H. Wang, M. Xu, S. Nie, and H. Lu, “A wavelet-based single-view reconstruction approach for cone beam x-ray luminescence tomography imaging,” Biomed. Opt. Express 5(11), 3848–3858 (2014). [CrossRef] [PubMed]

17. D. Chen, S. Zhu, X. Cao, F. Zhao, and J. Liang, “X-ray luminescence computed tomography imaging based on X-ray distribution model and adaptively split Bregman method,” Biomed. Opt. Express 6(7), 2649–2663 (2015). [CrossRef] [PubMed]

18. P. Gao, H. Pu, J. Rong, W. Zhang, T. Liu, W. Liu, Y. Zhang, and H. Lu, “Resolving adjacent nanophosphors of different concentrations by excitation-based cone-beam X-ray luminescence tomography,” Biomed. Opt. Express 8(9), 3952–3965 (2017). [CrossRef] [PubMed]

19. T. Liu, J. Rong, P. Gao, W. Zhang, W. Liu, Y. Zhang, and H. Lu, “Cone-beam x-ray luminescence computed tomography based on x-ray absorption dosage,” J. Biomed. Opt. 23(2), 1–11 (2018). [PubMed]

20. G. Zhang, S. Tzoumas, K. Cheng, F. Liu, J. Liu, J. Luo, J. Bai, and L. Xing, “Generalized Adaptive Gaussian Markov Random Field for X-ray Luminescence Computed Tomography,” IEEE Trans. Biomed. Eng. 65, 2130–2133 (2017). [PubMed]

21. G. Zhang, F. Liu, J. Liu, J. Luo, Y. Xie, J. Bai, and L. Xing, “Cone Beam X-ray Luminescence Computed Tomography Based on Bayesian Method,” IEEE Trans. Med. Imaging 36(1), 225–235 (2017). [CrossRef] [PubMed]

22. E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag. 25(2), 21–30 (2008). [CrossRef]

23. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]

24. J. Shi, X. Cao, F. Liu, B. Zhang, J. Luo, and J. Bai, “Greedy reconstruction algorithm for fluorescence molecular tomography by means of truncated singular value decomposition conversion,” J. Opt. Soc. Am. A 30(3), 437–447 (2013). [CrossRef] [PubMed]

25. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Imaging Sci. 2(1), 183–202 (2009). [CrossRef]

26. L. Feldkamp, L. Davis, and J. Kress, “Practical cone-beam algorithm,” J. Opt. Soc. Am. A 1(6), 612–619 (1984). [CrossRef]

27. J. Shi, F. Liu, G. Zhang, J. Luo, and J. Bai, “Enhanced spatial resolution in fluorescence molecular tomography using restarted L1-regularized nonlinear conjugate gradient algorithm,” J. Biomed. Opt. 19(4), 046018 (2014). [CrossRef] [PubMed]

28. X. Cao, B. Zhang, X. Wang, F. Liu, K. Liu, J. Luo, and J. Bai, “An adaptive Tikhonov regularization method for fluorescence molecular tomography,” Med. Biol. Eng. Comput. 51(8), 849–858 (2013). [CrossRef] [PubMed]

29. J. Verhaeghe, D. Van de Ville, I. Khalidov, Y. D’Asseler, I. Lemahieu, and M. Unser, “Dynamic PET reconstruction using wavelet regularization with adapted basis functions,” IEEE Trans. Med. Imaging 27(7), 943–959 (2008). [CrossRef] [PubMed]

30. Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, “Low-dose CT reconstruction via edge-preserving total variation regularization,” Phys. Med. Biol. 56(18), 5949–5967 (2011). [CrossRef] [PubMed]

31. A. Behrooz, H. M. Zhou, A. A. Eftekhar, and A. Adibi, “Total variation regularization for 3D reconstruction in fluorescence tomography: experimental phantom studies,” Appl. Opt. 51(34), 8216–8227 (2012). [CrossRef] [PubMed]

32. D. A. Boas, D. H. Brooks, E. L. Miller, C. A. DiMarzio, M. Kilmer, R. J. Gaudette, and Q. Zhang, “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag. 18(6), 57–75 (2001). [CrossRef]

33. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, A. Goy, C. Vonesch, M. Unser, and D. Psaltis, “Optical tomographic image reconstruction based on beam propagation and sparse regularization,” IEEE Trans. Comput. Imag. 2(1), 59–70 (2016). [CrossRef]

34. X. Liu, F. Liu, Y. Zhang, and J. Bai, “Unmixing dynamic fluorescence diffuse optical tomography images with independent component analysis,” IEEE Trans. Med. Imaging 30(9), 1591–1604 (2011). [CrossRef] [PubMed]

35. X. Cao, B. Zhang, F. Liu, X. Wang, and J. Bai, “Reconstruction for limited-projection fluorescence molecular tomography based on projected restarted conjugate gradient normal residual,” Opt. Lett. 36(23), 4515–4517 (2011). [CrossRef] [PubMed]

36. C. M. Carpenter, G. Pratx, C. Sun, and L. Xing, “Limited-angle x-ray luminescence tomography: methodology and feasibility study,” Phys. Med. Biol. 56(12), 3487–3502 (2011). [CrossRef] [PubMed]

Method	LE (mm)		DICE		CNR	SPI
Method	tube 1	tube 2	tube 1	tube 2	CNR	SPI
ART	4.62	3.42	0.07	0.07	1.64	0.56
Adaptik	4.44	3.65	0.08	0.08	1.80	0.65
FISTA	4.57	4.20	0.07	0.07	2.05	0.54
Proposed method	0.10	1.49	0.79	0.60	6.52	0.91

Method	LE (mm)		DICE		CNR	SPI
Method	tube 1	tube 2	tube 1	tube 2	CNR	SPI
ART	3.80	3.40	0.08	0.08	2.25	0.39
Adaptik	4.69	4.10	0.08	0.08	1.73	0.39
FISTA	4.84	3.97	0.08	0.08	2.34	0.36
Proposed method	0.20	0.32	0.55	0.69	9.21	0.72

Method	LE (mm)		DICE		CNR	SPI
Method	tube 1	tube 2	tube 1	tube 2	CNR	SPI
ART	2.93	2.04	0.05	0.05	2.36	0.39
Adaptik	2.88	4.08	0.05	0.05	1.93	0.34
FISTA	2.72	3.97	0.06	0.06	1.89	0.40
Proposed method	0.41	1.10	0.60	0.74	8.38	0.75

Method	LE (mm)		DICE		CNR	SPI
Method	tube 1	tube 2	tube 1	tube 2	CNR	SPI
ART	5.56	5.22	0.18	0.18	1.89	0.06
Adaptik	6.23	4.34	0.31	0.21	3.66	0.81
FISTA	5.53	5.12	0.34	0.27	4.31	0.60
Proposed method	1.28	0.96	0.68	0.61	6.33	1.00

Method	LE (mm)		DICE		CNR	SPI
Method	tube 1	tube 2	tube 1	tube 2	CNR	SPI
ART	4.22	3.34	0.13	0.28	1.88	0.29
Adaptik	3.14	4.11	0.20	0.33	2.82	0.20
FISTA	3.38	3.56	0.38	0.34	2.66	0.05
Proposed method	1.79	1.15	0.46	0.67	3.41	0.86

Sparse view cone beam X-ray luminescence tomography based on truncated singular value decomposition

Abstract

1. Introduction

2. Method

2.1 Forward problem

2.2 Reconstruction for CB-XLCT

2.2.1 The inverse problem of XLCT

2.2.2 Compressed sensing for CB-XLCT

2.2.3 Sparse view reconstruction with L1-norm regularization

2.3 Numerical simulation setup

2.4 Phantom experiments setup

2.5 Quantitative evaluation

3. Results

3.1 Numerical simulations

3.2 Phantom experiments

4. Discussion and conclusion

Funding

Disclosures

References

Cited By

Figures (10)

Tables (6)

Equations (22)

Optics Express