A CT reconstruction algorithm based on non-aliasing contourlet transform and compressive sensing.
Since computed tomography (CT)  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)  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  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  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) . 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  method, which is derived from Bregman iteration and can accelerate iteration convergence and produce better reconstruction results. Split-Bregman method uses an intermediate variable to split [L.sub.1] regularization and [L.sub.2] regularization into two equations, where [L.sub.2] and [L.sub.1] 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 . Furthermore, an iterative CT reconstruction is proposed using shearlet-based regularization . Chu et al. proposed multienergy CT reconstruction based on low rank and sparsity with the Split-Bregman method (MLRSS) . Chang et al. proposed a few-view reweighted sparsity hunting (FRESH) method for CT images reconstruction .
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:
Af = p, (1)
where A is the system matrix, f is the original image, and p is the projection data. Traditional CT reconstruction algorithms such as filtered backprojection (FBP)  and algebraic reconstruction technique (ART)  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:
min [[parallel]y[parallel].sub.0] s.t. p = Af = A[[PHI].sup.H]y, (2)
where [PHI] is a orthogonal transform, [[PHI].sup.H] is the corresponding inverse transform, and f is the CT image to be reconstructed and has a special relationship with [[PHI].sup.H]; that is, f = [[PHI].sup.H]y. p is the projection data of f through matrix A.
Inspired by CS theory, Sidky et al. proposed a total variation (TV-) based CT reconstruction algorithm using gradient operator as the sparse representation , in which TV is defined as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)
where [nabla]f represents gradient operator of an image f.
2.2. Non-Aliasing Contourlet Transform. The traditional sparse representation, such as gradient operator and Wavelet transform , cannot get ideal sparse representation of CT images. In 2002, Ni et al. proposed Contourlet transform  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: [L.sub.0]([omega]), [D.sub.0]([omega]) and [L.sub.1]([omega]), [D.sub.1]([omega]). [L.sub.0]([omega]) and [L.sub.1]([omega]) mean low-pass filters. [D.sub.0]([omega]) and [D.sub.1]([omega]) mean high-pass filters. The relationships of two different filter banks are as follows:
[D.sup.2.sub.1]([omega]) + [L.sup.2.sub.1]([omega]) = 1 [D.sup.2.sub.0]([omega]) + [L.sup.2.sub.0]([omega]) = 1. (4)
We assume that [[omega].sub.p,0] and [[omega].sub.s,0] represent pass-band frequency and stop-band frequency of [L.sub.0]([omega]), respectively. Accordingly, [[omega].sub.p,1] and [[omega].sub.s,1] represent pass-band frequency and stop-band frequency of [L.sub.1]([omega]), respectively. In order to eliminate frequency aliasing, the filter parameters should meet (1) [[omega].sub.s,1] < [pi]/2; (2) [[omega].sub.p,0] + [[omega].sub.s,0])/2 = [pi]/2 and ([[omega].sub.p,1] + [[omega].sub.s,1])/2 = [pi]/4; (3) [[omega].sub.s,0] [less than or equal to] [pi] - a and [[omega].sub.s,1] [less than or equal to] ([pi] - a)/2, where a is the maximum width of mixing ingredients in DFB .
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  in the study.
2.3. Split Bregman Method. In CS theory, [L.sub.0] norm is the most ideal regularization norm, but it is difficult to solve equations and easily interfered by noise in CT reconstruction, so [L.sub.0] norm is commonly replaced by [L.sub.1] norm. Then the reconstruction problem depicted by (2) can be converted into
min [[parallel]y[parallel].sub.0] s.t. p = Af = A[[PHI].sup.H]y (5a)
min [[parallel][PHI]f[parallel].sub.1] s.t. p = Af = A[[PHI].sup.H]y, (5b)
where y = [PHI]f, [PHI] is the sparse transform which is normally used as Wavelet transform, Curvelet transform, gradient operator, and so forth.
Furthermore, (5b) can be converted into
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (6)
where [lambda] is penalty function weight.
In order to solve (6), Goldstein and Osher proposed SplitBregman method , using an intermediate variable to split [L.sub.1] regularization and [L.sub.2] regularization into two equations; [L.sub.2] regularization equation can be solved by gradient descent method and [L.sub.1] regularization equation can be solved by thresholding algorithm. Split-Bregman method contains the following three iteration steps:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (7)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (8)
[b.sup.k+1] = [b.sup.k] + ([PHI][f.sup.k+1] - [d.sup.k+1]), (9)
where k is the Split-Bregman iteration index, p is convergence parameter, and d and b 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:
min [[parallel]f[parallel].sub.TV] + [[parallel][PHI]f[parallel].sub.1]. s.t. [[parallel]Af - p[parallel].sup.2.sub.2] < [[sigma].sup.2] (10)
min [[parallel]f[parallel].sub.TV] + [[parallel][PHI]f[parallel].sub.1] + [lambda] [[parallel]Af - p[parallel].sup.2.sub.2]. (11)
Applying the Split-Bregman method to (11), we have the following three iteration steps:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (12)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (13)
[b.sup.k+1] = [b.sup.k] + ([nabla][f.sup.k+1] - [d.sup.k+1]) [b.sup.k+1.sub.[phi]] = [b.sup.k.sub.[phi]] + ([PHI][f.sup.k+1] - [d.sup.k+1.sub.[phi]]), (14)
where [gamma] is convergence parameter and [d.sub.[phi]] and [b.sub.[phi]] are intermediate variables.
The steepest descent method is applied to solve (12). The derivative of (12) is calculated as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (15)
where n denotes the iteration index of the steepest descent method, m = 2, ..., [N.sub.data] denotes the projection angles, [A.sub.m] is mth row vector, and system matrix A includes [N.sub.data] row vector [A.sub.m]. Accordingly, [N.sub.data] row vectors [p.sub.m] compose the projection-data vector p, [alpha] 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:
[d.sup.k+1] = shrink ([nabla][f.sup.k+1] + [b.sup.k], [1/[lambda]]) [d.sup.k+1.sub.[phi]] = shrink ([PHI][f.sup.k+1] + [b.sup.k.sub.[phi]], [1/[mu]]). (16)
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 n and the inside loop is labeled by k. 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 y<ART-DATA) [n, m] to denote the mth DATA-step subiteration with the nth iteration and [y.sup.(ART-POS)] [n] to denote the POS-step with the nth iteration in the outside loop. We use [f.sup.(NACTTV-DATA)][k, m] to denote the mth DATA-step subiteration with the kth iteration and [f.sup.(NACTTV-POS)][k] to denote the POS-step with the kth iteration in the inside loop. The steps of the algorithm are as follows:
n = 1, [f.sup.(ART-DATA)][n, 1] = 0; (17)
(B) data projection iteration, for m = 2, ..., [N.sub.data]:
[f.sub.(ART-DATA)][n, m] = [f.sub.(ART-DATA)][n, m-1] + [A.sub.m] [[[[p.sub.m] - [A.sub.m] x [f.sub.(ART-DATA)][n, m-1]]/[[A.sub.m] x [A.sub.m]]]; (18)
(C) positivity constraint:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (19)
(D) initialization of Split-Bregman:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (20)
(E) iteration for m = 2, ..., [N.sub.data]:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (21)
(F) positivity constraint:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (22)
(G) update [d.sub.x], [d.sub.y], [d.sub.[phi]], [b.sub.x], [b.sub.y], [b.sub.[phi]], increase k, and return to step (E) until k = [K.sub.NACTTV] as follows:
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]; (23)
(H) initialize next loop:
[f.sup.(ART-DATA)][n + 1, 1] = [f.sup.(NACTTV-POS)] [K.sub.NACTTV], 1]; (24)
increase n and return to step (B). The iteration is stopped when [[parallel]Af - p[parallel].sup.2.sub.2] < [[sigma].sup.2]. In our study, we selected [lambda] = 1000, [gamma] = 30, [mu] = 30, a = 0.2 and [K.sub.NACTTV] = 10, 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)  to evaluate the quality of the reconstructed images.
RMSE is the most widely applied way to evaluate image quality, and RMSE is defined as
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (25)
where [f.sub.i,j] is the pixel value of original image and [f.sup.R.sub.i,j] 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
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (26)
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (27)
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 200 x 200. We assume that the CT system was viewed as in a typical pencil-beam geometry, and the scanning range was from 1[degrees] to 360[degrees] with a [theta] angular increment; projection angles can be indicated as
[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (28)
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 [N.sub.view] is 60, 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.
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 [N.sub.view]. 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 [N.sub.view] 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.
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.
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).
 G. Wang, H. Yu, and B. de Man, "An outlook on X-ray CT research and development," Medical Physics, vol. 35, no. 3, pp. 1051-1064, 2008.
 D. L. Donoho, "Compressed sensing," IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289-1306, 2006.
 M. N. Do and M. Vetterli, "The contourlet transform: an efficient directional multiresolution image representation," IEEE Transactions on Image Processing, vol. 14, no. 12, pp. 2091-2106, 2005.
 Y. Lu and M. N. Do, "A new contourlet transform with sharp frequency localization," in Proceedings of the IEEE International Conference on Image Processing (ICIP 06), vol. 2, pp. 1629-1632, Atlanta, Ga, USA, October 2006.
 P. Feng, B. Wei, Y. J. Pan, and D. L. Mi, "Construction of non-aliasing pyramidal transform," Acta Electronica Sinica, vol. 37, no. 11, pp. 2510-2514, 2009.
 T Goldstein and S. Osher, "The split Bregman method for L1-regularized problems," SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 323-343, 2009.
 L. Bregman, "The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex optimization," USSR Computational Mathematics and Mathematical Physics, vol. 7, pp. 200-217, 1967.
 B. Vandeghinste, B. Goossens, J. de Beenhouwer et al., "Split-Bregman-based sparse-view CT reconstruction," in Proceedings of the 11th International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine (Fully 3D '11), pp. 431-434, 2011.
 B. Vandeghinste, B. Goossens, R. van Holen et al., "Iterative CT reconstruction using shearlet-based regularization," in Medical Imaging 2012: Physics of Medical Imaging, vol. 8313 of Proceedings ofSPIE, p. 83133I, San Diego, Calif, USA, February 2012.
 J. Chu, L. Li, Z. Chen, G. Wang, and H. Gao, "Multi-energy CT reconstruction based on low rank and sparsity with the split-bregman method (MLRSS)," in Proceedings of the IEEE Nuclear Science Symposium and Medical Imaging Conference Record (NSS/MIC '12), pp. 2411-2414, Anaheim, Calif, USA, November 2012.
 M. Chang, L. Li, Z. Chen, Y. Xiao, L. Zhang, and G. Wang, "A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction," Journal of X-Ray Science and Technology, vol. 21, no. 2, pp. 161-176, 2013.
 S. L. Zhang, W. B. Li, and G. F. Tang, "Study on image reconstruction algorithm of filtered backprojection," Journal of Xianyang Normal University, vol. 23, no. 4, pp. 47-49, 2008.
 R. Gordon, R. Bender, and G. T Herman, "Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography," Journal of Theoretical Biology, vol. 29, no. 3, pp. 471-481, 1970.
 E. Y. Sidky, C. Kao, and X. Pan, "Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT," Journal of X-Ray Science and Technology, vol. 14, no. 2, pp. 119-139, 2006.
 M. Abramowitz and C. A. Stegun, A Wavelet Tour of Signal Processing, Academic Press, San Diego, Calif, USA, 3rd edition, 2008.
 X. Ni, H. L. Wang, L. Chen, and J.M. Wang, "Image compressed sensing based on sparse representation using Contourlet directional subbands," Application Research of Computers, vol. 30,no. 6, pp. 1889-1898, 2013.
 S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
 Y. R. Liu, Research on Objective Full-Reference Image Quality Evaluation Method, Computer Science & Technology, Nanjing, China, 2010.
Lu-zhen Deng, Peng Feng, Mian-yi Chen, Peng He, Quang-sang Vo, and Biao Wei
The Key Lab of Optoelectronic Technology and Systems of the Education Ministry of China, Chongqing University, Chongqing 400044, China
Correspondence should be addressed to Peng Feng; firstname.lastname@example.org and Peng He; email@example.com
Received 11 April 2014; Revised 6 June 2014; Accepted 9 June 2014; Published 30 June 2014
Academic Editor: Fenglin Liu
Table 1: RMSE and UQI of reconstructed images using three different algorithms. RMSE UQI Methods art art-tv SpBr-NACT art art-tv SpBr-NACT Noisy-free 0.0502 0.0321 0.0196 0.9869 0.9947 0.9980 Noisy 0.0554 0.0406 0.0318 0.9839 0.9914 0.9948
|Printer friendly Cite/link Email Feedback|
|Title Annotation:||Research Article|
|Author:||Deng, Lu-zhen; Feng, Peng; Chen, Mian-yi; He, Peng; Vo, Quang-sang; Wei, Biao|
|Publication:||Computational and Mathematical Methods in Medicine|
|Date:||Jan 1, 2014|
|Previous Article:||A new multistage medical segmentation method based on superpixel and fuzzy clustering.|
|Next Article:||Ultrasound image enhancement using structure-based filtering.|