Regularized reconstruction based on joint L<sub>1</sub> and total variation for sparse-view cone-beam X-ray luminescence computed tomography

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

doi:10.1364/BOE.10.000001

1. Introduction

With the advances of X-ray excitable nanophosphors, X-ray luminescence computed tomography (XLCT) has attracted more attention for its promising performance [1,2]. In XLCT, X-ray excitable nanophosphors are used as imaging probes and emit visible or near-infrared (NIR) light when irradiated by X-rays which penetrates the object to be imaged and can be detected by sensitive photon detectors. By solving an inverse problem using an appropriate imaging model of X-ray and photon transport, the three-dimensional (3-D) distribution of nanophosphors in the imaged object can be resolved. Compared with traditional bio-optical imaging modalities such as bioluminescence tomography (BLT) [3,4] and fluorescence molecular tomography (FMT) [5,6], XLCT has the following advantages. Firstly, due to the use of X-ray excitation, the interference of autofluorescence and background fluorescence can be avoided, which can improve the contrast and resolution of imaging. Secondly, because of the high penetrability and collimation of X-rays, the depth of imaging can be improved, and the deep tissue imaging is expected to be realized. Thirdly, the X-ray CT of high resolution imaging and the optical molecular tomography with high sensitivity can be obtained simultaneously. Therefore, XLCT has become a promising imaging technique for fundamental research, drug development, and clinical studies.

After the first demonstration of XLCT, systems with different imaging geometries have been proposed to improve the spatial and temporal resolution. Compared with the narrow-beam and fan-beam XLCT systems, cone-beam XLCT speed up the imaging process significantly at the cost of spatial resolution [7,8]. To improve its reconstruction quality, Zhang et al proposed a self-adaptive Bayesian method for CB-XLCT reconstruction and validated its superiority with numerical simulations and mouse experiments [9]. For fast XLCT imaging, Liu et al combined the compressive sensing (CS) technique with a wavelet transform to implement CB-XLCT imaging of a single target from single-view data [10]. However, due to highly scattered nature of light through biological tissues, fast reconstruction of multiple targets from few-view detection data is more ill-conditioned, limiting the application of CB-XLCT in in-vivo imaging.

Considering that the distribution of nanophosphors inside bodies preferentially accumulates in specific areas of interest, the reconstruction of XLCT image is usually sparse with some locally smoothed high-intensity regions. It indicates that a combination of regularizations enforcing the conditions of sparsity and smoothness could further improve the reconstruction of multi-targets CB-XLCT. For sparsity consideration, the most basic sparse index is the L₀ norm, which is the number of non-zero elements in all vectors. However, the normal inverse problem of L₀, at least in the case of undetermined, is NP-hard and difficult to be solved [11]. In contrast, L₁ norm is a convex relaxation of L₀ norm and often used to enhance the sparsity of images, especially in the field of compressed sensing [12,13]. Corresponding optimization problem can be solved by many algorithms, such as gradient projection (GP) algorithm [14], iterative shrinkage-thresholding (IST) algorithm [15], fast iterative shrinkage-thresholding (FISTA) algorithm [16–18] etc. In this study, we choose L₁ penalty to enhance the sparsity of XLCT images and use the FISTA algorithm, which could get fast and accurate L₁ solutions, to solve the reconstruction problem.

A series of smooth priors have been applied in the reconstruction of tomography [19]. Among them, the prior constraints of quadratic terms are widely used. However, the prior restraint of the quadratic term is more inclined smooth out edges in images. Therefore, in this study, the total variation (TV) penalty [20], is designed to promote the smoothness while preserving the edges in the image. A variety of methods can be used to deal with TV penalty, such as gradient descent (GD) [21,22] and separable paraboloidal surrogates (SPS) [23,24] algorithm. In this study, the traditional GD method is adopted to minimize the TV constraint in XLCT reconstruction.

In this paper, we propose a reconstruction approach based on joint L₁ and total variation regularization for the CB-XLCT reconstruction. The remainder of this paper is organized as follows. In Section 2, the imaging model, the proposed L₁-TV regularization and how to solve the algorithm are described in detail. In Section 3, both numerical simulations and phantom experiments are detailed. In Section 4, the numerical simulations and phantom experiments results are described for the performance evaluation of the proposed reconstruction approach. Finally, discussions and conclusions are given in Section 4.

2. Methods

2.1 Forward model of XLCT

For XLCT imaging, when irradiated by X-rays, nanophosphors in the object can emit visible or NIR light. Based on the previous studies, the number of optical photons emitted is proportional to the intensity distribution of the X-rays and the concentration of nanophosphor in the object, which can be expressed as [25]:

S (r) = Γ X (r) n (r)

where S(r) is the light emitted, n(r) is the concentration of nanophosphors, Γ is the light yield of the nanophosphors, and X(r) is the intensity of X-rays at position r, which can be given by the Lambert-Beer law [9].

In the visible and NIR spectral window, biological tissues have the characteristics of high scattering and low absorbing. Therefore, the propagation model of the emitted light in biological tissues can be established by the diffusion equation (DE) [26]:

- \nabla \cdot [D (r) \nabla Φ (r)] + μ_{a} (r) Φ (r) = S (r) (r \in Ω)

where Ω is the image domain, Φ(r) is the photon fluence, μ_a(r) is the absorption coefficient. D(r) represents the diffusion coefficient that can be calculated by

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

, in which

μ_{s}^{'} (r)

is the reduced scattering coefficient.

To solve the diffusion Eq. (2), the Robin boundary conditions are usually applied [27,28], as shown below:

Φ (r) + 2 κ D (r) [ν \nabla Φ (r)] = 0 (r \in \partial Ω)

where

\partial Ω

is the boundary of

Ω

,

κ

is the boundary mismatch parameter and

ν

represents the outward unit normal vector on the boundary.

With the finite element method (FEM), Eq. (2) and Eq. (3) can be discretized into a matrix equation as:

A Φ = Γ F N

with

\begin{array}{l} a_{i j} = \int_{Ω} [D (r) \nabla ψ_{i} (r) \nabla ψ_{j} (r) + μ_{a} (r) ψ_{i} (r) ψ_{j} (r)] d r \\ + \frac{1}{2 κ} \int_{\partial Ω} ψ_{i} (r) ψ_{j} (r) d r \end{array}

f_{i j} = \int_{Ω} X (r) ψ_{i} (r) ψ_{j} (r) d r,

Here N is the distribution vector of nanophosphors, a_ij and f_ij are the elements of matrix A and F, respectively. ψ_i(r) and ψ_j(r) are the corresponding elements of discretized geometrical meshes of the imaging domain, and X(r) is the intensity of X-rays at position r.

Since the matrix A is positive definite, Eq. (4) can be further recast into

Φ = M N

where

M = Γ A^{- 1} F

, Φ represents the distribution vector of photon fluence. For optical tomography, only intensity values of Φ on the object surface could be measured, then Eq. (7) becomes

Φ_{m e a s} = W N

where Φ_meas is the vector of photon fluence acquired on the object surface, and W consists of rows of the weight matrix M corresponding to surface measurements.

2.2 Inverse problem

The goal of the XLCT reconstruction is to estimate the nanophosphor distribution N from Φ_meas. In practical application of XLCT, noise of the XLCT imaging system needs to be considered, and Eq. (8) becomes:

b = Φ_{m e a s} + ς = W x + ς

where b represents the actual fluorescence signals measured on the surface of the imaging object,

ς

is the noise of the system, x = N represents the unknown distribution of nanophosphors in the imaging object.

The objective function of the reconstructing method based on the joint L₁ and total variation regularization can be expressed as:

x = \arg \min_{x \geq 0} \frac{1}{2} {‖ W x - b ‖}_{2}^{2} + λ_{L_{^{1}}} {‖ x ‖}_{1}^{} + λ_{T V} {| x |}_{T V}

where

λ_{L_{^{1}}}

and

λ_{T V}

are the regularization parameters.

2.3 Solution of the inverse problem

The solution to Eq. (10) is divided into two parts. First, the L₁ regular solution of reconstructed image is obtained by the FISTA algorithm. Then the gradient descent (GD) strategy is applied to deal with the TV constraint in the reconstruction of XLCT. The flowchart of the proposed method is recapitulated in Fig. 1. Considering that x represents the distribution of nanophosphors, the non-negative constraint is applied before the TV minimization of x. The implementation process of these two parts is described as below, respectively.

Fig. 1 The flowchart of the proposed method.

Download Full Size | PDF

2.3.1 FISTA for solving the L₁ regularization problem

Let $f (x) = \frac{1}{2} {‖ W x - b ‖}_{2}^{2}$ , $g (x) = {‖ x ‖}_{1}^{}$ , the optimization problem of the L₁ regularization term is changed to:

\min {F (x) \equiv f (x) + g (x) : x \in R^{n}}

In this study, the FISTA algorithm is used to solve Eq. (11) and a simple flowchart is summarized in Algorithm 1.

Algorithm 1. FISTA-L₁

View Table | View all tables in this article

2.3.2 GD for TV constraint in the reconstructed image

The second part of the optimized problem (Eq. (10) can be expressed as the TV minimization problem, which is shown as:

\min_{x} ({‖ x ‖}_{T V}) s . t . \frac{1}{2} {‖ W x - b ‖}_{2}^{2} + λ_{L_{^{1}}} {‖ x ‖}_{1}^{} \leq ε x \geq 0

{‖ x ‖}_{T V} = \sum_{i, j} \sqrt{{(x_{i, j} - x_{i, j - 1})}^{2} + {(x_{i, j} - x_{i - 1, j})}^{2}}

After getting the current L₁ regularization solution in the first part, the GD method is used to solve Eq. (12), in which the derivative of the image TV with respect to each pixel is approximately written as:

\begin{array}{l} \frac{\partial {‖ x ‖}_{T V}}{\partial x_{i, j}} \approx \frac{2 (x_{i, j} - x_{i, j - 1}) + 2 (x_{i, j} - x_{i - 1, j})}{\sqrt{δ + {(x_{i, j} - x_{i, j - 1})}^{2} + {(x_{i, j} - x_{i - 1, j})}^{2}}} \\ \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} - \frac{2 (x_{i + 1, j} - x_{i, j})}{\sqrt{δ + {(x_{i + 1, j} - x_{i, j})}^{2} + {(x_{i + 1, j} - x_{i + 1, j - 1})}^{2}}} \\ \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} - \frac{2 (x_{i, j + 1} - x_{i, j})}{\sqrt{δ + {(x_{i, j + 1} - x_{i, j})}^{2} + {(x_{i, j + 1} - x_{i - 1, j + 1})}^{2}}} \end{array}

where

δ

is a small positive number in order to avoid the denominator to zero.

3. Experimental design

Numerical simulations and phantom experiments were performed to evaluate the performance of the proposed method based on the custom-developed CB-XLCT system in our laboratory. For comparison, three traditional methods, algebra reconstruction technique (ART), adaptive tikhonov regularization (ADAPTIK), and fast iterative shrinkage-thresholding (FISTA) algorithm, were also implemented for solving Eq. (9) with no regularization, L₂ regularization and L₁ regularization, respectively. The relaxation factor was set as 0.1 for ART algorithm. For ADAPTIK and FISTA algorithm, the regularization parameters were empirically set to 10⁻⁴ and 10⁻⁵, respectively. In this study, to make the calculating process consistent for a fair comparison between the FISTA algorithm and the proposed algorithm, the L₁ regularization parameter was set as 10⁻⁵ in the proposed algorithm. For the TV regularization, the gradient descent (GD) method was used to solve the minimization problem, in which the step length and total descent numbers were set as 0.1 and 10, respectively. For ART, ADAPTIK, FISTA and the proposed L₁-TV algorithm, the iterative numbers were empirically set as 3000, 20, 3000, 30, respectively.

3.1 Numerical simulations setup

Firstly, numerical simulations were carried out with a cylinder phantom. The cylinder phantom (Fig. 2) was composed of a large cylinder tank (3.0cm in diameter and 2.3cm in height) and two small tubes (4mm in diameter and 4mm in height) filled with Y₂O₃:Eu³⁺ (50mg/ml), which were placed inside the cylinder as the targets. The tank was filled with a mixture of water and intralipid and the optical properties were set as $μ_{a} = 0.03 c m^{- 1}$ and $μ_{s}^{1} = 10 c m^{- 1}$ [29,30]. To evaluate the reconstruction results of two targets with different distances, numerical simulations were performed with the targets positioned at distances of 3, 2 and 1mm (edge-to-edge distance between the two targets), respectively.

Fig. 2 The cylinder phantom used in simulation studies. (a), (b), (c) A 3D view of the phantom, (d), (e), (f) the overhead view of the phantom. edge-to-edge distance between the two targets: (a),(d):3mm, (b),(e):2mm, (c),(f):1mm.

Download Full Size | PDF

For phantom simulations, the imaging model is discretized into 2,695 nodes and 12,285 tetrahedral elements in a 3D region of 3.0 × 3.0 × 2.3 cm³. To make the results comparable with phantom experiments, in the numerical simulations, the distance from the X-ray source to the rotation center of the imaging system was set as 26.3cm and the EMCCD camera was positioned perpendicularly to the X-ray source-detector axis, with a distance of 45.0cm between the EMCCD and the rotation center. The voltage and current of the cone beam X-ray source were set as 50kVp and 1mA, respectively. The simulated 24 projections were obtained every 15° during a 360° scan. After optical luminescence was obtained at different angles, white Gaussian noise was added into all projections with a zero-mean and signal-to-noise ratio (SNR) of 30DB to simulate noisy measurements.

3.2 Phantom experiments setup

In order to further validate the performance of the proposed method with real luminescence measurements, phantom experiments were performed by using a custom-developed CB-XLCT system, which was shown in Fig. 3. The system mainly included a micro-focus X-ray source (Oxford Instrument, U.K.), a rotation stage, a flat-panel X-ray detector (2923, Dexela, U.K.) for high-resolution CT imaging, and an electron-multiplying charge-coupled device (EMCCD) camera (iXon DU-897, Andor, U.K.) for optical imaging. The maximum voltage of the X-ray source was 80kVp with the maximum power of 80W. The distances from the X-ray source to the rotation center of the imaging system and to the flat-panel detector were 26.3cm and 86.3cm, respectively. The EMCCD camera coupled with a Nikon 50-mm f/1.8D lens was positioned at 90° towards the X-ray source- detector axis, with a distance of 45.0cm to the rotation center. The minimum cooling temperature of EMCCD camera was −80°, which could effectively reduce the dark noise.

Fig. 3 The schematic diagram of the CB-XLCT system.

Download Full Size | PDF

The configuration of the physical phantom used in imaging experiments was shown in Fig. 4, which was based on observations from the simulation studies. A glass cylinder (4.0cm in diameter, 4.0cm in height) containing a mixture of water and intralipid was fixed on the rotation stage. Two small glass tubes (3mm in diameter) filled with Y₂O₃: Eu³⁺ (50mg/ml) were symmetrically placed in the cylinder to simulate two targets. The edge-to-edge distances (EED) between the two tubes were 5.5mm and 1.8mm.

Fig. 4 The physical phantom used in imaging experiments. (a), (c) X-ray projections of the phantom corresponding to the two tubes with edge-to-edge distances of 5.5mm, 1.8mm. The region between two red lines is used for reconstruction in this study. (b), (d) CT slices of the phantom, corresponding to the transverse slices indicated by the blue lines in (a), (c).

Download Full Size | PDF

During imaging experiments, the voltage and current of the X-ray source were set as 50kVp and 1mA, respectively. The phantom was rotated from 0° to 360° and the optical images were obtained every 15° by the EMCCD camera. The exposure time of the EMCCD camera was set as 1s, with the EM gain and binning set as 260 and 1 × 1. For the CT imaging, the projections were acquired with an angular increment of 1° on a 360° circular orbit, while the voltage and current of the X-ray source were set as 50kVp and 1mA, and the acquisition time for each projection was 150ms. The Feldkamp-Davis-Kress (FDK) algorithm [31,32] was used for CT reconstruction.

3.3 Quantitative evaluation

The quality of reconstructed CB-XLCT images was evaluated quantitatively by several indexes including the location error (LE), dice similarity coefficient (DICE) and contrast-to-noise ratio (CNR) [30].

LE evaluates the localization accuracy of the reconstructed target, which is defined as the Euclidean distance error between the centers of true and reconstructed targets:

L E = {‖ L_{r} - L_{t} ‖}_{2}

where L_r and L_t denote the centers of the reconstructed and true targets, respectively.

DICE reflects the similarity of the true and reconstructed targets and can be calculated by:

D I C E = \frac{2 | R O I_{r} \cap R O I_{t} |}{| R O I_{r} | + | R O I_{t} |}

where ROI_t and ROI_r denote the regions of true and reconstructed targets, respectively, and |⋅| defines the number of voxels in a region.

CNR is used for quantitative evaluation of noise and artifacts in reconstructed images, as shown below:

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

where ROI and BCK denote the target and background regions of the imaged object,

w_{R O I}

and

w_{B C K}

are weighting factors determined by the relative volumes of the target and background,

μ_{R O I}

and

μ_{B C K}

are the mean intensity values of the ROI and BCK, and

σ_{R O I}^{2}

and

σ_{B C K}^{2}

represent the variances of the ROI and BCK, respectively.

4. Results

4.1 Numerical simulations

Firstly, the XLCT tomographic images of the targets positioned at different distances were reconstructed using 24 projections with different algorithms, as shown in Fig. 5. All the reconstruction results are normalized based on their maximum values. It can be seen for ART and ADAPTIK algorithms, the distribution of two targets can be effectively distinguished when the EED is 3 mm, but the task became difficult when the EED is 2 mm and 1 mm, with relatively large location errors, as shown in Figs. 5(d)-(i). Comparatively, the distribution of two targets can be effectively distinguished with FISTA and the L₁-TV methods, when the EED is 2 mm. However, when the EED is 1 mm, no algorithms except for the proposed L₁-TV algorithm can distinguish them, as shown in Figs. 5(m)-(o).

Fig. 5 Tomographic images of the targets positioned at different distance were reconstructed from 24 projections. 1^st row (from left to right): the true locations of the target edge to edge distances of 3mm, 2mm and 1mm respectively. 2^nd row (from left to right): reconstructions with the ART algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively. 3^rd row (from left to right): reconstructions with the ADAPTIK algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively. 4^th row (from left to right): reconstructions with the FISTA algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively 0.5^th row (from left to right): reconstructions with the proposed L₁-TV algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively.

Download Full Size | PDF

To evaluate the performance of the proposed method with fewer projections, the XLCT images of the targets positioned at different distance were reconstructed using 4 projections, as shown in Fig. 6. It can be seen that based on ART and ADAPTIK algorithms, the distribution of two targets can be effectively distinguished when the EED is 3 mm, but it is difficult to be distinguished when the EED is 2 mm and 1 mm using 4 projections and compared with the reconstruction results using 24 projections, the reconstruction results using 4 projections yield worse target shape and larger location error, as shown in Figs. 6(d)-(i). Besides, there is background noise in the reconstruction results because the ill-conditioned of reconstruction is more serious with the decrease of the projection angle. Correspondingly, the distribution of two targets can be effectively distinguished when the EED is 2 mm, but it is difficult to be distinguished when the EED is 1 mm using 4 projections based on the FISTA method, as shown in Figs. 6(j)-(l). Compared with the method of ART, ADAPTIK and FISTA, the distribution of two targets can be effectively distinguished when the EED is 1 mm based on the proposed L₁-TV algorithm although only 4 projections can be used, as shown in Figs. 6(m)-(o).

Fig. 6 Tomographic images of the targets positioned at different distance were reconstructed from 4 projections. 1^st row (from left to right): the true locations of the target edge to edge distances of 3mm, 2mm and 1mm respectively. 2^nd row (from left to right): reconstructions with the ART algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively. 3^rd row (from left to right): reconstructions with the ADAPTIK algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively. 4^th row (from left to right): reconstructions with the FISTA algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively 0.5^th row (from left to right): reconstructions with the proposed L₁-TV algorithm for the two targets with edge to edge distances of 3 mm, 2mm and 1mm respectively.

Download Full Size | PDF

Table 1 summarizes the quantitative evaluation of the reconstructions using different methods. Here LE1 and LE2 represent for the localization error of the reconstructed target 1 and target 2, respectively. For considering the prior information of the targets sparsity and local smoothness, the reconstruction results based on the proposed L₁-TV method yield highest reconstruction quality in terms of the target shape and localization accuracy, among the four methods, which further confirm the observation in Fig. 6.

Table 1. Quantitative evaluation of numerical simulations with 4 projections

View Table | View all tables in this article

4.2 Phantom experiments

In order to compare the performance of different algorithms, the XLCT tomographic images of the targets positioned at different distance for phantom experiments were reconstructed using 24 projections, as shown in Fig. 7. It can be seen that for ART and ADAPTIK algorithms, the distribution of two targets can be effectively distinguished when the EED is 5.5 mm, but the task became difficult when the EED is 1.8 mm, with relatively large location errors, as shown in Figs. 7 (c)-(k). Compared with the method of ART and ADAPTIK, FISTA method can reduce background noise. However, it is also difficult to distinguish the distribution of the two targets when the EED is 1.8mm, as shown in Figs. 7(l)-(o). Comparatively, when the EED is 1.8 mm, no algorithms except for the proposed L₁-TV algorithm can distinguish them, as shown in Figs. 7(p)-(s).

Fig. 7 Tomographic images of the targets positioned at different distance were reconstructed from 24 projections for phantom experiments. 1^st row (from left to right): the true locations of the target edge to edge distances of 5.5mm, 1.8mm respectively. 2nd row (from left to right): reconstructions with the ART algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 3^nd row (from left to right): reconstructions with the ADAPTIK algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 4^nd row (from left to right): reconstructions with the FISTA algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 5^nd row (from left to right): reconstructions with the proposed L₁-TV algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively.

Download Full Size | PDF

In order to further validate the performance of the proposed method for the fast imaging, the XLCT tomographic images of the targets positioned at different distance were reconstructed using 4 projection angles, as shown in Fig. 8. It can be seen that based on ART, it is difficult to be distinguished when the EED is 5.5 and 1.8 mm using 4 projection angles, as shown in Figs. 8(c)-(f). Based on the ADAPTIK method, the distribution of two targets can be effectively distinguished when the EED is 5.5 mm, but it is difficult to be distinguished when the EED is 1.8 mm using 4 projection angles. With sparse regularization, the distribution of two targets can be effectively distinguished when the EED is 5.5 mm, but it is difficult to be distinguished when the EED is 1.8 mm using 4 projection angles based on the FISTA method, as shown in Figs. 8(l)-(o). Figures 8(p)-(s) show that with our proposed L₁-TV algorithm, however, even the targets with a distance of 1.8 mm could be separated successfully, demonstrating its advantage in improving spatial resolution in real imaging experiments.

Fig. 8 Tomographic images of the targets positioned at different distance were reconstructed from 4 projections for phantom experiments. 1^st row (from left to right): the true locations of the target edge to edge distances of 5.5mm, 1.8mm respectively. 2^nd row (from left to right): reconstructions with the ART algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 3^nd row (from left to right): reconstructions with the ADAPTIK algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 4^nd row (from left to right): reconstructions with the FISTA algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively. 5^nd row (from left to right): reconstructions with the proposed L₁-TV algorithm and fusion of XLCT and XCT images for the two targets with edge to edge distances of 5.5mm, 1.8mm respectively.

Download Full Size | PDF

Table 2 summarizes the quantitative evaluation of the reconstructions using different methods with 4 projections. The results indicate that even in phantom experiments, the proposed L₁-TV algorithm still performs better in target location, shape recovery and image contrast, when compared to the conventional reconstruction methods, which further confirm the observation in Fig. 8.

Table 2. Quantitative evaluation of phantom experiments with 4 projections

View Table | View all tables in this article

CB-XLCT imaging using 4 projections requires less data acquisition time than that using 24 projections. In addition, due to reduced measured data, the dimension of the system matrix is reduced greatly, which is 14237 *8810 for 24 projections and 14237*1469 for four projections. Therefore, Compared with XLCT imaging with 24 projections, reconstructions with 4 projections could implement fast imaging with reduced scanning time and memory footprint, as shown in Table 3.

Table 3. The data acquisition and reconstruction time of phantom experiments with different projections

View Table | View all tables in this article

5. Discussion and conclusions

In this study, a reconstruction approach based on joint L₁ and total variation regularization is proposed for the CB-XLCT inverse problem. The implementation of this method is split into the L₁ penalty and the TV penalty. The L₁ constraint is used to enhance the target sparsity and the TV constraint is used to maintain the local smoothness and preserve the edges of the targets of the reconstructed image.

Numerical simulations and phantom experiments results confirm the superiority of the proposed method over the conventional ART, ADAPTIK and the FISTA methods. In numerical simulations, two targets could be effectively distinguished by the proposed method when the EED is 1 mm, while in phantom experiments, two targets could be resolved when the EED is 1.8 mm, demonstrating its advantage in improving spatial resolution. With the consideration of the prior information on target sparsity and local smoothness, the proposed method indicates its potential in improving the XLCT image quality in terms of target shape and localization accuracy, and in reducing imaging time by reconstructions with fewer scanning angles.

For the phantom experiments, the whole scanning time consist of the rotation time and the acquisition time of the EMCCD camera. In our imaging experiments, the speed of the rotation was 6 degrees per second and the exposure time of the EMCCD camera was set as 1s per view. That resulted in more time spent on the rotation process, where the acquisition time of 4 projections was quite limited. The current configuration of the imaging system makes the reduction of imaging time using the proposed method not obvious in phantom experiments. With the increase of rotation speed and the exposure time in in-vivo experiments, the benefit of the proposed method on imaging time reduction would be further recognized.

Although the system setup for numerical simulations and phantom experiments appeared similar, the results were not identical, especially when the target distance was 1.8 mm in phantom experiments. The possible reason may be that in numerical simulations, inaccuracies were mainly caused by simulation noise and numerical errors due to the ill-posedness of the system matrix. In contrast, in phantom experiments, inaccurate modeling of optical properties and measurement noise, geometric error of the imaging system, and other factors, would further deteriorate the reconstruction results.

In this study, the fast iterative shrinkage-thresholding (FISTA) algorithm is used to solve L₁ regular constraint. Besides FISTA algorithm, it can be solved based on the L₁-Ls algorithm [33] and the separable paraboloidal surrogates algorithm [23,24]. In addition, the traditional TV model is adopted in the proposed L₁-TV method. To further improve the performance of the proposed method, the edge-preserving total variation (EPTV) [34] and the adaptive-weighted total variation (AWTV) [35] can be used in the future studies. Nevertheless, due to the uncertainty of the CB-XLCT forward model, the reconstruction is very ill-conditioned. It is difficult to obtain high quality reconstructions in the phantom experiments, especially when targets are closer. In addition to studies on reconstruction methods, further investigations on constructing more accurate propagation model of photons, such as by using the radiation transfer equation, may provide additional benefit for CB-XLCT imaging.

In summary, a reconstruction approach based on joint L₁ and total variation regularization was proposed in this study. Compared to the traditional ART, ADAPTIK and FISTA method, the proposed method demonstrated its advantage in improving spatial resolution and reducing imaging time, which can promote the widely use of CB-XLCT in vivo.

Funding

National Key Research and Development Program of China (2017YFC0107400); National Natural Science Foundation of China (NSFC) (81230035, 11805274); Natural Science Foundation of Shaanxi Province (2016JQ1012); Key project supported by Military Science and Technology Foundation (BWS14C030); NIH (CA206171).

Disclosures

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

References

1. 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]

2. 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]

3. J. Feng, C. Qin, K. Jia, S. Zhu, X. Yang, and J. Tian, “Bioluminescence tomography imaging in vivo: recent advances,” IEEE J Sel Top Quant 18(4), 1394–1402 (2012). [CrossRef]

4. C. Darne, Y. Lu, and E. M. Sevick-Muraca, “Small animal fluorescence and bioluminescence tomography: a review of approaches, algorithms and technology update,” Phys. Med. Biol. 59(1), R1–R64 (2014). [CrossRef] [PubMed]

5. Y. Zhu, A. K. Jha, J. K. Dreyer, H. N. Le, J. U. Kang, P. E. Roland, D. F. Wong, and A. Rahmim, “A three-step reconstruction method for fluorescence molecular tomography based on compressive sensing,” SPIE Bios. Proc SPIE Int Soc Opt Eng, 1005911 (2017).

6. R. Baikejiang, Y. Zhao, B. Z. Fite, K. W. Ferrara, and C. Li, “Anatomical image-guided fluorescence molecular tomography reconstruction using kernel method,” J. Biomed. Opt. 22(5), 55001 (2017). [CrossRef] [PubMed]

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

8. 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]

9. 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]

10. 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]

11. B. K. Natarajan, “Sparse approximate solutions to linear systems,” SIAM J. Comput. 24(2), 227–234 (1995). [CrossRef]

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

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

14. M. A. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE J. Sel. Top. Signal Process. 1(4), 586–597 (2007). [CrossRef]

15. I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Commun. Pure Appl. Math. 57(11), 1413–1457 (2004). [CrossRef]

16. 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]

17. A. Y. Yang, S. S. Sastry, A. Ganesh, and Y. Ma, “Fast ℓ 1-minimization algorithms and an application in robust face recognition: A review,” Image Processing (ICIP), 17th IEEE International Conference on. IEEE, 1849–1852 (2010).

18. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm with application to wavelet-based image deblurring,” in Acoustics, Speech and Signal Processing, 2009. ICASSP 2009. IEEE International Conference on, 693–696 (2009). [CrossRef]

19. T. Hebert and R. Leahy, “A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors,” IEEE Trans. Med. Imaging 8(2), 194–202 (1989). [CrossRef] [PubMed]

20. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D 60(1-4), 259–268 (1992). [CrossRef]

21. E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” J. XRay Sci. Technol. 14, 119–139 (2006).

22. E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53(17), 4777–4807 (2008). [CrossRef] [PubMed]

23. M. Defrise, C. Vanhove, and X. Liu, “An algorithm for total variation regularization in high-dimensional linear problems,” Inverse Probl. 27(6), 065002 (2011). [CrossRef]

24. J. Dutta, S. Ahn, C. Li, S. R. Cherry, and R. M. Leahy, “Joint L1 and total variation regularization for fluorescence molecular tomography,” Phys. Med. Biol. 57(6), 1459–1476 (2012). [CrossRef] [PubMed]

25. 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]

26. 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]

27. M. Schweiger, S. R. Arridge, M. Hiraoka, and D. T. Delpy, “The finite element method for the propagation of light in scattering media: boundary and source conditions,” Med. Phys. 22(11 Pt 1), 1779–1792 (1995). [CrossRef] [PubMed]

28. Y. Lv, J. Tian, W. Cong, G. Wang, J. Luo, W. Yang, and H. Li, “A multilevel adaptive finite element algorithm for bioluminescence tomography,” Opt. Express 14(18), 8211–8223 (2006). [CrossRef] [PubMed]

29. 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]

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

31. D. Dong, S. Zhu, C. Qin, V. Kumar, J. V. Stein, S. Oehler, C. Savakis, J. Tian, and J. Ripoll, “Automated recovery of the center of rotation in optical projection tomography in the presence of scattering,” IEEE J. Biomed. Health Inform. 17(1), 198–204 (2013). [CrossRef] [PubMed]

32. S. Zhu, J. Tian, G. Yan, C. Qin, and J. Feng, “Cone beam micro-CT system for small animal imaging and performance evaluation,” Int. J. Biomed. Imaging 2009, 960573 (2009). [CrossRef] [PubMed]

33. H. Gao and H. Zhao, “Multilevel bioluminescence tomography based on radiative transfer equation Part 1: l1 regularization,” Opt. Express 18(3), 1854–1871 (2010). [CrossRef] [PubMed]

34. 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]

35. Y. Liu, J. Ma, Y. Fan, and Z. Liang, “Adaptive-weighted total variation minimization for sparse data toward low-dose x-ray computed tomography image reconstruction,” Phys. Med. Biol. 57(23), 7923–7956 (2012). [CrossRef] [PubMed]

Targets distances /mm	Indexes	ART	ADAPTIK	FISTA	L₁-TV
3	LE1(mm)	0.56	0.63	0.59	0.34
	LE2(mm)	0.84	0.81	0.67	0.35
	Dice	0.46	0.39	0.58	0.61
	CNR	7.7	6.9	10.2	11.7
2	LE1(mm)	0.74	0.81	0.49	0.31
	LE2(mm)	0.62	0.65	0.59	0.41
	Dice	0.32	0.18	0.55	0.63
	CNR	7.9	6.3	9.2	10.8
1	LE1(mm)	0.63	0.8	0.79	0.54
	LE2(mm)	1	1.2	1.1	0.75
	Dice	0.37	0.19	0.41	0.67
	CNR	7.6	6.8	8.4	9.1

Targets distances /mm	Indexes	ART	ADAPTIK	FISTA	L₁-TV
5.5	LE1(mm)	3.6	2.8	2.6	2.2
	LE2(mm)	2.1	1.9	1.1	0.78
	Dice	0.22	0.35	0.4	0.46
	CNR	1.4	2.1	2.5	3.7
1.8	LE1(mm)	4.8	4.5	2.8	2.3
	LE2(mm)	2.9	2.3	0.93	0.67
	Dice	0.29	0.32	0.41	0.45
	CNR	1.6	1.7	2.5	3.1

Projections	data acquisition time /s	reconstruction time /s
24	84	343
4	49	298

Targets distances /mm	Indexes	ART	ADAPTIK	FISTA	L₁-TV
3	LE1(mm)	0.56	0.63	0.59	0.34
	LE2(mm)	0.84	0.81	0.67	0.35
	Dice	0.46	0.39	0.58	0.61
	CNR	7.7	6.9	10.2	11.7
2	LE1(mm)	0.74	0.81	0.49	0.31
	LE2(mm)	0.62	0.65	0.59	0.41
	Dice	0.32	0.18	0.55	0.63
	CNR	7.9	6.3	9.2	10.8
1	LE1(mm)	0.63	0.8	0.79	0.54
	LE2(mm)	1	1.2	1.1	0.75
	Dice	0.37	0.19	0.41	0.67
	CNR	7.6	6.8	8.4	9.1

Targets distances /mm	Indexes	ART	ADAPTIK	FISTA	L₁-TV
5.5	LE1(mm)	3.6	2.8	2.6	2.2
	LE2(mm)	2.1	1.9	1.1	0.78
	Dice	0.22	0.35	0.4	0.46
	CNR	1.4	2.1	2.5	3.7
1.8	LE1(mm)	4.8	4.5	2.8	2.3
	LE2(mm)	2.9	2.3	0.93	0.67
	Dice	0.29	0.32	0.41	0.45
	CNR	1.6	1.7	2.5	3.1

Regularized reconstruction based on joint L₁ and total variation for sparse-view cone-beam X-ray luminescence computed tomography

Abstract

1. Introduction