Home | About | Journals | Submit | Contact Us | Français |

**|**Int J Biomed Imaging**|**v.2010; 2010**|**PMC2860338

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Mathematical Principles
- 3. Algorithm Development
- 4. Numerical Simulation
- 5. Discussions and Conclusion
- References

Authors

Related links

Int J Biomed Imaging. 2010; 2010: 934847.

Published online 2010 April 26. doi: 10.1155/2010/934847

PMCID: PMC2860338

SBES Division and ICTAS Center for Biomedical Imaging, VT-WFU School of Biomedical Engineering and Sciences, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA

*Hengyong Yu: Email: gro.eeei@uy-gnoygneh and

*Ge Wang: Email: gro.eeei@gnaw-eg

Academic Editor: Yibin Zheng

Received 2009 December 31; Accepted 2010 February 10.

Copyright © 2010 H. Yu and G. 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.

This article has been cited by other articles in PMC.

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.

Worldwide there are growing concerns on radiation induced genetic, cancerous, and other diseases [1–3]. 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
_{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
_{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.

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} ^{N} be an object function and **g** = [*g*
_{1},*g*
_{2},…,*g*
_{M}]^{T} ^{M} a set of measurements. Usually, they are linked by

$$\mathbf{g}=\mathbf{A}\mathbf{f}+\mathbf{e},$$

(1)

where **A** = (*a*
_{mn}) ^{M} × ^{N} is the linear measurement matrix, and **e** ^{M} the measurement noise. Let us define the
_{p} norm of the vector **g** as

$${||\mathbf{g}||}_{p}={\left({\displaystyle \sum}_{m=1}^{M}{g}_{m}^{p}\right)}^{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**)

$$\mathrm{\Delta}\left(\mathbf{f}\right)={||\mathbf{g}-\mathbf{A}\mathbf{f}||}^{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 (
_{γ})_{γΓ} of the space ^{N} satisfying **f** = ∑_{γΓ}**f**,
_{γ}
_{γ}, and a sequence of strictly positive weights **w** = (*w*
_{γ})_{γΓ}, we define the functional Φ_{w,p}(**f**) by

$${\mathrm{\Phi}}_{\mathbf{w},p}\left(\mathbf{f}\right)=\mathrm{\Delta}\left(\mathbf{f}\right)+{\displaystyle \sum}_{\gamma \in \mathrm{\Gamma}}2{w}_{\gamma}{\left|\langle \mathbf{f},{\mathbf{\phi}}_{\gamma}\rangle \right|}^{p},$$

(4)

where ·, · represents the inner product and 1 ≤ *p* ≤ 2. Using the
_{p} norm definition (2), let us define the
_{p} norm of a matrix operator **A** as

$${||\mathbf{A}||}_{p}=\underset{\mathbf{f}\ne 0}{\mathrm{max}\hspace{0.17em}}\left(\frac{{||\mathbf{A}\mathbf{f}||}_{p}}{{||\mathbf{f}||}_{p}}\right).$$

(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
_{p} norm regularization term ∑_{γΓ}2*w*
_{γ}|**f**,
_{γ}|^{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:

$${\mathbf{f}}^{k}={\mathbb{S}}_{\mathbf{w},p}\left({\mathbf{f}}^{k-1}+{\mathbf{A}}^{T}\left(\mathbf{g}-\mathbf{A}{\mathbf{f}}^{k-1}\right)\right),$$

(6)

where *k* = 1,2,… is the iteration number, **f**
^{0} the initial value in ^{N}, and

$${\mathbb{S}}_{\mathbf{w},p}\left(\mathbf{f}\right)={\displaystyle \sum}_{\gamma \in \mathrm{\Gamma}}{S}_{{w}_{\gamma},p}\left(\langle \mathbf{f},{\mathbf{\phi}}_{\gamma}\rangle \right){\mathbf{\phi}}_{\gamma}$$

(7)

with *S*
_{w,p} = (*F*
_{w,p})^{−1} being a one-to-one map from to its self for *p* > 1 with

$${F}_{w,p}\left(x\right)=x+wp\mathrm{sgn}\hspace{0.17em}\left(x\right){\left|x\right|}^{p-1}.$$

(8)

Particularly,

$${S}_{w,3/2}\left(x\right)=\{\begin{array}{cc}x-\frac{3w\left(\sqrt{9{w}^{2}+16\left|x\right|}-3w\right)}{8}\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}x\ge 0,\hfill \\ x+\frac{3w\left(\sqrt{9{w}^{2}+16\left|x\right|}-3w\right)}{8}\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}x<0.\hfill \end{array}$$

(9)

When *p* = 1, we can set [13]

$${S}_{w,1}\left(x\right)=\{\begin{array}{cc}x-w\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}x\ge w,\hfill \\ 0\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}\left|x\right|<w,\hfill \\ x+w\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}x\le -w.\hfill \end{array}$$

(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 *γ* Γ, Φ_{w,p}(**f**) can be rewritten as

$${\mathrm{\Phi}}_{\mathbf{w},p}\left(\mathbf{f}\right)={\mathrm{\Phi}}_{\tau ,p}\left(\mathbf{f}\right)=\mathrm{\Delta}\left(\mathbf{f}\right)+{\displaystyle \sum}_{\gamma \in \mathrm{\Gamma}}2\tau {\left|\langle \mathbf{f},{\mathbf{\phi}}_{\gamma}\rangle \right|}^{p}.$$

(11)

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

$$R\left({\mathbf{f}}^{\ast},p\right)={\left({\displaystyle \sum}_{\gamma \in \mathrm{\Gamma}}{\left|\langle {\mathbf{f}}^{\ast},{\mathbf{\phi}}_{\gamma}\rangle \right|}^{p}\right)}^{1/p},$$

(12)

which is the
_{p} norm radius of **f*** in the sparse space, we have the accelerated projected gradient algorithm

$${\mathbf{f}}^{k}={\mathbb{P}}_{R\left({\mathbf{f}}^{\ast},p\right)}\left({\mathbf{f}}^{k-1}+{\beta}^{k}{\mathbf{r}}^{k}\right),$$

(13)

where

$$\begin{array}{c}{\mathbf{r}}^{k}={\mathbf{A}}^{T}\left(\mathbf{g}-\mathbf{A}{\mathbf{f}}^{k-1}\right),\hspace{1em}\hspace{1em}{\beta}^{k}={||{\mathbf{r}}^{k}||}^{2}/{||\mathbf{A}{\mathbf{r}}^{k}||}^{2},\\ {\mathbb{P}}_{R\left({\mathbf{f}}^{\mathbf{\ast}},p\right)}\left(\mathbf{f}\right)={\mathbb{S}}_{\mu ,p}\left(\mathbf{f}\right)={\displaystyle \sum}_{\gamma \in \mathrm{\Gamma}}{S}_{\mu ,p}\left(\langle \mathbf{f},{\mathbf{\phi}}_{\gamma}\rangle \right){\mathbf{\phi}}_{\gamma},\end{array}$$

(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
_{R(f*,p)}(**f**) = **f**. When *R*(**f**, *p*) > *R*(**f***, *p*), the adapted threshold *μ* should be chosen to satisfy

$$R\left({\mathbb{P}}_{R\left({\mathbf{f}}^{\mathbf{\ast}},p\right)}\left(\mathbf{f}\right),p\right)=R\left({\mathbb{S}}_{\mu ,p}\left(\mathbf{f}\right),p\right)=R\left({\mathbf{f}}^{\mathbf{\ast}},p\right).$$

(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.

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 *n*th pixel can be viewed as a rectangular region with a constant value *f*
_{n}, the *m*th 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 *n*th pixel and the *m*th 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]

$${f}_{n}^{k}={f}_{n}^{k-1}+{\lambda}^{k}\frac{1}{{a}_{+n}}{\displaystyle \sum}_{m=1}^{M}\frac{{a}_{mn}}{{a}_{m+}}\left({g}_{m}-{\mathbf{A}}_{m}{\mathbf{f}}^{k-1}\right),$$

(16)

where *a*
_{+n} = ∑_{m=1}
^{M}
*a*
_{mn} > 0, *a*
_{m+} = ∑_{n=1}
^{N}
*a*
_{mn} > 0, **A**
_{m} is the *m*th row of **A**, *k* the iteration index, and 0 < *λ*
^{k} < 2 a free relaxation parameter. Let Λ^{+N} ^{N} × ^{N} be a diagonal matrix with Λ_{nn}
^{+N} = 1/*a*
_{+n} and let Λ^{M+} ^{M} × ^{M} be a diagonal matrix with Λ_{mm}
^{M+} = 1/*a*
_{m+}, then (16) can be rewritten as

$${\mathbf{f}}^{k}={\mathbf{f}}^{k-1}+{\lambda}^{k}{\tilde{\mathbf{r}}}^{k},$$

(17)

with

$${\tilde{\mathbf{r}}}^{k}={\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}\left(\mathbf{g}-\mathbf{A}{\mathbf{f}}^{k-1}\right).$$

(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

$${\mathbf{r}}^{k}=\frac{||{\mathbf{A}}^{T}||}{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}||}{\tilde{\mathbf{r}}}^{k}=\alpha {\tilde{\mathbf{r}}}^{k}.$$

(19)

$${\mathbf{f}}^{k}={\mathbb{P}}_{R\left({\mathbf{f}}^{\mathbf{\ast}},p\right)}\left({\mathbf{f}}^{k-1}+\alpha {\beta}^{k}{\tilde{\mathbf{r}}}^{k}\right),$$

(20)

with *β*
^{k} = ${||{\tilde{\mathbf{r}}}^{k}||}^{2}/{||\mathbf{A}{\tilde{\mathbf{r}}}^{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

$${\alpha}^{2}=\frac{{||{\mathbf{A}}^{T}||}^{2}}{{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}||}^{2}}=\frac{||{\mathbf{A}}^{T}||\xb7||\mathbf{A}||}{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}||\xb7||{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}||}.$$

(21)

Let **I** ^{N} be the vector with all whose components being “1”. We have ||**A**
^{T}
**A**||_{1} = ||**A**
^{T}
**A**||_{∞} = max_{1≤n≤N}(**A**
^{T}
**A**
**I**) because **A**
^{T}
**A** is a symmetric matrix. Hence, we have

$$||{\mathbf{A}}^{T}\mathbf{A}||\le \sqrt{{||{\mathbf{A}}^{T}\mathbf{A}||}_{1}{||{\mathbf{A}}^{T}\mathbf{A}||}_{\infty}}=\underset{1\le n\le N}{\mathrm{max}\hspace{0.17em}}\left({\mathbf{A}}^{T}\mathbf{A}\mathbf{I}\right).$$

(22)

Similarly, we have

$$\begin{array}{cc}\hfill & ||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}||\hfill \\ \hfill & \le \sqrt{{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}||}_{1}{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}||}_{\infty}}\hfill \\ \hfill & =\underset{1\le n\le N}{\mathrm{max}\hspace{0.17em}}\left({\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}\mathbf{I}\right).\hfill \end{array}$$

(23)

That is,

$${\alpha}^{2}=\frac{{\mathrm{max}\hspace{0.17em}}_{1\le n\le N}\hspace{0.17em}\left({\mathbf{A}}^{T}\mathbf{A}\mathbf{I}\right)}{\mathrm{max}\hspace{0.17em}\left(||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}\mathbf{I}||\right)}{\alpha}_{0}^{2},$$

(24)

with

$${\alpha}_{0}^{2}=\frac{{\mathrm{max}\hspace{0.17em}}_{1\le n\le N}\left({\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}\mathbf{I}\right)}{||{\mathrm{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathrm{\Lambda}}^{M+}||\xb7||{\mathrm{\Lambda}}^{M+}\mathbf{A}{\mathrm{\Lambda}}^{+N}||}\frac{||{\mathbf{A}}^{T}||\xb7||\mathbf{A}||}{{\mathrm{max}\hspace{0.17em}}_{1\le n\le N}\left({\mathbf{A}}^{T}\mathbf{A}\mathbf{I}\right)}.$$

(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 (
_{γ})_{γΓ} of the space ^{N}, in which **f** has a sparse representation. Our SART-type CT algorithm regularized by sparsity can be summarized in the following pseudocode:

- S1Initialize
*α*_{0},**f**^{(0)}and*k*; - S2Estimate
*R*(**f***,*p*); - S3Precompute
*α*,*a*_{+n}and*a*_{m+}; - S4Update the current estimation iteratively:
- S4.1
*k*: =*k*+ 1; - S4.2${\tilde{\mathbf{r}}}^{k}:={\mathbf{\Lambda}}^{+N}{\mathbf{A}}^{T}{\mathbf{\Lambda}}^{M+}(\mathbf{g}-\mathbf{A}{\mathbf{f}}^{k-1})$;
- S4.3${\beta}^{k}:={||{\tilde{\mathbf{r}}}^{k}||}^{2}/{||\mathbf{A}{\tilde{\mathbf{r}}}^{k}||}^{2}$;
- S4.4${\tilde{\mathbf{f}}}^{k}:={\mathbf{f}}^{k-1}+\alpha {\beta}^{k}{\tilde{\mathbf{r}}}^{k}$;
- S4.5Compute the sparse transform ${\varphi}_{\gamma}:=\langle {\tilde{\mathbf{f}}}^{k},{\mathbf{\phi}}_{\gamma}\rangle $ for
*γ*Γ; - S4.6Estimate the adapted threshold
*μ*; - S4.7Perform the soft-thresholding ${\tilde{\varphi}}_{\gamma}:={S}_{\mu ,p}({\varphi}_{\gamma})$;
- S4.8Perform the inverse sparse transform ${\mathbf{f}}^{k}:={\sum}_{\gamma \in \Gamma}{\tilde{\varphi}}_{\gamma}{\mathbf{\phi}}_{\gamma}$;

- 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]:

$${E}^{k}=\frac{||{\mathbf{f}}^{k}-{\mathbf{f}}^{\mathbf{\ast}}||}{||{\mathbf{f}}^{\mathbf{\ast}}||}\times 100.$$

(26)

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.0cm in fan-beam geometry. The object was a 128×128 modified Shepp-Logan phantom in a compact support with a radius of 10.0cm. We used an equispatial virtual detector array of length 20.0cm. 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
_{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.

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 **...**

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, **...**

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.

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
_{1}-norm would be small.

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
_{1}-norm seems giving the best performance than any other
_{p} norm with 1 < *p* ≤ 2 in a Haar space. However, it does not imply that the
_{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
_{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.

This work is partially supported by NIH/NIBIB Grants (EB002667, EB004287, and EB007288).

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**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |