Abstract
Compressive sensing (CS) theory has great potential for reconstructing CT images from sparse-views projection data. Currently, total variation (TV-) based CT reconstruction method is a hot research point in medical CT field, which uses the gradient operator as the sparse representation approach during the iteration process. However, the images reconstructed by this method often suffer the smoothing problem; to improve the quality of reconstructed images, this paper proposed a hybrid reconstruction method combining TV and non-aliasing Contourlet transform (NACT) and using the Split-Bregman method to solve the optimization problem. Finally, the simulation results show that the proposed algorithm can reconstruct high-quality CT images from few-views projection using less iteration numbers, which is more effective in suppressing noise and artefacts than algebraic reconstruction technique (ART) and TV-based reconstruction method.
1. Introduction
Since computed tomography (CT) [1] technique was born in 1973, CT has been widely applied in medical diagnose, industrial nondestructive detection, and so forth. In medical CT field, how to reconstruct high-quality CT images from few-views or sparse-views data is a significant research problem. Recently, compressive sensing (CS) [2] theory has been applied in CT images reconstruction which makes it possible to reconstruct high-quality images from few-views data. In CS theory, CT images can be sparsely represented in an appropriate domain, such as gradient transform and Wavelet transform, and the quality of CT reconstructed images will be improved by some appropriate sparse representations in CT images reconstruction.
Contourlet transform [3] is proposed by Do and Vetterli in 2002, which is a sparse representation for 2D images with some properties such as multiresolution, multiscale, and multidirection. Contourlet transform can also get important smooth contour features of the image with few data, but there is frequency aliasing in Contourlet transform. Sharp frequency localization Contourlet transform [4] is firstly proposed by Lu and Do in 2006 and Feng et al. introduced a detailed explanation and construction in 2009 which is named as non-aliasing Contourlet transform (NACT) [5]. NACT which can eliminate the frequency aliasing in Contourlet transform is more efficient in capturing geometrical structure and can represent image sparser than traditional Contourlet transform.
To solve the optimization problem in CT images reconstruction based on CS, Goldstein and Osher proposed Split-Bregman [6] method, which is derived from Bregman [7] iteration and can accelerate iteration convergence and produce better reconstruction results. Split-Bregman method uses an intermediate variable to split regularization and regularization into two equations, where and regularization equation can be solved by steepest descent method and thresholding algorithm, respectively. Based on Split-Bregman method, Vandeghinste et al. proposed Split-Bregman-based sparse-view CT reconstruction approach [8]. Furthermore, an iterative CT reconstruction is proposed using shearlet-based regularization [9]. Chu et al. proposed multienergy CT reconstruction based on low rank and sparsity with the Split-Bregman method (MLRSS) [10]. Chang et al. proposed a few-view reweighted sparsity hunting (FRESH) method for CT images reconstruction [11].
In this paper, we propose a CT reconstruction algorithm based on NACT and compressive sensing which tries to explore the sparse capability of NACT in order to reconstruct high-quality CT images. In the following section, the proposed algorithm will be introduced. In the third section, we will analyze the experimental results and discuss relevant issues. In the last section, Section 4, we will conclude the paper.
2. Theory and Method
2.1. CT Reconstruction Theory Based on Compressive Sensing
Theoretically, the mathematical CT model can be expressed as: where is the system matrix, is the original image, and is the projection data. Traditional CT reconstruction algorithms such as filtered backprojection (FBP) [12] and algebraic reconstruction technique (ART) [13] cannot reconstruct high quality CT images with the sparse sampling or limited projection data.
In 2006, Candes and Donoho put forward the CS theory which makes it possible to get high quality CT images with sparse projection data. The main idea of CS is that a signal can be reconstructed with far less sampled frequency than required by conventional Nyquist-Shannon sampling frequency, if the image has a sparse/compressible representation in a transform domain.
Compressive sensing theory can be expressed by the following equation: where is a orthogonal transform, is the corresponding inverse transform, and is the CT image to be reconstructed and has a special relationship with ; that is, . is the projection data of through matrix .
Inspired by CS theory, Sidky et al. proposed a total variation (TV-) based CT reconstruction algorithm using gradient operator as the sparse representation [14], in which TV is defined as follows: where represents gradient operator of an image .
2.2. Non-Aliasing Contourlet Transform
The traditional sparse representation, such as gradient operator and Wavelet transform [15], cannot get ideal sparse representation of CT images. In 2002, Ni et al. proposed Contourlet transform [16] which can utilize intrinsic structure information of image to represent images more efficiently compared with Wavelet transform. However, suffering from frequency aliasing, Contourlet transform does not show good performance in image denoising, fusion, and enhancement. In order to solve this problem, a new multiscale analysis method, named non-aliasing Contourlet transform (NACT) was proposed. NACT consists of non-aliasing pyramidal filter banks (NPFB) and directional filter banks (DFB). NPFB contains two different filter banks: and . and mean low-pass filters. and mean high-pass filters. The relationships of two different filter banks are as follows:
We assume that and represent pass-band frequency and stop-band frequency of , respectively. Accordingly, and represent pass-band frequency and stop-band frequency of , respectively. In order to eliminate frequency aliasing, the filter parameters should meet ; and ; and , where is the maximum width of mixing ingredients in DFB [3].
As a sparse representation approach, NACT integrate NPFB and DFB which can decompose image into multidirection and multiresolution. NPFB decomposes image into an approximation subband and several detail subbands with different resolutions; DFB decomposes the detail subbands into directional subbands. The process of decomposition with 3 levels is shown in Figure 1. We will use “9-7” filter and “pkva” directional filter bank [17] in the study.
2.3. Split Bregman Method
In CS theory, norm is the most ideal regularization norm, but it is difficult to solve equations and easily interfered by noise in CT reconstruction, so norm is commonly replaced by norm. Then the reconstruction problem depicted by (2) can be converted intowhere , is the sparse transform which is normally used as Wavelet transform, Curvelet transform, gradient operator, and so forth.
Furthermore, (5b) can be converted into where is penalty function weight.
In order to solve (6), Goldstein and Osher proposed Split-Bregman method [6], using an intermediate variable to split regularization and regularization into two equations; regularization equation can be solved by gradient descent method and regularization equation can be solved by thresholding algorithm. Split-Bregman method contains the following three iteration steps:
Step 1.
Step 2.
Step 3. where is the Split-Bregman iteration index, is convergence parameter, and and are intermediate variables, with which each subproblem can be solved easily.
2.4. Proposed Algorithm
According to aforementioned methods, we propose a CT reconstruction algorithm based on NACT and compressive sensing method which can be defined as a constrained form (10) or an unconstrained form (11) as follows: Applying the Split-Bregman method to (11), we have the following three iteration steps:
Step 1.
Step 2.
Step 3. where is convergence parameter and and are intermediate variables.
The steepest descent method is applied to solve (12). The derivative of (12) is calculated as follows: where denotes the iteration index of the steepest descent method, denotes the projection angles, is mth row vector, and system matrix includes row vector . Accordingly, row vectors compose the projection-data vector , is an appropriate step size. The ART method is used to get initial image of iteration. Equation (13) can be explicitly computed as (16) using the shrinkage operator as follows:
We now describe the iterative steps of the proposed algorithm. The iteration process contains two loops, the outside loop operate ART and the inside loop solve the optimization problem which is constrained by TV and NACT. The outside loop is labeled by and the inside loop is labeled by . The steps comprising each loop are the DATA-step, which enforces consistency with the projection data; the POS-step, which ensures a nonnegative image. We use to denote the th DATA-step subiteration with the th iteration and to denote the POS-step with the th iteration in the outside loop. We use to denote the th DATA-step subiteration with the th iteration and to denote the POS-step with the th iteration in the inside loop. The steps of the algorithm are as follows:(A)initialization: (B)data projection iteration, for : (C)positivity constraint: (D)initialization of Split-Bregman: (E)iteration for : (F)positivity constraint: (G)update , , , , , , increase , and return to step (E) until as follows: (H)initialize next loop: increase and return to step . The iteration is stopped when . In our study, we selected , , , and , which can strike a good balance in the steepest descent and generate good reconstruction results in the experiments.
3. Experimental Results
3.1. The Image Quality Evaluation
This paper uses the root mean square errors (RMSE) and universal quality index (UQI) [18] to evaluate the quality of the reconstructed images.
RMSE is the most widely applied way to evaluate image quality, and RMSE is defined as where is the pixel value of original image and is the pixel value of reconstructed image.
Wang and Bovic proposed UQI mode which evaluates images distortion problem including correlation distortion, brightness distortion, and contrast distortion. The value of UQI is between −1 and 1. When the reconstructed image is the same as the original image, the value of UQI is 1. UQI is defined as where
3.2. Numerical Simulation
In this section, a head phantom as shown in Figure 2 is used to reconstruct and compare by 3 different methods: ART, ART-TV, and our proposed algorithm (SpBr-NACT). The size of phantom image is . We assume that the CT system was viewed as in a typical pencil-beam geometry, and the scanning range was from 1° to 360° with a angular increment; projection angles can be indicated as
In the simulation, we reconstruct the head phantom from noise-free and noisy projection data. To obtain noisy projection data, we add 10 dB Gaussian noise into noise-free projection data. Projection number is , and iteration numbers for all reconstruction algorithms are 50. The reconstructed images are shown in Figure 3 and the profile of line 140 in different reconstructed images is plotted in Figure 4.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
From Figures 3 and 4, we can see that the reconstructed images using ART and ART-TV methods contain a lot of noise and artifacts, while the reconstructed images using SpBr-NACT method contain less noise and artifacts and have clearer edges.
Table 1 lists all the RMSE and UQI calculated from reconstructed images. It is obvious that the RMSE of reconstructed images using SpBr-NACT method is much smaller than that of reconstructed images using ART and ART-TV methods; the UQI is much bigger. Thus SpBr-NACT method can reconstruct higher quality images.
Figure 5 plots the change of RMSE and UQI with respect to iteration number. Figure 6 plots the change of RMSE and UQI with respect to projection number . In both figures, ART, ART-TV, and the proposed SpBr-NACT approach are used to reconstruct images from noise-free and noisy data. The blue-solid line is for ART, the green-dashed line is for ART-TV, and the red dashed line is SpBr-NACT. For Figure 5, the projection number is fixed and is 60. For figure 6, the iteration number is fixed and equals 50. From both Figures, it is easy to find that with the increase of projection number or iteration number, SpBr-NACT approach can always get the minimum RMSE and maximum UQI which means that the quality of reconstructed images with SpBr-NACT is better than those with ART and ART-TV. And also we see from Figure 5 when the iteration number is relatively small, 3 methods that almost have the same RMSE and UQI, which implies that our proposed method has no advantage if the iteration step does not converge.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4. Conclusion
In this study, we proposed a CT reconstruction algorithm based on NACT and compressive sensing. The experimental results demonstrate that the proposed method can reconstruct high-quality images from few-views data and has a potential for reducing the radiation dose in clinical application. In the further research, we will try to explore more directional information from NACT so as to improve the performance of SpBr-NACT algorithm, especially when the projection number is far more below what we setup in the current experiment.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China (no. 61201346) and the Fundamental Research Funds for the Central Universities (no. 106112013CDJZR120020 and no. CDJZR14125501).