PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of ijbiInternational Journal of Biomedical Imaging
 
Int J Biomed Imaging. 2010; 2010: 934847.
Published online 2010 April 26. doi:  10.1155/2010/934847
PMCID: PMC2860338

SART-Type 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 Shepp-Logan 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 [13]. Computed tomography (CT) is considered as a radiation-intensive procedure, yet it becomes more and more common. In the mid-1990s, CT scans only accounted for 4% of the total X-ray procedures but they contributed 40% of the collective dose [4]. The introduction of helical, multislice, and cone-beam 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 per-capita dose of ionizing radiation from clinical imaging exams in the U.S. increased almost 600% from 1980 to 2006.” More recently, in a high-profile 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 well-known 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 high-quality 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 [ell] 1 norm minimization. For example, the discrete gradient sparsifying transform has been widely utilized in CS-inspired 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 tumor-like structures in the TV-minimization-based 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 [ell] 1 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 SART-type 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 f = [f 1,f 2,…,f N]T [set membership] RN be an object function and g = [g 1,g 2,…,g M]T [set membership] RM a set of measurements. Usually, they are linked by

g=Af+e,
(1)

where A = (a mn) [set membership] RM × RN is the linear measurement matrix, and e [set membership] RM the measurement noise. Let us define the [ell] p norm of the vector g as

||g||p=(m=1Mgmp)1/p.
(2)

In practical applications, we usually omit the subscript p when p = 2. To estimate f from g, one can minimize the discrepancy Δ(f)

Δ(f)=||gAf||2.
(3)

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 ([var phi] γ)γ[set membership]Γ of the space RN satisfying f = ∑γ[set membership]Γleft angle bracketf, [var phi] γright angle bracket[var phi] γ, and a sequence of strictly positive weights w = (w γ)γ[set membership]Γ, we define the functional Φw,p(f) by

Φw,p(f)=Δ(f)+γΓ2wγ|f,φγ|p,
(4)

where left angle bracket·, ·right angle bracket represents the inner product and 1 ≤ p ≤ 2. Using the [ell] p norm definition (2), let us define the [ell] p norm of a matrix operator A as

||A||p=maxf0(||Af||p||f||p).
(5)

Let A T be the transpose matrix of A, the operator A in (1) is bounded, and ||A T A|| < C. In the following, we will assume C = 1 because A can always be renormalized. To find an estimate of f from g under the [ell] p norm regularization term ∑γ[set membership]Γ2w γ|left angle bracketf,[var phi] γright angle bracket|p, we can minimize Φw,p(f) defined in (4). The minimizer of Φw,p(f) can be recursively determined by the soft-thresholding algorithm:

fk=𝕊w,p(fk1+AT(gAfk1)),
(6)

where k = 1,2,… is the iteration number, f 0 the initial value in RN, and

𝕊w,p(f)=γΓSwγ,p(f,φγ)φγ
(7)

with S w,p = (F w,p)−1 being a one-to-one map from R to its self for p > 1 with

Fw,p(x)=x+wpsgn(x)|x|p1.
(8)

Particularly,

Sw,3/2(x)={x3w(9w2+16|x|3w)8ifx0,x+3w(9w2+16|x|3w)8ifx<0.
(9)

When p = 1, we can set [13]

Sw,1(x)={xwifxw,0if|x|<w,x+wifxw.
(10)

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 w γ = τ for all γ [set membership] Γ, Φw,p(f) can be rewritten as

Φw,p(f)=Φτ,p(f)=Δ(f)+γΓ2τ|f,φγ|p.
(11)

Denote the minimizer of (11) as f* and define

R(f,p)=(γΓ|f,φγ|p)1/p,
(12)

which is the [ell] p norm radius of f* in the sparse space, we have the accelerated projected gradient algorithm

fk=R(f,p)(fk1+βkrk),
(13)

where

rk=AT(gAfk1),βk=||rk||2/||Ark||2,R(f,p)(f)=𝕊μ,p(f)=γΓSμ,p(f,φγ)φγ,
(14)

with an adapted soft-threshold μ = μ(R(f*, p), f) depending on R(f*, p) and f. When R(f, p) ≤ R(f*, p), μ(R(f*, p), f) = 0 and P R(f*,p)(f) = f. When R(f, p) > R(f*, p), the adapted threshold μ should be chosen to satisfy

R(R(f,p)(f),p)=R(𝕊μ,p(f),p)=R(f,p).
(15)

Regarding algorithm (13), we have several comments in order. First, although Daubechies et al. only proved the convergence for the case p = 1 [14], we believe that it should stand for 1 ≤ p ≤ 2. Second, while we have previously assumed that ||A T A|| < C and C = 1, it can be proved that the algorithm (13) holds for any positive C. Third, it is generally impossible to know the exact value of R(f*, p) but we can have an approximate estimate.

3. Algorithm Development

In the context of image reconstruction, each component of the function f in (1) is interpreted as an image pixel with N being the total pixel number. Each component of the function g is a measured datum with M being the product of the number of projections and the number of detector elements. In fan-beam geometry with a discrete image grid, the nth pixel can be viewed as a rectangular region with a constant value f n, the mth measured datum g m as an integral of areas of pixels partially covered by a narrow beam from an X-ray source to a detector element and respectively weighted by the corresponding X-ray linear attenuation coefficients. Thus, the component a mn in (1) can be understood as the interaction area between the nth pixel and the mth fan-beam path (Figure 1). While the whole matrix A represents the forward projection, A T implements the back projection. The SART-type solution to (1) can be written as [15]

Figure 1
Projection model of a discrete image in fan-beam geometry.

fnk=fnk1+λk1a+nm=1Mamnam+(gmAmfk1),
(16)

where a +n = ∑m=1 M a mn > 0, a m+ = ∑n=1 N a mn > 0, A m is the mth row of A, k the iteration index, and 0 < λ k < 2 a free relaxation parameter. Let Λ+N [set membership] RN × RN be a diagonal matrix with Λnn +N = 1/a +n and let ΛM+ [set membership] RM × RM be a diagonal matrix with Λmm M+ = 1/a m+, then (16) can be rewritten as

fk=fk1+λkr˜k,
(17)

with

r˜k=Λ+NATΛM+(gAfk1).
(18)

Due to the introduction of Λ +N and Λ M+, (18) cannot be directly applied to solve (13). However, we can modify (18) to obtain a new r k defined as

rk=||AT||||Λ+NATΛM+||r˜k=αr˜k.
(19)

Substituting (19) into (13), we have a SART-type algorithm

fk=R(f,p)(fk1+αβkr˜k),
(20)

with β k = ||r˜k||2/||Ar˜k||2. The heuristic rationale for the above modification is to incorporate the SART-type weighting scheme for a more uniform convergence behavior. Now, our task is to estimate α. Since ||A T|| = ||A||, we have

α2=||AT||2||Λ+NATΛM+||2=||AT||·||A||||Λ+NATΛM+||·||ΛM+AΛ+N||.
(21)

Let I [set membership] RN be the vector with all whose components being “1”. We have ||A T A||1 = ||A T A|| = max 1≤nN(A T A I) because A T A is a symmetric matrix. Hence, we have

||ATA||||ATA||1||ATA||=max1nN(ATAI).
(22)

Similarly, we have

||Λ+NATΛM+ΛM+AΛ+N||||Λ+NATΛM+ΛM+AΛ+N||1||Λ+NATΛM+ΛM+AΛ+N||=max1nN(Λ+NATΛM+ΛM+AΛ+NI).
(23)

That is,

α2=max1nN(ATAI)max(||Λ+NATΛM+ΛM+AΛ+NI||)α02,
(24)

with

α02=max1nN(Λ+NATΛM+ΛM+AΛ+NI)||Λ+NATΛM+||·||ΛM+AΛ+N||||AT||·||A||max1nN(ATAI).
(25)

In practical applications, we can set α 0 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 α 0 until the algorithm converges.

For a basis ([var phi] γ)γ[set membership]Γ of the space RN, in which f has a sparse representation. Our SART-type CT algorithm regularized by sparsity can be summarized in the following pseudocode:

  • S1
    Initialize α 0, f (0) and k;
  • S2
    Estimate R(f*, p);
  • S3
    Precompute α, a +n and a m+;
  • S4
    Update the current estimation iteratively:
    • S4.1
      k : = k + 1;
    • S4.2
      r˜k:=Λ+NATΛM+(g-Afk-1);
    • S4.3
      βk:=||r˜k||2/||Ar˜k||2;
    • S4.4
      f˜k:=fk-1+αβkr˜k;
    • S4.5
      Compute the sparse transform ϕγ:=f˜k,φγ for γ [set membership] Γ;
    • S4.6
      Estimate the adapted threshold μ;
    • S4.7
      Perform the soft-thresholding ϕ˜γ:=Sμ,p(ϕγ);
    • S4.8
      Perform the inverse sparse transform fk:=γΓϕ˜γφγ;
  • S5
    Go to S.4 until certain convergence criteria are satisfied.

In the above pseudocode, S.4.5 represents a sparse transform in a basis ([var phi] γ)γ[set membership]Γ. 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]:

Ek=||fkf||||f||×100.
(26)

4. Numerical Simulation

The above-proposed 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 fan-beam geometry. The object was a 128 × 128 modified Shepp-Logan 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 X-ray source. The detector array consisted of 128 elements. The well-known “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 [ell] 1-norm was focused as suggested by the CS theory. For each of our selected numbers of projections over a full-scan 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 “Scheme-A.” For comparison, we also run our codes without S.4.6–S.4.8. This strategy is an adapted SART-type reconstruction without regularization with the sparsity constraint, which is called “Scheme-B.” In reality, the real solution f* 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 R(f*, p) as R k = (0.4 + 0.6(k/20000)0.05)R(f*, p) in each iteration step k to implement the corresponding version of interior algorithm, which is called “Scheme-C.”

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 Figures2 and2 and and3,3, we have several observations. When the number of projections was 55, Scheme-A 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 Scheme-A or Scheme-C performed far better than Scheme-B, which confirmed that the sparse regularization did help improve the reconstructed image quality. Initially, the convergence speed of Scheme-A was faster than Scheme-C. However, after a number of iterations, the convergence speed of Scheme-A became slower than Scheme-C. If the ill posedness of the problem was not too bad, such as the cases of 55 and 45 projections, both Scheme-A and Scheme-C could perform well. When the problem was rather ill posed, such as the cases of 35 and 25 projections, Scheme-C would perform better than Scheme-A.

Figure 2
Relative error curves for various reconstructions from noise-free projections. The numbers of projections were set to (a) 55, (b) 45, (c) 35, and (d) 25, respectively, where the horizontal and vertical axes represent the iteration index k and the corresponding ...
Figure 3
Reconstructed 128 × 128 images from noise-free projection datasets. The 1st, 2nd, and 3rd rows are reconstructed using Scheme-A, Scheme-B, and Scheme-C, respectively, and the 1st, 2nd, 3rd, and 4th columns are from 55, 45, 35, ...
Table 1
Maximum iteration numbers and relative errors associated with each of three representative reconstruction schemes for different numbers of noise-free projections.

Compared to the original algorithm proposed by Daubechies et al., one unique feature of the proposed SART-type algorithm is the weighting functions ΛM+, Λ+N 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 α = 1.0 and setting ΛM+ and Λ+N to unit diagonal matrices. The aforementioned reconstruction strategies were named as “Scheme-AD,” “Scheme-BD,” and “Scheme-CD” 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 under-determined, 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., Scheme-AD was actually better than Scheme-A in the case of 55 projections). When the problem was seriously under-determined, such as in the cases of 35 and 25 projections, the proposed algorithms performed better than their direct implementation counterparts.

Table 2
Same as Table 1 but for the direct implementation counterparts of the proposed algorithms.

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 “Scheme-A,” “Scheme-B,” and “Scheme-C” with projections bearing 0.1% Gaussian noise, which are denoted as “Scheme-AN,” “Scheme-BN,” and “Scheme-CN”, 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 noise-free counterparts. Due to noise in the projections, all the reconstructed images from noisy datasets generally had larger RREs than those from noise-free datasets given an iteration number. Note that Scheme-AN had a better performance than “Scheme-A” in the initial iterations in the case of 55 views. Our interpretation to this phenomenon is that in the initial iterations the discrepancy Δ(f) may be relatively large, which implies that Δ(f) may dominate the total cost functional Φw,p(f), and the regularization effect of the [ell] 1-norm would be small.

Figure 4
Relative error curves for various reconstructions from noisy projection. The numbers of projections were set to (a) 55, (b) 45, (c) 35, and (d) 25, respectively, where the horizontal and vertical axes represent the iteration index k and the corresponding ...
Figure 5
Counterpart of Figure 3 in the case of noisy projection datasets.
Table 3
Same as Table 1 but for projection data bearing 0.1% Gaussian noise.

5. Discussions and Conclusion

Although the above CS-based reconstruction algorithms have been accelerated relative to the previous benchmark [14], the convergence speed is still slow for large-scale images and/or very ill-posed conditions. In the future, we could use the state-of-the-art computing techniques to speedup the convergence, such as ordered-subset [15] and multiscale computing [16] techniques. At the same time, we should optimize the reconstruction parameters and imaging protocols as well.

For the modified Shepp-Logan phantom, the [ell] 1-norm seems giving the best performance than any other [ell] p norm with 1 < p ≤ 2 in a Haar space. However, it does not imply that the [ell] 1-norm is the best option for any application. In fact, our algorithm was implemented for any 1 ≤ p ≤ 2. For a specific application, the optimal p 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 [ell] p-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 CS-based reconstruction methods. The heuristic rule is to achieve a minimal compression ratio.

In conclusion, we have developed a SART-type 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 1 ≤ p ≤ 2 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

1. Brenner DJ, Elliston CD, Hall EJ, Berdon WE. Estimated risks of radiation-induced fatal cancer from pediatric CT. American Journal of Roentgenology. 2001;176(2):289–296. [PubMed]
2. de Gonzalez AB, Darby S. Risk of cancer from diagnostic X-rays: estimates for the UK and 14 other countries. The Lancet. 2004;363(9406):345–351. [PubMed]
3. Brenner DJ, Hall EJ. Current concepts—computed tomography—an increasing source of radiation exposure. The New England Journal of Medicine. 2007;357(22):2277–2284. [PubMed]
4. Linton OW, Mettler FA., Jr. National conference on dose reduction in CT, with an emphasis on pediatric patients. American Journal of Roentgenology. 2003;181(2):321–329. [PubMed]
5. Wang G, Ye Y, Yu H. Approximate and exact cone-beam reconstruction with standard and non-standard spiral scanning. Physics in Medicine and Biology. 2007;52(6):R1–R13. [PubMed]
6. Wang G, Yu H, De Man B. An outlook on x-ray CT research and development. Medical Physics. 2008;35(3):1051–1064. [PubMed]
7. Donoho DL. Compressed sensing. IEEE Transactions on Information Theory. 2006;52(4):1289–1306.
8. Candes EJ, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory. 2006;52(2):489–509.
9. Chen G-H, Tang J, Leng S. Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets. Medical Physics. 2008;35(2):660–663. [PMC free article] [PubMed]
10. Yu H, Wang G. Compressed sensing based interior tomography. Physics in Medicine and Biology. 2009;54(9):2791–2805. [PMC free article] [PubMed]
11. Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D. 1992;60(1–4):259–268.
12. Herman GT, Davidi R. Image reconstruction from a small number of projections. Inverse Problems. 2008;24(4):17 pages. Article ID 045011. [PMC free article] [PubMed]
13. Daubechies I, Defrise M, De Mol C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics. 2004;57(11):1413–1457.
14. Daubechies I, Fornasier M, Loris I. Accelerated projected gradient method for linear inverse problems with sparsity constraints. Journal of Fourier Analysis and Applications. 2008;14(5-6):764–792.
15. Wang G, Jiang M. Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART) Journal of X-Ray Science and Technology. 2004;12(3):169–177.
16. Oh S, Milstein AB, Bouman CA, Webb KJ. A general framework for nonlinear multigrid inversion. IEEE Transactions on Image Processing. 2005;14(1):125–140. [PubMed]

Articles from International Journal of Biomedical Imaging are provided here courtesy of Hindawi Publishing Corporation