Research Article  Open Access
SARTType Image Reconstruction from a Limited Number of Projections with the Sparsity Constraint
Abstract
Based on the recent mathematical findings on solving the linear inverse problems with sparsity constraints by Daubechiesx et al., here we adapt a simultaneous algebraic reconstruction technique (SART) for image reconstruction from a limited number of projections subject to a sparsity constraint in terms of an invertible compression transform. The algorithm is implemented with an exemplary Haar wavelet transform and tested with a modified SheppLogan phantom. Our preliminary results demonstrate that the sparsity constraint helps effectively improve the quality of reconstructed images and reduce the number of necessary projections.
1. Introduction
Worldwide there are growing concerns on radiation induced genetic, cancerous, and other diseases [1–3]. Computed tomography (CT) is considered as a radiationintensive procedure, yet it becomes more and more common. In the mid1990s, CT scans only accounted for 4% of the total Xray procedures but they contributed 40% of the collective dose [4]. The introduction of helical, multislice, and conebeam technologies has increased and continues the increasing usage of CT [5, 6]. In US, the number of CT examinations has been estimated as high as nearly 60 million in 2002, which account for 15% of imaging procedures and 75% of the radiation exposure [4]. In June 2007, the New York Times reported that “the percapita dose of ionizing radiation from clinical imaging exams in the U.S. increased almost 600% from 1980 to 2006.” More recently, in a highprofile article on the rapid growth in CT use and its associated radiation risks [3], Brenner and Hall estimated that “on the basis of such risk estimates and data on CT use from 1991 through 1996, it was estimated that about 0.4% of all cancers in the United States may be attributable to the radiation from CT studies. By adjusting this estimate for current CT use, this estimate might now be in the range of 1.5 to 2.0%.” Facing the increasing radiation risk, the wellknown As Low As Reasonably Achievable (ALARA) principle is widely accepted in the medical community. One of the practical strategies is to reduce the number of necessary projection.
Very interestingly, an elegant theory of compressive sampling or compressive sensing (CS) has recently emerged which shows that highquality signals and images can be reconstructed from far fewer measurements than what is usually considered necessary according to the Nyquist sampling theorem [7, 8]. The main idea of CS is that most signals are sparse in an appropriate orthonormal system; that is, a majority of their coefficients are close or equal to zero, when represented in the proper domain. Typically, CS starts with taking a limited amount of samples in a much less correlated basis, and then the signal is exactly recovered with an overwhelming probability from the limited amount of data via the norm minimization. For example, the discrete gradient sparsifying transform has been widely utilized in CSinspired CT reconstruction [9, 10], which was also referred to as the total variation minimization [11]. However, because the discrete gradient transform does not satisfy the restricted isometry property (RIP) required by the CS theory and is not invertible in general, such a reconstruction does not always convey the medically precise information. In particular, when a small number of projections are collected by a CT scanner, data noise may hide tumorlike structures in the TVminimizationbased reconstruction [12].
The above problem can be overcome using an invertible sparsifying transform such as a wavelet transform for image compression. For an object of interest such as a medical image, we can find an orthonormal basis (in a more general setting, a frame) to make the object sparse in terms of significant transform coefficients. Then, we can perform image reconstruction from a limited number of projections by minimizing the corresponding norm. Based on the recent mathematical findings made by Daubechies et al. [13, 14], here we will adapt a simultaneous algebraic reconstruction technique (SART) [15] for image reconstruction from a limited number of projections subject to a sparsity constraint in terms of an invertible sparsifying transform.
This paper is organized as follows. In the next section, the mathematical principles will be summarized. In the third section, a SARTtype reconstruction algorithm will be developed with a sparsity constraint. In the fourth section, preliminary numerical simulation results will be presented. In the last section, the related issues will be discussed.
2. Mathematical Principles
Daubechies and her collaborators proposed a general iterative thresholding algorithm to solve linear inverse problems regularized by a sparsity constraint and proved its convergence [13, 14]. Their approach can be directly applied for the CT reconstruction from a limited number of projections. Their major results can be summarized as follows.
Let = be an object function and = a set of measurements. Usually, they are linked by where is the linear measurement matrix, and the measurement noise. Let us define the norm of the vector as In practical applications, we usually omit the subscript when . To estimate from , one can minimize the discrepancy When system (1) is ill posed, the solution to (3) is not satisfactory, and additional constraints are required to regularize the solution. Particularly, given a complete basis of the space satisfying , and a sequence of strictly positive weights , we define the functional by where represents the inner product and . Using the norm definition (2), let us define the norm of a matrix operator as Let be the transpose matrix of , the operator in (1) is bounded, and . In the following, we will assume because can always be renormalized. To find an estimate of from under the norm regularization term , we can minimize defined in (4). The minimizer of can be recursively determined by the softthresholding algorithm: where is the iteration number, the initial value in , and with being a onetoone map from to its self for with Particularly, When , we can set [13] The main result of Daubechies et al. in [13] is that the solution of (6) is convergent.
Unfortunately, the convergence speed of (6) is very slow. To facilitate practical applications, an accelerated projected gradient method was then developed [14]. When for all , can be rewritten as Denote the minimizer of (11) as and define which is the norm radius of in the sparse space, we have the accelerated projected gradient algorithm where with an adapted softthreshold = depending on and . When , = and = . When , the adapted threshold should be chosen to satisfy Regarding algorithm (13), we have several comments in order. First, although Daubechies et al. only proved the convergence for the case [14], we believe that it should stand for . Second, while we have previously assumed that and , it can be proved that the algorithm (13) holds for any positive . Third, it is generally impossible to know the exact value of but we can have an approximate estimate.
3. Algorithm Development
In the context of image reconstruction, each component of the function in (1) is interpreted as an image pixel with being the total pixel number. Each component of the function is a measured datum with being the product of the number of projections and the number of detector elements. In fanbeam geometry with a discrete image grid, the th pixel can be viewed as a rectangular region with a constant value , the th measured datum as an integral of areas of pixels partially covered by a narrow beam from an Xray source to a detector element and respectively weighted by the corresponding Xray linear attenuation coefficients. Thus, the component in (1) can be understood as the interaction area between the th pixel and the th fanbeam path (Figure 1). While the whole matrix represents the forward projection, implements the back projection. The SARTtype solution to (1) can be written as [15] where , , is the th row of , the iteration index, and a free relaxation parameter. Let be a diagonal matrix with and let be a diagonal matrix with , then (16) can be rewritten as with
Due to the introduction of and , (18) cannot be directly applied to solve (13). However, we can modify (18) to obtain a new defined as Substituting (19) into (13), we have a SARTtype algorithm with = . The heuristic rationale for the above modification is to incorporate the SARTtype weighting scheme for a more uniform convergence behavior. Now, our task is to estimate . Since = , we have Let be the vector with all whose components being “1”. We have because is a symmetric matrix. Hence, we have Similarly, we have That is, with In practical applications, we can set to a reasonably large constant such as 2.0 in our simulation in the next section. If the algorithm does not converge, we can reduce until the algorithm converges.
For a basis of the space , in which has a sparse representation. Our SARTtype CT algorithm regularized by sparsity can be summarized in the following pseudocode:S1Initialize , and ;S2Estimate ;S3Precompute , and ;S4Update the current estimation iteratively:S4.1;S4.2;S4.3;S4.4;S4.5Compute the sparse transform for ;S4.6Estimate the adapted threshold ;S4.7Perform the softthresholding ;S4.8Perform the inverse sparse transform ;S5Go to S.4 until certain convergence criteria are satisfied.
In the above pseudocode, S.4.5 represents a sparse transform in a basis . It can be either orthonormal (e.g., Fourier transform) or nonorthonormal, and is the corresponding coefficient in the sparse space. In S.4.6, the adapted threshold can be estimated by a dichotomy searching method. S.4.7 performs the inverse sparse transform. Finally, the stopping criteria for S.5 can be either the maximum iteration number is reached or the relative reconstruction error (RRE) comes under a predefined threshold [14]:
4. Numerical Simulation
The aboveproposed algorithm was implemented in MatLab. To demonstrate its validity, we performed several numerical tests assuming a circular scanning locus of radius 57.0 cm in fanbeam geometry. The object was a 128 128 modified SheppLogan phantom in a compact support with a radius of 10.0 cm. We used an equispatial virtual detector array of length 20.0 cm. The detector was centered at the system origin and made perpendicular to the direction from the system origin to the Xray source. The detector array consisted of 128 elements. The wellknown “Haar” wavelet transform was selected to derive a sparse representation. While the pixel number of the original phantom image was 16384, there were only 1708 nonzero wavelet coefficients. In our preliminary study, the norm was focused as suggested by the CS theory. For each of our selected numbers of projections over a fullscan range, we equiangularly acquired the corresponding projection dataset based on the discrete projection model as shown in Figure 1. The stopping criterion in S.5 was defined as either the maximum iteration number 20,000 was reached or the RRE became less than 0.1%.
From each acquired dataset, we first reconstructed the image using the algorithm developed in Section 3, which is called “SchemeA.” For comparison, we also run our codes without S.4.6–S.4.8. This strategy is an adapted SARTtype reconstruction without regularization with the sparsity constraint, which is called “SchemeB.” In reality, the real solution is usually unknown, Daubechies et al. suggested an interior algorithm that could slowly increase the radius of the solution in each iteration step [14]. Thus, we also modified our algorithm described in Section 3 by replacing as in each iteration step to implement the corresponding version of interior algorithm, which is called “SchemeC.”
In the numerical simulation, after the stopping criteria were met, the iteration numbers and relative errors were listed in Table 1 with respect to different numbers of projections. The corresponding relative error convergence curves were plotted in Figure 2. The reconstructed images were shown in Figure 3. From the results in Table 1, Figures 2 and 3, we have several observations. When the number of projections was 55, SchemeA reached a 0.1% RRE after 19040 iterations. Because 0.1% was really small, the corresponding reconstructed image would serve as a gold standard for all other reconstructed images. First of all, in any tested cases either SchemeA or SchemeC performed far better than SchemeB, which confirmed that the sparse regularization did help improve the reconstructed image quality. Initially, the convergence speed of SchemeA was faster than SchemeC. However, after a number of iterations, the convergence speed of SchemeA became slower than SchemeC. If the ill posedness of the problem was not too bad, such as the cases of 55 and 45 projections, both SchemeA and SchemeC could perform well. When the problem was rather ill posed, such as the cases of 35 and 25 projections, SchemeC would perform better than SchemeA.

(a) 55 views
(b) 45 views
(c) 35 views
(d) 25 views
Compared to the original algorithm proposed by Daubechies et al., one unique feature of the proposed SARTtype algorithm is the weighting functions , and the associated constant for a more uniform converging behavior. To demonstrate this advantage, we modified our codes into a direct implementation of the algorithm described in [14] by forcing and setting and to unit diagonal matrices. The aforementioned reconstruction strategies were named as “SchemeAD,” “SchemeBD,” and “SchemeCD” and tested, respectively. The corresponding stopping conditions were listed in Table 2 with respect to the number of projections. The relative error curves were plotted in Figure 2. It can be observed in Figure 2 that when the problem was not too underdetermined, such as in the cases of 55 and 45 projections, the proposed methods did not perform significantly better, and might do even worse (e.g., SchemeAD was actually better than SchemeA in the case of 55 projections). When the problem was seriously underdetermined, such as in the cases of 35 and 25 projections, the proposed algorithms performed better than their direct implementation counterparts.

In practical applications, measurement noise is unavoidable. It is always important to use a stable algorithm for noisy data. To test the noise characteristic and stability of the proposed algorithms, we repeated the aforementioned reconstruction tests using “SchemeA,” “SchemeB,” and “SchemeC” with projections bearing 0.1% Gaussian noise, which are denoted as “SchemeAN,” “SchemeBN,” and “SchemeCN”, respectively. The corresponding stopping conditions were listed in Table 3 with respect to the number of projections. The converging curves were plotted in Figure 4. The reconstructed images were in Figure 5. It can be seen from the above results that the proposed algorithms produced similar relative error curves for noisy datasets compared to the noisefree counterparts. Due to noise in the projections, all the reconstructed images from noisy datasets generally had larger RREs than those from noisefree datasets given an iteration number. Note that SchemeAN had a better performance than “SchemeA” in the initial iterations in the case of 55 views. Our interpretation to this phenomenon is that in the initial iterations the discrepancy may be relatively large, which implies that may dominate the total cost functional , and the regularization effect of the norm would be small.

(a) 55 views
(b) 45 views
(c) 35 views
(d) 25 views
5. Discussions and Conclusion
Although the above CSbased reconstruction algorithms have been accelerated relative to the previous benchmark [14], the convergence speed is still slow for largescale images and/or very illposed conditions. In the future, we could use the stateoftheart computing techniques to speedup the convergence, such as orderedsubset [15] and multiscale computing [16] techniques. At the same time, we should optimize the reconstruction parameters and imaging protocols as well.
For the modified SheppLogan phantom, the norm seems giving the best performance than any other norm with in a Haar space. However, it does not imply that the norm is the best option for any application. In fact, our algorithm was implemented for any . For a specific application, the optimal may be studied.
Furthermore, this orthonormality of the wavelet transform used in this study is not necessary. If an image can be sparsely expanded in a certain basis or frame, the norm minimization can be in principle performed to regularize the reconstruction process. Since there exist many compression methods for medical images, we should evaluate representative bases and frames for sparse representations and CSbased reconstruction methods. The heuristic rule is to achieve a minimal compression ratio.
In conclusion, we have developed a SARTtype reconstruction algorithm based on the recent mathematical findings by Daubechies et al. Our preliminary simulation results have confirmed its merits and suggested research directions. Because the approach accommodates any and any sparse expansion, there should be a large room for further improvements of the algorithm performance.
Acknowledgment
This work is partially supported by NIH/NIBIB Grants (EB002667, EB004287, and EB007288).
References
 D. J. Brenner, C. D. Elliston, E. J. Hall, and W. E. Berdon, “Estimated risks of radiationinduced fatal cancer from pediatric CT,” American Journal of Roentgenology, vol. 176, no. 2, pp. 289–296, 2001. View at: Google Scholar
 A. B. de Gonzalez and S. Darby, “Risk of cancer from diagnostic Xrays: estimates for the UK and 14 other countries,” The Lancet, vol. 363, no. 9406, pp. 345–351, 2004. View at: Publisher Site  Google Scholar
 D. J. Brenner and E. J. Hall, “Current concepts—computed tomography—an increasing source of radiation exposure,” The New England Journal of Medicine, vol. 357, no. 22, pp. 2277–2284, 2007. View at: Google Scholar
 O. W. Linton and F. A. Mettler Jr., “National conference on dose reduction in CT, with an emphasis on pediatric patients,” American Journal of Roentgenology, vol. 181, no. 2, pp. 321–329, 2003. View at: Google Scholar
 G. Wang, Y. Ye, and H. Yu, “Approximate and exact conebeam reconstruction with standard and nonstandard spiral scanning,” Physics in Medicine and Biology, vol. 52, no. 6, pp. R1–R13, 2007. View at: Publisher Site  Google Scholar
 G. Wang, H. Yu, and B. De Man, “An outlook on xray CT research and development,” Medical Physics, vol. 35, no. 3, pp. 1051–1064, 2008. View at: Publisher Site  Google Scholar
 D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006. View at: Publisher Site  Google Scholar
 E. J. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp. 489–509, 2006. View at: Publisher Site  Google Scholar
 G.H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Medical Physics, vol. 35, no. 2, pp. 660–663, 2008. View at: Publisher Site  Google Scholar
 H. Yu and G. Wang, “Compressed sensing based interior tomography,” Physics in Medicine and Biology, vol. 54, no. 9, pp. 2791–2805, 2009. View at: Publisher Site  Google Scholar
 L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D, vol. 60, no. 1–4, pp. 259–268, 1992. View at: Google Scholar
 G. T. Herman and R. Davidi, “Image reconstruction from a small number of projections,” Inverse Problems, vol. 24, no. 4, Article ID 045011, 17 pages, 2008. View at: Publisher Site  Google Scholar
 I. Daubechies, M. Defrise, and C. De Mol, “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint,” Communications on Pure and Applied Mathematics, vol. 57, no. 11, pp. 1413–1457, 2004. View at: Publisher Site  Google Scholar
 I. Daubechies, M. Fornasier, and I. Loris, “Accelerated projected gradient method for linear inverse problems with sparsity constraints,” Journal of Fourier Analysis and Applications, vol. 14, no. 56, pp. 764–792, 2008. View at: Publisher Site  Google Scholar
 G. Wang and M. Jiang, “Orderedsubset simultaneous algebraic reconstruction techniques (OSSART),” Journal of XRay Science and Technology, vol. 12, no. 3, pp. 169–177, 2004. View at: Google Scholar
 S. Oh, A. B. Milstein, C. A. Bouman, and K. J. Webb, “A general framework for nonlinear multigrid inversion,” IEEE Transactions on Image Processing, vol. 14, no. 1, pp. 125–140, 2005. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2010 Hengyong Yu and Ge Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.