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

**|**Int J Biomed Imaging**|**v.2011; 2011**|**PMC2943093

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Imaging Model
- 3. Reconstruction Algorithm
- 4. Numerical Simulation
- 5. Discussions and Conclusion
- References

Authors

Related links

Int J Biomed Imaging. 2011; 2011: 549537.

Published online 2010 September 5. doi: 10.1155/2011/549537

PMCID: PMC2943093

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

Academic Editor: Yangbo Ye

Received 2010 May 7; Accepted 2010 July 29.

Copyright © 2011 Hengyong Yu et al.

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.

To maximize the time-integrated X-ray flux from multiple X-ray sources and shorten the data acquisition process, a promising way is to allow overlapped projections from multiple sources being simultaneously on without involving the source multiplexing technology. The most challenging task in this configuration is to perform image reconstruction effectively and efficiently from overlapped projections. Inspired by the single-source simultaneous algebraic reconstruction technique (SART), we hereby develop a multisource SART-type reconstruction algorithm regularized by a sparsity-oriented constraint in the soft-threshold filtering framework to reconstruct images from overlapped projections. Our numerical simulation results verify the correctness of the proposed algorithm and demonstrate the advantage of image reconstruction from overlapped projections.

Since the first computed tomography (CT) scanner was made [1], all the commercial scanners have been employing the X-ray source with a single small focal spot, which can be modeled as a point source. In micro-CT and even nano-CT applications, the reduced X-ray focal spot size becomes increasingly a limiting factor to achieve contrast and temporal resolution targets. To address this issue, our group recently proposed to use a line-shaped X-ray source so that more photons can be generated in a given data acquisition interval [2]. In this context, the X-ray source can be modeled as a line-segment, which can be further discretely modeled as an array of points [3]. In single source CT scanners, contrast resolution is limited by the finite focal-spot size necessary to generate a sufficient number of X-ray photons, and temporal resolution is limited by the time taken to acquire sufficiently many projections over a full-scan or half-scan angular range. Since a line source covers a wide angular range per view, the number of photons is increased to radiate an object to be reconstructed. Therefore, use of a line-shaped X-ray source or a multiple source array is a candidate scheme to balance spatial, contrast and temporal resolution.

Interestingly, the technology of field-emission X-ray sources based on carbon nanotubes (CNT) is a recent invention with several intrinsic advantages over conventional X-ray tubes. To maximize the time-integrated X-ray flux from multiple sources and improve the signal-to-noise ratio (SNR), multiplexing was used, where multi-source are excited with different temporal modulations [4]. If two or more sources are simultaneously fired, X-ray photons reach the same detector together for any single measurement, and one would not be able to identify which photon comes from which source. To unmix the signals from various sources, several projections can be collected at different time instants for the same view [5]. Due to the limited readout speed and mixed signals, there seems little advantage to use multiplexing configuration in terms of contrast and temporal resolution. In collaboration with Dr. Otto Zhou's group at University of North Carolina, we are developing a novel multi-source micro-CT system. In our system, we plan to fire a multi-source system simultaneously and acquire multiple projections at any viewing angle on nonoverlapped segments of a shared detector array [6]. With recent developed compressive sampling (CS) techniques [7, 8], we are also working to improve temporal resolution and reduce radiation dose from a limited number of nonoverlapped projections [9].

Here, we consider how to reconstruct an image from overlapped projections. Previously, our group developed a generalized simultaneous algebraic reconstruction technique (SART) algorithm to reconstruct an image from data collected with an X-ray line source [2]. Assuming that the differences between measured and predicted projections from various source points have equal weights, these differences are then backprojected to different X-ray source points in the SART framework. However, this algorithm converges to the least square solution that is not necessarily the true image. In our simulation, the reconstructed images sometimes suffered from blurring [2]. Then, our group developed a CS-based algorithm to solve this problem [3]. The algorithm is implemented in a projection onto convex sets (POCS) framework and employed a steepest gradient searching strategy. Although this algorithm often converges to the true image, its convergence speed is rather slow. Because the SART framework has an excellent convergence behavior especially when the ordered subset (OS) format is applied, in this paper we will develop a SART-type algorithm for image reconstruction from overlapped projections.

The paper is organized as follows. In the next section, we will formulate the overlapped projection model, which is not a sum of line integrals and quite different from that in the single-source case. In the third section, we will design a SART-type reconstruction algorithm for image reconstruction from overlapped projections. In the fourth section, we will report numerical results. In the last section, we will discuss related issues and conclude the paper.

For CT reconstruction, a two-dimensional digital image can be expressed as **f** = (*f*_{i,j}) ^{I} × ^{J}, where the index 1 ≤ *i* ≤ *I* and 1 ≤ *j* ≤ *J* are integers. Define

$${f}_{n}={f}_{i,j},\hspace{1em}n=\left(i-1\right)\times J+j,$$

(1)

with 1 ≤ *n* ≤ *N* and *N* = *I* × *J*, we have the image in a vector representation **f** = [*f*_{1},*f*_{2},...,*f*_{N}]^{T} ^{N}. In this paper, we will use both the signs *f*_{i,j} and *f*_{n} for convenience. Let **g** = [*g*_{1},*g*_{2},...,*g*_{M}]^{T} ^{M} be a measured vector with *M* being the product of the number of projections and the number of detector elements. They are linked by the following linear system:

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

(2)

where **A** = (*a*_{m,n}) ^{M} × ^{N} is a measurement matrix. In a typical fan-beam geometry, 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 partially covered pixel areas by a narrow beam from an X-ray source to a detector element which are weighted by the corresponding X-ray linear attenuation coefficients. Thus, the component *a*_{m,n} in (2) denotes the intersection area between the *n*th pixel and the *m*th fan-beam ray (Figure 1). While the whole matrix **A** represents the forward projection, **A**^{T} implements the backprojection.

While the imaging model (2) is valid for a single-source system, it cannot be used for multi-source-generated projections. The reason is that the measured data in (2) has been postprocessed by a logarithmic operation. In other words, we must model the raw data directly. For a multi-source system with *Q* sources, assuming that the *q*th source emits *I*_{q} photons towards each detector element. According to Beer's law, we have the following imaging model:

$$\mathbf{p}={\displaystyle \sum}_{q=1}^{Q}{I}_{q}{e}^{-s{\mathbf{A}}^{q}\mathbf{f}},$$

(3)

where **A**^{q} is the system matrix defined in Section 2.1 for the *q*th source, *s* a constant to normalize the area model (Figure 1) for the measurement matrix to the line integral model, and *e*^{−s Aqf} is a vector whose element is the exponential function of the corresponding element of −*s ***A**^{q}**f**. In practice, we can approximate *s* as the reciprocal of the average width of an X-ray path in an object to be reconstructed. If we assume that all the X-ray sources have the same intensities, *I*_{1} = *I*_{2} = = *I*_{q} = *I*_{0}, we have

$$\tilde{\mathbf{p}}=\frac{\mathbf{p}}{{I}_{0}}={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}\tilde{\mathbf{f}}},$$

(4)

where the constant *s* is absorbed by $\tilde{\mathbf{f}}=s\mathbf{f}$, and 1/*I*_{0} absorbed by $\tilde{\mathbf{p}}=\tilde{\mathbf{p}}/{I}_{0}$. The key task in this paper is to reconstruct $\tilde{\mathbf{f}}$ from $\tilde{\mathbf{p}}$ and known measurement matrices **A**^{q}.

In the past decade, our group studied a block-iterative (BI) or ordered-subset (OS) version of a general Landweber scheme [10], of which the SART and OS-SART [11] are special cases, for minimization of a weighted least square functional in the real/complex space, and proved its convergence under quite general conditions [12, 13]. The SART or OS-SART technology has been widely used in the CT field. Particularly, for the system (2) the SART-type solution is given by [11]

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

(5)

where ${a}_{+n}=\sum _{m=1}^{M}{a}_{m,n}>0$, ${a}_{m+}=\sum _{n=1}^{N}{a}_{m,n}>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 Λ_{n,n}^{+N} = 1/*a*_{+n} and Λ^{M+} ^{M} × ^{M} be a diagonal matrix with Λ_{m,m}^{M+} = 1/*a*_{m+}, (5) can be rewritten as

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

(6)

Note that for all the iteration index *k*, Λ^{M+}, Λ^{+N} and **A** remain unchanged.

Since (4) is a non-linear equation, there is no analytic solution to this problem. Here, we will find a solution in the SART framework summarized in the proceeding subsection. Our strategy is to linearize the equation and approximate the solution successively. Denote the approximate solution for (4) as ${\tilde{\mathbf{f}}}^{k}$ after *k* iterations, we have

$$\tilde{\mathbf{p}}={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}\left({\tilde{\mathbf{f}}}^{k}+\tilde{\mathbf{f}}-{\tilde{\mathbf{f}}}^{k}\right)}={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}\left({\tilde{\mathbf{f}}}^{k}+\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}\right)}={\displaystyle \sum}_{q=1}^{Q}\left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot {e}^{-{\mathbf{A}}^{q}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}}\right),$$

(7)

where $\Delta {\tilde{\mathbf{f}}}^{k}$ is the error image, and “” represents the component-wise multiplication of vectors of the same size. Then, ${e}^{-{\mathbf{A}}^{q}\Delta {\tilde{\mathbf{f}}}^{k}}$in (7) can be expanded in a Taylor series

$${e}^{-{\mathbf{A}}^{q}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}}={\displaystyle \sum}_{t=0}^{\infty}\frac{{\left(-1\right)}^{t}{\left({\mathbf{A}}^{q}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}\right)}^{t}}{t!},$$

(8)

where ${({\mathbf{A}}^{q}\Delta {\tilde{\mathbf{f}}}^{k})}^{t}=({\mathbf{A}}^{q}\Delta {\tilde{\mathbf{f}}}^{k})\odot ({\mathbf{A}}^{q}\Delta {\tilde{\mathbf{f}}}^{k})\cdots .$ Substituting the 1st order approximation of ${e}^{-{\mathbf{A}}^{q}\Delta {\tilde{\mathbf{f}}}^{k}}$ into (7), we have

$$\begin{array}{cc}\hfill \tilde{\mathbf{p}}& \approx {\displaystyle \sum}_{q=1}^{Q}\left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot \left(1-{\mathbf{A}}^{q}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}\right)\right)\hfill \\ \hfill & ={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}-{\displaystyle \sum}_{q=1}^{Q}\left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot \left({\mathbf{A}}^{q}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}\right)\right)\hfill \\ \hfill & ={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}-\left({\displaystyle \sum}_{q=1}^{Q}\left({\mathbf{E}}^{q,k}{\mathbf{A}}^{q}\right)\right)\left(\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k}\right),\hfill \end{array}$$

(9)

where **E**^{q,k} ^{M} × ^{M} be a diagonal matrix with ${E}_{m,m}^{q,k}={e}^{-{\mathbf{A}}_{m}^{q}{\tilde{\mathbf{f}}}^{k}}$. Equation (9) can be rewritten as

$${\tilde{\mathbf{g}}}^{k}={\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k},$$

(10)

where

$$\begin{array}{c}{\tilde{\mathbf{g}}}^{k}={\displaystyle \sum}_{q=1}^{Q}{e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}-\tilde{\mathbf{p}},\\ {\mathbf{B}}^{k}={\displaystyle \sum}_{q=1}^{Q}\left({\mathbf{E}}^{q,k}{\mathbf{A}}^{q}\right).\end{array}$$

(11)

Since the projection errors ${\tilde{\mathbf{g}}}^{k}$ from the involved sources are known, (10) can be understood as a linear system of $\Delta {\tilde{\mathbf{f}}}^{k}$ with a measurement matrix **B**^{k}.

Because (10) has the same structure as that of (2), we can use the SART-type formula (6) to solve for $\Delta {\tilde{\mathbf{f}}}^{k}$. Let Λ^{+N,k} ^{N} × ^{N} be a diagonal matrix with

$${\mathrm{\Lambda}}_{n,n}^{+N,k}=\frac{1}{{\sum}_{m=1}^{M}{b}_{m,n}^{k}},$$

(12)

and Λ^{M+,k} ^{M} × ^{M}be a diagonal matrix with

$${\mathrm{\Lambda}}_{m,m}^{M+,k}=\frac{1}{\sum _{n=1}^{N}{b}_{m,n}^{k}},$$

(13)

we have the SART-type solution for $\Delta {\tilde{\mathbf{f}}}^{k}$

$$\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l}=\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}+{\lambda}^{k,l}\left({\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{B}}^{k}\right)}^{T}{\mathrm{\Lambda}}^{M+,k}({\tilde{\mathbf{g}}}^{k}-{\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1})\right),$$

(14)

where *l* is the iteration index and *λ*^{k,l} the relax parameter. Since

$$\begin{array}{cc}\hfill {\left({\mathbf{B}}^{k}\right)}^{T}& ={\left({\displaystyle \sum}_{q=1}^{Q}\left({\mathbf{E}}^{q,k}{\mathbf{A}}^{q}\right)\right)}^{T}={\displaystyle \sum}_{q=1}^{Q}{\left({\mathbf{E}}^{q,k}{\mathbf{A}}^{q}\right)}^{T}\hfill \\ \hfill & ={\displaystyle \sum}_{q=1}^{Q}\left({\left({\mathbf{A}}^{q}\right)}^{T}{\left({\mathbf{E}}^{q,k}\right)}^{T}\right)={\displaystyle \sum}_{q=1}^{Q}\left({\left({\mathbf{A}}^{q}\right)}^{T}{\mathbf{E}}^{q,k}\right),\hfill \end{array}$$

(15)

Equation (14) can be rewritten as

$$\begin{array}{cc}\hfill \mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l}& =\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}+{\lambda}^{k,l}({\mathrm{\Lambda}}^{+N,k}\left({\displaystyle \sum}_{q=1}^{Q}\left({\left({\mathbf{A}}^{q}\right)}^{T}{\mathbf{E}}^{q,k}\right)\right)\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times {\mathrm{\Lambda}}^{M+,k}\left({\tilde{\mathbf{g}}}^{k}-{\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}\right))\hfill \\ \hfill & =\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}+{\lambda}^{k,l}({\displaystyle \sum}_{q=1}^{Q}(({\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{A}}^{q}\right)}^{T}{\mathbf{E}}^{q,k}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times {\mathrm{\Lambda}}^{M+,k}\left({\tilde{\mathbf{g}}}^{k}-\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}\right)))\hfill \\ \hfill & =\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}+{\lambda}^{k,l}({\displaystyle \sum}_{q=1}^{Q}({\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{A}}^{q}\right)}^{T}{\mathrm{\Lambda}}^{M+,k}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times {\mathbf{E}}^{q,k}({\tilde{\mathbf{g}}}^{k}\mathrm{\hspace{0.17em}\hspace{0.17em}}-\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1})))\hfill \\ \hfill & =\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1}+{\lambda}^{k,l}({\displaystyle \sum}_{q=1}^{Q}({\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{A}}^{q}\right)}^{T}{\mathrm{\Lambda}}^{M+,k}\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\times \left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot ({\tilde{\mathbf{g}}}^{k}\mathrm{\hspace{0.17em}\hspace{0.17em}}-\mathrm{\hspace{0.17em}\hspace{0.17em}}{\mathbf{B}}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,l-1})\right)))\hfill \end{array}$$

(16)

Once we have an approximation solution $\Delta {\tilde{\mathbf{f}}}^{k,L}$ for $\Delta {\tilde{\mathbf{f}}}^{k}$ after *L* iterations,we can update the reconstructed image by

$${\tilde{\mathbf{f}}}^{k+1}={\tilde{\mathbf{f}}}^{k}+\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,L}.$$

(17)

For example, we can choose $\Delta {\tilde{\mathbf{f}}}^{k,0}=0$ as the initial image and set *L* = 1 for one step iteration to approximate $\Delta {\tilde{\mathbf{f}}}^{k}$, which results in a simplified algorithm

$${\tilde{\mathbf{f}}}^{k+1}={\tilde{\mathbf{f}}}^{k}+{\lambda}^{k}\left({\displaystyle \sum}_{q=1}^{Q}\left({\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{A}}^{q}\right)}^{T}{\mathrm{\Lambda}}^{M+,k}\left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot {\tilde{\mathbf{g}}}^{k}\right)\right)\right).$$

(18)

For numerical implementation, our SART-type reconstruction algorithm for image reconstruction from overlapped projections can be summarized in the following pseudocode:

- (S.1.) initialize ${\tilde{\mathbf{f}}}^{k}$ and
*λ*^{k}for*k*: = 0; - (S.2.) compute ${\tilde{\mathbf{g}}}^{k}$, Λ
^{+N,k}and Λ^{M+,k}; - (S.3.) initialize the error image $\Delta {\tilde{\mathbf{f}}}^{k,1}:=0$;
- (S.4.) for
*q*= 1 to*Q*backproject the projection error;- (S.4.1.) weight the projection error ${\tilde{\mathbf{g}}}^{k}$ by ${e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}$;
- (S.4.2.) weight the projection error ${\tilde{\mathbf{g}}}^{k}$ by Λ
^{M+,k}; - (S.4.3.) backproject the weighted projection error towards the
*q*th X-ray source; - (S.4.4.) weight the backprojected error image by Λ
^{+N,k}; - (S.4.5.) add the backprojected image to $\Delta {\tilde{\mathbf{f}}}^{k,1}$ by
$$\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,1}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{:=}\mathrm{\hspace{0.17em}\hspace{0.17em}}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,1}+{\mathrm{\Lambda}}^{+N,k}{\left({\mathbf{A}}^{q}\right)}^{T}{\mathrm{\Lambda}}^{M+,k}\left({e}^{-{\mathbf{A}}^{q}{\tilde{\mathbf{f}}}^{k}}\odot {\tilde{\mathbf{g}}}^{k}\right);$$(19)

- (S.5.) update the current estimated image by

$${\tilde{\mathbf{f}}}^{k+1}\mathrm{\hspace{0.17em}\hspace{0.17em}}\text{:=}\mathrm{\hspace{0.17em}\hspace{0.17em}}{\tilde{\mathbf{f}}}^{k}+{\lambda}^{k}\mathrm{\Delta}{\tilde{\mathbf{f}}}^{k,1};$$

(20)

- (S.6.) set
*k*: =*k*+ 1; - (S.7.) go to (S.2) until the convergence criteria are satisfied;
- (S.8.) scale the reconstructed image ${\tilde{\mathbf{f}}}^{\ast}$ to obtain the final result ${\mathbf{f}}^{\ast}={\tilde{\mathbf{f}}}^{\ast}/s$.

In the above pseudo-code, (S.1.) initializes the iteration index *k*, the relax parameter *λ*, and the initial image ${\tilde{\mathbf{f}}}^{0}$. In our numerical simulation, we always set *λ*^{k} = 1 and ${\tilde{\mathbf{f}}}^{0}=0$. The outer loop (S.2)–(S.7) solves for $\tilde{\mathbf{f}}$ successively. (S.2) precomputes several important intermediate variables to update an reconstructed image in the iteration step *k*. (S.3)-(S.4) computes the current error image for one step iteration according to (16). Because the backprojection operation for different X-ray sources in the inner loop (S.4) has the same structure, it can be implemented by calling a common procedure. (S.5) updates an reconstructed image. (S.6) updates the iteration index. In (S.7), the convergence criteria are checked. The stopping criteria for (S.7) can be the maximum iteration number is reached and/or the relative reconstruction error comes under a predefined threshold [14]. Finally, the reconstructed image **f*** is scaled by dividing with the constant *s*.

The conventional data acquisition is based on the Nyquist sampling theory, which states that for accurate reconstruction of a band-limited signal or image the sampling rate must be at least twice the highest frequency of the signal or image. However, the recently developed CS theory shows that a high-quality signal or image can be reconstructed from far fewer measurements than what is usually required by the Nyquist sampling theorem [7, 8]. In light of the work on solving the linear inverse problems with sparsity constraints by Daubechies et al. [14, 15], we recently adapted the single source SART to reconstruct an image from a limited number of projections subject to a sparsity constraint [16], and demonstrated that the sparsity constraint helped improve the quality of reconstructed images effectively and reduce the number of projections significantly. Using the same strategy described in our previous papers [16, 17], here we use the sparisty constraint to regularize the proposed multi-source SART algorithm. This can be done by adding a soft-threshold filtering step between (S.5) and (S.6) in the pseudo-code given in Section 3.2. Particularly, we have the following pseudo-code segment:

- (SS.5.) perform a soft-threshold filtering operation for ${\tilde{\mathbf{f}}}^{k+1}$;
- (SS.5.1.) compute the sparse transform;
- (SS.5.2.) estimate the optimal threshold;
- (SS.5.3.) perform the soft-threshold filtering;
- (SS.5.4.) perform the inverse sparse transform.

In the above pseudo-code, the sparse transform in (SS.5.1) can be either any invertible lossless compressible transform such as wavelet transform [16] and Fourier transform or uninvertible transforms such as discrete gradient transform (DGT) and discrete difference transform (DDT) [17]. For an uninvertible transform, the inverse sparse transform in (SS.5.4.) is in terms of pseudoinversion as we performed for DGT and DDT [17]. (SS.5.2.) determines an optimal threshold automatically using the projected gradient method for fast convergence [14]. In fact, we can omit (SS.5.2.) and specify any fixed filtering threshold. However, both the convergence speed and final result depend on the choice of the threshold.

To verify the proposed SART-type algorithm for image reconstruction from overlapped projections, we implemented it in Matlab on a PC, with the computationally intensive segments coded in C and linked via the MEX mechanism. As illustrated in Figure 2, we simulated a triple-source fan-beam micro-CT system. In the system, the source XS0 is rotated on a circular scanning locus of radius 120mm. The object was a modified Shepp-Logan phantom in a compact support with a radius of 35mm. We used an equi-distance detector array of length 120mm. The detector was perpendicular to the direction from the origin to the X-ray source XS0 that is 40mm from the system origin. The detector array consisted of 500 elements. On the line through the X-ray source XS0 and parallel to the detector, we put two sources XS1 and XS2 being 25mm apart from XS0 on its right and left sides, respectively. With this triple-source configuration, we simulated single-source (only XS0 was fired), dual-source (XS1 and XS2 fired simultaneously), and triple-source cases (XS0, XS1 and XS2 fired simultaneously).

For each of our selected numbers of projections over a full-scan range, we first equiangularly acquired the corresponding projection dataset based on the aforementioned projection model in the single, dual, and triple source cases, respectively. Then, we reconstructed the images using our algorithm described in Section 3.2. In our simulation, the parameter *λ*^{k} in the SART formula (18) was set to 1.0, and the stopping criterion was defined as reaching the maximum iteration number 5000. Because the Shepp-Logan phantom is a piecewise constant function, its DGT and DDT are sparse. Hence, we also employed the sparsity regularization in terms of total difference minimization [17] and the threshold for filtering was automatically computed using the projected gradient method [14]. Figure 3 shows the reconstructed 256 × 256 images from 9, 11, 13, and 15 projections, respectively. For real-world applications, measurement noise is unavoidable. To test the stability of the proposed algorithm against data noise, we repeated the above reconstructions from projections corrupted by Poisson noise, assuming *I*_{0} = 5 × 10^{4} photons per detector element [18]. The results are in Figure 4, which indicate the stability of the proposed algorithm.

In the CT field, the line integral model along an X-ray path has been widely used in consistency with Beer's law. However, it does not reflect the divergence due to the combination of the finite detector size and the source focal spot. As shown in Figure 1, we have assumed an area model for the X-ray path and normalized it for the multi-source imaging model. Our area model treats the X-ray path as a narrow fan-beam from the X-ray point source to the detector element, and we believe that the area model works better than the line model. Note that the proposed algorithm views the projection procedure as a matrix transform, both the area and line models can be handled by our algorithm. In other words, the proposed algorithm is independent of the imaging model as long as it is linear or can be transformed into a linear one. Additionally, to simplify the derivation and demonstrate our idea, we have assumed that the photon numbers emitting from all the sources to each detector are the same. In fact, these numbers may be different, and can be easily incorporated into our algorithm.

As far as the convergence of the proposed algorithm is concerned, it should converge to the least square solution in the cases of either noise-free and noisy projections. The reason is that the proposed method is in the framework of the general Landweber scheme, whose convergence has been well studied under quite general conditions [12, 13]. When only a small number of projections are available, we can use some sparse constraints to steer the solution to the truth. However, the convergence speed of the current soft-threshold filtering technology is still slow although it has been accelerated using the projected gradient method [14]. In the future, we will employ more advanced techniques for a faster speed, which may include but not limited to optimizing the code, employing parallel computation, and developing new algorithms.

In conclusion, we have developed a SART-type algorithm for image reconstruction from overlapped projections. The algorithm has been verified and demonstrated in the numerical simulation. Our methodology has a potential to support more flexible designs of multi-source CT/micro-CT systems for better contrast and temporal resolution.

This work is partially supported by the NSF/MRI Program (CMMI-0923297) and NIH/NIBIB Grants (EB009275 and EB011785). The authors are grateful for constructive comments from the reviewers.

1. Hounsfield GN. Computerized transverse axial scanning (tomography) .1. Description of system (reprinted from British-Journal-of-Radiology, vol 46, pg. 1016–1022, 1973) *The British Journal of Radiology*. 1995;68(815):H166–H172. [PubMed]

2. Bharkhada D, Wang G, Yu H, Liu H, Plemmons R. Line-source based X-ray tomography. *International Journal of Biomedical Imaging*. 2009;2009:8 pages. Article ID 534516. [PMC free article] [PubMed]

3. Yang L, Lu Y, Wang G. Compressed sensing based image reconstruction from overlapped proejctions. *International Journal of Biomedical Imaging*. 2010;2010:8 pages. Article ID 284073.

4. Zhang J, Yang G, Lee Y, et al. Multiplexing radiography for ultra-fast computed tomography: a feasibility study. *Medical Physics*. 2007;34(6):2527–2527.

5. De Man B, Pelc NJ, Dumoulin C, Bernstein T. Propagation of quantum noise in multiplexed x-ray imaging. In: Hsieh J, Samei E, editors. In: Medical Imaging 2008: Physics of Medical Imaging; 2008; San Diego, Calif, USA. *Proceedings of SPIE*. Paper no. 69131U.

6. Lu Y. Multi-beam field emission x-ray system and its reconstruction algorithm. *Medical Physics*. 2010;37(7):3773–3781. [PubMed]

7. Donoho DL. Compressed sensing. *IEEE Transactions on Information Theory*. 2006;52(4):1289–1306.

8. Candès 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. Ji CG, Yu HY, Wang G. Compressive sensing CT reconstruction with distributed x-ray sources and very few projection views. in draft, 2010.

10. Landweber L. An iteration formula for Fredholm integral equations of the first kind. *American Journal of Mathematics*. 1951;73:615–624.

11. Wang G, Jiang M. Ordered-subset simultaneous algebraic reconstruction techniques (OS-SART) *Journal of X-Ray Science and Technology*. 2004;12(3):169–177.

12. Jiang M, Wang G. Convergence studies on iterative algorithms for image reconstruction. *IEEE Transactions on Medical Imaging*. 2003;22(5):569–579. [PubMed]

13. Jiang M, Wang G. Development of iterative algorithms for image reconstruction. *Journal of X-Ray Science and Technology*. 2002;10(1-2):77–86.

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

16. Yu H, Wang G. SART-type image reconstruction from a limited number of projections with the sparsity constraint. *International Journal of Biomedical Imaging*. 2010;2010:9 pages. Article ID 934847. [PMC free article] [PubMed]

17. Yu HY, Wang G. A soft-threshold filtering approach for reconstruction from a limited number of projections. *Physics in Medicine and Biology*. 2010;55(13):3905–3916. [PMC free article] [PubMed]

18. Yu H, Ye Y, Zhao S, Wang G. Local ROI reconstruction via generalized FBP and BPF algorithms along more flexible curves. *International Journal of Biomedical Imaging*. 2006;2006:7 pages. Article ID 14989. [PMC free article] [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. |