Sparse-View CT Image Recovery Using Two-Step Iterative Shrinkage-Thresholding Algorithm
Sparse-View CT Image Recovery Using Two-Step Iterative Shrinkage-Thresholding Algorithm
ETRI Journal. 2015. Dec, 37(6): 1251-1258
Copyright © 2015, Electronics and Telecommunications Research Institute (ETRI)
  • Received : May 07, 2015
  • Accepted : October 29, 2015
  • Published : December 01, 2015
Export by style
Cited by
About the Authors
Byung Gyu, Chae
Sooyeul, Lee

We investigate an image recovery method for sparse-view computed tomography (CT) using an iterative shrinkage algorithm based on a second-order approach. The two-step iterative shrinkage-thresholding (TwIST) algorithm including a total variation regularization technique is elucidated to be more robust than other first-order methods; it enables a perfect restoration of an original image even if given only a few projection views of a parallel-beam geometry. We find that the incoherency of a projection system matrix in CT geometry sufficiently satisfies the exact reconstruction principle even when the matrix itself has a large condition number. Image reconstruction from fan-beam CT can be well carried out, but the retrieval performance is very low when compared to a parallel-beam geometry. This is considered to be due to the matrix complexity of the projection geometry. We also evaluate the image retrieval performance of the TwIST algorithm using measured projection data.
I. Introduction
Image reconstruction technology for sparse-view computed tomography (CT) has become increasingly important in obtaining medical images without any harmful effects [1] [3] . Given a sufficient number of projection views, a high-quality CT image can be retrieved; however, human cell damage from too much radiation exposure is unavoidable, particularly in image-guided surgery and therapy [4] . Therefore, image restoration from less projection data is required to avoid radiation damage during the measurement process.
Compressed sensing (CS) technology [5] [7] has recently shown that tomographic images can be well retrieved from far less sample data than the Shannon-Nyquist criterion. It has been reported that an original image can be exactly reconstructed even when using under-sampled data of a Fourier space [5] ; this led to the development of a new reconstruction technique — a CT image retrieval algorithm based on a Radon space [8] . Various optimization algorithms on the basis of CS theory have been developed for use in image reconstruction [9] [13] . These methods are related to an L 1 - norm minimization of an underdetermined linear system. The use of a total variation (TV) norm as a constrained term is an effective approach for solving such an optimization problem [14] .
TV regularization aims to minimize the L 1 -norm of a gradient image showing the sparsifying domain. A CT image from a limited number of projections has been accurately reconstructed by using a TV minimization and an algebraic reconstruction technique [8] . Other first-order gradient approaches using forward and backward projection operators have been studied for improving algorithmic performance and convergence speed [9] [13] .
The iterative shrinkage-thresholding (IST) algorithm is derived from a consideration of the L 1 -norm of the proximal gradient method [15] and is now a common tool for image recovery — one that is based on the principle of CS [16] [17] . A soft-thresholding filtering algorithm using a pseudo-inverse of a discrete difference transform demonstrates a good image recovery [13] .
Recently, to improve convergence performance, a two-step iterative shrinkage-thresholding (TwIST) algorithm based on a second-order approach was developed [18] . In an iterative optimization process based on a Fourier space, the TwIST algorithm showed a convergence rate better than that of other first-order methods. However, there have been no reports that this algorithm is applicable to CT image reconstruction based on a Radon space.
In this research, we evaluate CT image reconstruction from a small amount of projection data using the TwIST algorithm, the results of which are compared with those of some first-order methods. The reconstruction properties of numerical phantom images using a parallel-beam geometry and a divergent-beam geometry are investigated, and the experimentally measured projection data are also used to carry out the analysis.
II. TwIST Algorithm for CT Image Recovery
Image recovery can be carried out by using the integral of the back projection of observed Radon projection data [19] , where the integral kernel is represented as a convolution of the projection data, p ( t , θ ), and filtering function, h ( t ).
u(x,y)= 0 π [ p(t,θ)h(t) ] t=xcosθ+ysinθ dθ
For a discrete matrix form, the pixel data of a 2D image, u ( x , y ), and projection data can be represented as a 1D column vector. Then, the above equation can be rewritten as
u p = A p T F 1 HFp ,
where F is a 2D matrix representation of a 1D discrete Fourier transform, and H is a Fourier-transformed filter term. Projection system matrix A p is related to the CT geometry.
If we define the filtered projection
p ^ (t,θ)
as F −1 HFp , especially in a noiseless projection data, then a simple linear equation is obtained as follows:
p= A p u .
The CT image recovery process is then mathematically represented by the above equation. In particular, for an insufficient amount of measurement data, a CS-based algorithm can successfully restore a sparse signal from an L 1 -norm minimization subject in accordance with the constraint in (3). The sparsity of an unknown signal and the incoherency of a system matrix are the main measure of a successful application of CS theory [5] , [7] . First, the desired image should have a sparse representation in a known transform domain. If the sparse signal, u sd , is represented by a sparsifying operator, Ψ , then (3) can then be written as
p= A p ψ u sd = A p u sd .
We also find that the coherence of system matrix
A p ′
should be very low for a good image recovery.
A CS recovery method that makes use of TV regularization becomes an alternative, where a sparsifying operator, such as a wavelet transform, is not required. The TV term can play the role of a sparsifying operator of a discrete gradient transform.
If the pixel values of the image are denoted by u m, n , then the TV norm is as follows [14] :
u m,n TV = m,n | u m,n | = m,n ( u m,n u m1,n ) 2 + ( u m,n u m,n1 ) 2 .
The TV norm is an L 1 -norm of the gradient image; the gradient image space is known as the sparse domain. TV regularization has been proven to be an effective method for image reconstruction from under-sampled data.
Image recovery can be conducted using a proximal gradient method for solving a convex optimization problem,
min  Aup 2 2 +λ u TV ,
where the objective function consists of the least squares together with the TV term, and λ is a regularized parameter. Generally, under an L 1 -norm penalty, the proximal gradient solution becomes the IST algorithm. The IST algorithm can be derived from the Majorize-Minimization technique or forward-backward splitting method [15] , [17] .
Considering a second-order regime, the TwIST algorithm was developed, which has been known to be a robust and stable image reconstruction method [18] . The iterative formula of the TwIST algorithm is represented as follows:
u k+1 =(1α) u k1 +(αβ) u k +β Γ λ ( u k ),
Γ λ (u)= Ψ λ ( u A T (Aup) ),
where Ψ becomes a TV operator. In our reconstruction algorithm, we adopt projection operators A and A T as forward and filtered back projections (FBP) based on (2), respectively. The parameters α and β are calculated from the spectral radius, ρ .
Let ξ be a real number such that 0 < ξ 1 Ki ( A T A ) ≤ ξm , where Ki (·) is the i th eigenvalue of its argument, and κ is defined as ξ 1 / ξm . Then, we have
ρ= 1 κ 1+ κ <1.
The performance of the algorithm is determined by parameters α and β , which can be expressed as follows:
α= ρ 2 +1,
β=2α/( ξ m + ξ 1 ).
The above parameters are defined values from an analysis of a system matrix eigenvalue; however, in a real simulation environment, the values are arbitrarily tuned for the optimal conditions.
III. Simulation Results
- 1. Parallel-Beam CT Image Reconstruction
We first evaluate the image recovery of a parallel-beam CT geometry by using the TwIST algorithm, and analyze the properties as compared to other first-order methods. The algorithms were written in MATLAB. Figure 1 shows an image reconstruction of a Shepp–Logan phantom with a size of 256 × 256 pixels. Projection data, which is acquired from 32 projection views, is uniformly distributed in the range [0, 2π], and used for a CS image reconstruction. A one-dimensional detector, having a pixel size of 367, can cover the projection data in a diagonal direction of an image according to a rotational angle. The pixels of both image and detector were set to have the same resolution, and this value does not affect the simulations.
The sinogram in the Radon space and the recovered image through the FBP prior to the application of CS theory are shown in Figs. 1(a) and 1(b) , respectively. The image retrieved by the FBP algorithm is in an inadequate state and shows many streaky artifacts owing to the projected geometry. As is well known, several-hundred projection views should be used to acquire a reliable image for a medical examination. However, we note that there is a fundamental difficulty in restoring the original information exactly using only the FBP algorithm. As shown in the inset of Fig. 1(c) , the boundary of the restored image through 360 projection views is smoothened. In a CT geometry, a smoothing of a restored image edge at a fine resolution is inevitable even in an image recovery using a greater number of projection views.
CS theory demonstrates that an original image value can be recovered accurately even under a relatively low number of projection views if the CS condition is satisfied [5] . Figure 1(d) displays the reconstructed image through the IST algorithm. As previously stated, the main iteration part of the algorithm is composed of forward and FBPs between a real space and a Radon space. A Hann filter was used for backward projection of the algorithm. The algorithm overshot in not using the filter, which will be discussed later. The magnified inset shows a better shape of the image edge compared to the FBP-recovered image shown in Fig. 1(c) despite the use of only 32 projection views. However, the reconstructed image does not show as sharp a boundary as the original image.
We carried out the image retrieval process using the fast iterative shrinkage-thresholding (FIST) algorithm, which is known as a CS algorithm with a higher performance level [16] . The projection views for this algorithm were chosen over π. A better image is obtained, as shown in Fig. 1(e) , but when applied to a CT image, the image quality is not largely enhanced. Figure 1(f) shows an image reconstructed through the TwIST algorithm, where the inset shows boundaries as sharp as those in the original image.
PPT Slide
Lager Image
Image recovery of parallel-beam CT geometry using iterative shrinkage algorithms: (a) sinogram, (b) FBP reconstruction for 32 projection views, (c) FBP reconstruction for 360 projection views, (d) IST reconstruction image, (e) FIST reconstruction image, and (f) TwIST reconstruction image for 32 projection views.
Figure 2 shows an image profile along a vertical line in the middle of the image. The restored image pixel values by the TwIST algorithm perfectly match those of the original image. It is very interesting that even when capturing only a few projection views, the image information can be retrieved within one-pixel resolution during CT image reconstruction.
PPT Slide
Lager Image
1D profiles along vertical line in middle of original image, FBP reconstruction image, and TwIST reconstruction image.
Figure 3 shows the convergence performance of the TwIST algorithm. The convergence is observed by the root-mean-square error (RMSE) as a function of the number of iterations. To increase the convergence performance, determining the optimal values of parameters α and β is crucial [18] . As described previously, the parameters are calculated by considering the eigenvalues of a system matrix. Under optimal conditions, α and β are 1.9 and 0.42, respectively. The acceptable range of these parameters is narrower than that of a Fourier-based system matrix, which will be discussed later by analyzing the property of the projection system matrix related to the Radon transform. The image reconstruction drastically converges until about 70 iterations, and then continues to converge gently. At near 70 iterations, a sharp boundary does not appear in the restored image.
The TwIST algorithm shows a higher convergence compared to other iterative methods shown in Fig. 3(b) . The convergence speed is several-orders of magnitude faster than the IST and FIST algorithms. To compare the convergence ability, all algorithms were optimized by varying the convergence variables. From this, we found that in CT image reconstruction, the TwIST algorithm is increasingly more stable and robust.
PPT Slide
Lager Image
Convergence rate of (a) TwIST algorithm and (b) comparison of convergence performance with other first-order methods.
- 2. Fan-Beam CT Image Reconstruction
Figure 4 shows the reconstructed results of fan-beam CT images through the previously described iterative shrinkage algorithms. The distances from the detector to fan-beam source and object image are 800 mm and 600 mm, respectively. The phantom image has a size of 256 × 256 pixels, and the one-dimensional detector array is 256 bins in size. The resolution of an image pixel is set to be 0.5 mm with respect to 1 mm of detector resolution and thus, the fan angle in this geometry is enough to cover the full detection of image projection data. The filtered back projection algorithm was composed from modifying the Feldkamp–Davis–Kress (FDK) algorithm — widely utilized for a cone-beam geometry [20] . Figure 4(a) shows an image reconstructed using the fan-beam FBP algorithm, where 64 projection views distributed over 2π are used. A Hann filter was used to prevent algorithm overshoot.
The IST and FIST algorithms reveal a similar retrieval behavior as the parallel-beam CT image recovery, as illustrated in Figs. 4(b) and 4(c) . However, in an image reconstructed using the TwIST algorithm, as shown in the magnified inset of Fig. 4(d) , there is no clear edge to the image. Therefore, we found that even the TwIST algorithm does not search the accurate information of the original image in a fan-beam geometry. Figure 5 shows the profile data of the vertical line of the reconstructed image. After 300 iterations, the image profile is recovered similarly to that of the original image, but the edge information is still insufficient.
PPT Slide
Lager Image
Reconstruction results of fan-beam CT image using iterative shrinkage algorithms: (a) fan-beam FBP reconstruction for 64 projection views, (b) IST reconstruction image, (c) FIST reconstruction image, and (d) TwIST reconstruction image.
PPT Slide
Lager Image
1D profiles along vertical line in middle of original image, fan-beam FBP reconstruction image, and TwIST reconstruction image.
As shown in Fig. 6 , which displays the variations of the RMSE as a function of the number of iterations, the minimum value is larger than that of the parallel-beam CT image for all other algorithms considered. The TwIST algorithm shows that the RMSE of the restored image sharply decreases until about 30 iterations, at which the value smoothly decreases. Parameters α and β were optimally found to be 1.6 and 0.5, respectively. Both parameters have a narrow window, similar to that of the parallel-beam geometry. Here, the IST and FIST algorithms were also optimized to compare their performances. Although the TwIST algorithm has a better convergence speed than that of other algorithms, a large enhancement of the algorithm performance does not appear.
PPT Slide
Lager Image
Comparison of convergence performance between TwIST algorithm and first-order methods for fan-beam CT image recovery.
Figure 7 shows the retrieval property of the anthropomorphic head phantom image using the TwIST algorithm, where α and β are reset to 1.6 and 0.3, respectively. This phantom image is much more complex in the spatial frequency domain and therefore indicates that the sparsity of the image is larger than that of the Shepp–Logan phantom. It can be seen that the image is well recovered when using 64 projection views, but the convergence rate is relatively slow. The algorithm also reveals a limitation in precisely finding the information of the original image, as shown in Fig. 7(b) of the central 1D profile data.
PPT Slide
Lager Image
Reconstruction property of anthropomorphic head phantom image using TwIST algorithm: (a) reconstructed image using fan-beam FBP and TwIST, and (b) image profiles along vertical line.
- 3. Experimental Phantom Study
Figure 8 shows the reconstruction characteristics of a cylindrical phantom image from the measured projection data using a cone-beam CT apparatus. The specifications of the cone-beam CT apparatus are briefly as follows: a flat-panel detector of 960 × 960 pixels is placed at a distance of 1,436 mm from the cone-beam source, and the distance between the cone-beam source and the object image is 1,000 mm. The pixel size is 0.444 mm. The projection data of 65 views measured from the cylindrical phantom are shown in Fig. 8(a) . Figure 8(b) shows a 512 × 512 pixels image reconstructed using the fan-beam FBP algorithm. To apply the proposed fan-beam algorithm, the central slice projection data was chosen. A high-quality image reconstructed from 320 projection views is displayed in Fig. 8(c) , where the interval of the rotational angle of the cone-beam source is set to 1.12°.
Figure 8(d) shows the image reconstructed from 65 projection views using the TwIST algorithm. The algorithm removes streaky artifacts and retrieves a clear image. We found that the TwIST algorithm shows a high image recovery from the projection data of a sparse number of views measuring the divergent-beam geometry.
PPT Slide
Lager Image
Reconstruction characteristics of cylindrical phantom image from measured projection data using cone-beam CT apparatus: (a) measured projection data, (b) fan-beam FBP reconstruction for 64 projection views, (c) fan-beam FBP reconstruction for 320 projection views, and (d) TwIST reconstruction image.
IV. Discussion
The performance of an image reconstruction using CS theory is determined through a condition known as the restricted isometry property (RIP) [6] , which is as follows:
(1 δ s ) u 2 2 || Au || 2 2 (1+ δ s ) u 2 2 .
For an effective retrieval of the image information, the isometry constant δs should be less than 1. A low δs indicates that the column vectors of system matrix A are nearly orthonormal. The following is another expression where the coherence of column vectors a i, j is very low [21]:
μ(A)= max ij |〈 a i , a j 〉|.
A Gaussian random matrix and Fourier matrix satisfy the above conditions well. A perfect image reconstruction is known to be achievable from the sparse Fourier-transform data, which demonstrates an exact reconstruction principle (ERP) [5] .
In CT geometry, the projection system matrix A p in (3) has a large condition number and an ill-posed problem is thus evident [8] . An interpolation process inducing data noise is inevitable because of the coordinate mapping between the image space and the Radon space. Therefore, appropriate regularization methods to prevent algorithm overshooting are required. In the iterative shrinkage algorithm, we confirmed that a filtering process of the back projection together with the TV regularization plays an important role in restraining this effect. Careful observation of algorithm performance according to filter species is needed — something we will look to address in future works.
A projection operator of a parallel-beam CT is assumed to be a simpler form compared with other geometries. The TwIST algorithm recovers the accurate information of an original image in a parallel-beam geometry. From this, we found that the coherence of the system matrix is sufficiently low to satisfy the ERP condition. On the other hand, the IST algorithm has difficulty in exactly recovering the original image even under optimal conditions. The algorithms based on a first-order approach showed a limitation for controlling overshoot facts. In fan-beam CT image recovery, an accurate image is not retrieved even by the TwIST algorithm. The fan-beam CT has a divergent beam geometry, and thus a rebinning process is imposed on the reconstruction process as compared to a parallel-beam geometry [19] . This is considered to lead to a poor performance in solving the measurement matrix.
The TwIST algorithm is well suited to the reconstruction of a CT image from sparse-view projection data. However, its ability to find exact image values does not reach that of a CS process using a simple Fourier transform. The incoherence of the system matrix is less than that of a simple Fourier sensing matrix. This fact is due to the complexity of the projection system matrix. In our parallel-beam geometry, the smallest and largest eigenvalues of the system matrix calculated from the optimal values of α and β appear to be 0.0063 and 9.041, respectively. Here, the maximum eigenvalue is greater than 1, which is anomalous in conventional convergence matrix systems. This may be due to an asymmetry in the forward and backward projection matrix. As mentioned, a filtering process in a backward projection is inevitably added. For a detailed understanding of this effect, a further analysis of the projection matrix is required. We also found that a retrieval method based on the Fourier basis of the system matrix will enhance the performance of algorithm convergence.
V. Conclusion
CT images were successfully recovered from sparse-view projection data by using the TwIST algorithm. The TwIST algorithm, based on a projection and filtered back projection, shows a better performance than other first-order methods. A parallel-beam CT image can be restored accurately, but in a divergent-beam geometry, a perfect recovery of the image was not accomplished, which is considered to be due to the complexity of the system matrix.
This work was supported by the R&D Convergence Program of NST (National Research Council of Science & Technology) of Republic of Korea (Development of Convergent Radio Therapy Equipment with O-arm CT, Grant CAP-12-9-KITECH).
Corresponding Author
Byung Gyu Chae received his BS, MS, and PhD degrees in physics from Pusan National University, Busan, Rep. of Korea, in 1993, 1995, and 1999, respectively. He joined ETRI in 2000. His main research interests are solid-state devices, holographic optics, and biomedical imaging. He currently studies next-generation biomedical imaging technology.
Sooyeul Lee received his BS, MS, and PhD degrees in physics from Seoul National University, Rep. of Korea, in 1990, 1992, and 1997, respectively. Since 1998, he has been with ETRI. His research interests include medical imaging and medical image processing.
Sidky E.Y. , Pan X.C. 2008 “Image Reconstruction in Circular Cone-Beam Computed Tomography by Constrained, Total-Variation Minimization,” Phys. Med. Biol. 52 (17) 4777 - 4807
Kudo H. , Suzuki T. , Rashed E.A. 2013 “Image Reconstruction for Sparse-View CT and Interior CT-Introduction to Compressed Sensing and Differentiated Backprojection,” Quant Imag. Med. Surgery 3 (3) 147 - 161
Chen G.H. , Tang J. , Leng S. 2008 “Prior Image Constrained Compressed Sensing (PICCS): a Method to Accurately Reconstruct Dynamic CT Images from Highly Undersampled Projection Data Sets,” Med. Phys. 35 (2) 660 - 663    DOI : 10.1118/1.2836423
Jaffray D.A. 2005 “Emergent Technologies for 3-Dimensional Image-Guided Radiation Delivery,” Seminars Radiation Oncology 15 (3) 208 - 216    DOI : 10.1016/j.semradonc.2005.01.003
Candes E.J. , Romberg J. , Tao T. 2006 “Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information,” IEEE Trans. Inf. Theory 52 (2) 489 - 509    DOI : 10.1109/TIT.2005.862083
Candes E.J. , Romberg J. , Tao T. 2006 “Stable Signal Recovery from Incomplete and Inaccurate Measurements,” Commun. Pure Appl. Math. 59 (8) 1207 - 1223    DOI : 10.1002/cpa.20124
Donoho D.L. 2006 “Compressed Sensing,” IEEE Trans. Inf. Theory 52 (4) 1289 - 1306    DOI : 10.1109/TIT.2006.871582
Sidky E.Y. , Kao C.M. , Pan X.C. 2006 “Accurate Image Reconstruction from Few-Views and Limited-Angle Data in Divergent-Beam CT,” J. X-Ray Sci. Technol. 14 (2) 119 - 139
Bian J. 2010 “Evaluation of Sparse-View Reconstruction from Flat-Panel-Detector Cone-Beam CT,” Phys. Med. Biol. 55 (22) 6575 - 6599    DOI : 10.1088/0031-9155/55/22/001
Park J.C. 2012 “Fast Compressed Sensing-Based CBCT Reconstruction Using Barzilai-Borwein Formulation for Application to on-Line IGRT,” Med. Phys. 39 (3) 1207 - 1217    DOI : 10.1118/1.3679865
Choi K. 2010 “Compressed Sensing Based Cone-Beam Computed Tomography Reconstruction with a First-Order Method,” Med. Phys. 37 (9) 5113 - 5125    DOI : 10.1118/1.3481510
Herman G.T. , Davidi R. 2008 “Image Reconstruction from a Small Number of Projections,” Inverse Problem 24 (4) 45011 - 45028    DOI : 10.1088/0266-5611/24/4/045011
Yu H. , Wang G. 2010 “A Soft-Threshold Filtering Approach for Reconstruction from a Limited Number of Projections,” Phys. Med. Biol. 55 (13) 3905 - 3916    DOI : 10.1088/0031-9155/55/13/022
Rudin L.I. , Osher S. , Fatemi E. 1992 “Nonlinear Total Variation Based Noise Removal Algorithm,” Physica D 60 259 - 268    DOI : 10.1016/0167-2789(92)90242-F
Daubechies I. , Defrise M. , De Mol C. 2004 “An Iterative Thresholding Algorithm for Linear Inverse Problems with a Sparsity Constraint,” Commun. Pure Appl. Math. 57 (11) 1413 - 1457    DOI : 10.1002/cpa.20042
Beck A. , Teboulle M. 2009 “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM J. Imag. Sci. 2 (1) 183 - 202    DOI : 10.1137/080716542
Combettes P.L. , Wajs V.R. 2006 “Signal Recovery by Proximal Forward-Backward Splitting,” Multiscale Modeling Simulation 4 (4) 1168 - 1200
Bioucas-Dias J. , Figueiredo M. 2007 “A New Twist: Two-Step Iterative Shrinkage/Thresholding Algorithms for Image Restoration,” IEEE Trans. Imag. Process. 16 (12) 2992 - 3004    DOI : 10.1109/TIP.2007.909319
Kak A.C. , Slaney M. 1988 “Principles of Computerized Tomographic Imaging,” IEEE Press New York, USA 49 - 112
Feldkamp L.A. , Davis L.C. , Kress J.W. 1984 “Practical Cone-Beam Algorithm,” J. Opt. Soc. America A 1 (6) 612 - 619    DOI : 10.1364/JOSAA.1.000612