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

**|**HHS Author Manuscripts**|**PMC2911484

Formats

Article sections

- Abstract
- I. Introduction
- II. Local Impulse Reponse
- III. Frequency Domain Analysis
- IV. The AIMA Approach
- V. Results
- VI. Discussion
- VII. Conclusion and Future Work
- References

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 July 28.

Published in final edited form as:

Published online 2008 October 31. doi: 10.1109/TMI.2008.2007366

PMCID: PMC2911484

NIHMSID: NIHMS217495

EECS Dept., 1301 Beal Ave., The University of Michigan, Ann Arbor, MI 48109-2122

Hugo R. Shi: ude.hcimu@ihsoguh: ude.hcimu.scee@ihsoguh; Jeffrey A. Fessler: ude.hcimu.scee@relssef

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

Statistical methods for tomographic image reconstruction have improved noise and spatial resolution properties that may improve image quality in X-ray CT. Penalized weighted least squares (PWLS) methods using conventional quadratic regularization lead to nonuniform and anisotropic spatial resolution due to interactions between the *weighting*, which is necessary for good noise properties, and the regularizer. Previously, we addressed this problem for parallel-beam emission tomography using matrix algebra methods to design data-dependent, shift-variant regularizers that improve resolution uniformity. This paper develops a fast Angular Integral Mostly Analytical, *AIMA*, regularization design method for 2D fan-beam X-ray CT imaging, for which parallel-beam tomography is a special case. Simulation results demonstrate that the new method for regularization design requires very modest computation and leads to nearly uniform and isotropic spatial resolution in transmission tomography when using quadratic regularization.

Statistical image reconstruction methods for X-ray Computed Tomography (CT) imaging have the potential to reduce patient dose, reduce artifacts from beam hardening and metal objects, and accommodate scanning geometries that are poorly suited for conventional FBP reconstruction. Unregularized image reconstruction leads to excessively noisy images, which necessitates noise control such as a penalized weighted least squares^{1} (PWLS) method, or a penalized-likelihood (PL) method [1].

Although PL and PWLS methods can control noise, interactions between a conventional quadratic roughness penalty and the *weighting* that is explicit in PWLS methods [2], and is implicit in PL methods, results in images with nonuniform and anisotropic spatial resolution, even for idealized shift-invariant imaging systems [3]. Fan-beam geometries used in X-ray CT are shift-variant and contain additional variations in spatial resolution over the field of view (FOV) [4]. The resulting non-uniformities and anisotropy could be avoided in part by using a conventional penalized *unweighted* least-squares (PULS) estimation method. However without the *weighting*, PULS yields poor noise properties (akin to FBP).

Much previous work on regularization design focuses on matrix-based approaches to fit the local impulse response of the estimator to a target impulse response. A shift-variant regularizer based on the aggregate certainty of measurement rays intersecting each pixel was developed that yielded uniform but anisotropic spatial resolution [3]. Stayman parameterized the quadratic regularizer to produce uniform and isotropic spatial resolution [5] and generalized regularization design to other non-Poisson noise models [6]. [7] presents a regularization design method for uniform and isotropic spatial resolution that is not based on an explicit target point spread function (PSF), but rather focuses on circular symmetry and uniformity. Extensions to 3D PET have also been proposed [8], [9].

We previously proposed an analytical approach to regularization design for 2D parallel-beam emission reconstruction that uses continuous space analogs to simplify the regularization design problem [10]. This paper extends [10] by developing a mostly analytical approach to regularization design for fan-beam geometries [11]. §II reviews the concept of a local impulse response (LIR) and discusses the design of regularizers that yield approximately uniform and isotropic spatial resolution. §III uses the frequency domain to gain some insight into the problem and derives a minimization problem to solve for the appropriate regularization coefficients. §IV describes efficient techniques for computing the regularization coefficients. §V presents the results of our simulations using simulated and real X-ray CT data, and §VI analyzes them.

Let ** y** = (

$${\overline{y}}_{\mathit{\text{i}}}(\mathit{x})=E[{y}_{\mathit{\text{i}}}]={b}_{\mathit{\text{i}}}{e}^{-{[\mathit{\text{Ax}}]}_{\mathit{\text{i}}}}+{r}_{\mathit{\text{i}}},$$

(1)

where ** A** is the system matrix,

$$\Psi (\mathit{\text{y}},\mathit{\text{x}})\triangleq {\Vert \ell (\mathit{\text{y}})-\mathit{\text{Ax}}\Vert}_{\mathit{\text{W}}}^{2}+\zeta \text{R}(\mathit{x})$$

(2)

$$\widehat{\mathit{\text{x}}}(\mathit{\text{y}})=\underset{\mathit{\text{x}}}{\text{arg}\phantom{\rule{0.2em}{0ex}}min}\Phi (\mathit{\text{y}},\mathit{x}),$$

(3)

where
$\ell (\mathit{\text{y}})\triangleq -\text{ln}\left(\frac{\mathit{\text{y}}-\mathit{\text{r}}}{\mathit{\text{b}}}\right)$ using element-wise division, *R*(** x**) is a regularizer that controls noise,

The goal of this work is to design *R*(** x**) so that the reconstructed image has approximately uniform and isotropic spatial resolution. The tools developed for this purpose could be modified to achieve other resolution design goals. Analyzing the local impulse response (LIR) is an essential part of assessing the resolution. Here we use the following definition of local impulse response for the

$$\begin{array}{ll}{\mathit{\text{l}}}^{j}\hfill & \triangleq \underset{\epsilon \to 0}{lim}\frac{\widehat{\mathit{\text{x}}}(\overline{\mathit{\text{y}}}({\mathit{\text{x}}}_{\text{true}}+\epsilon {\mathit{\delta}}^{j}))-\widehat{\mathit{\text{x}}}(\overline{\mathit{\text{y}}}({\mathit{\text{x}}}_{\text{true}}))}{\epsilon}\hfill \\ \hfill & ={\nabla}_{\mathit{\text{y}}}\widehat{\mathit{\text{x}}}(\mathit{y}){|}_{\mathit{y}=\overline{\mathit{y}}({\mathit{\text{x}}}_{\text{true}})}{\nabla}_{\mathit{\text{x}}}\overline{\mathit{y}}(\mathit{x}){|}_{\mathit{x}={\mathit{x}}_{\text{true}}}{\mathit{\delta}}^{j},\hfill \end{array}$$

(4)

where *δ** ^{j}* is a Kronecker impulse at the

$${\mathit{\text{l}}}^{j}={\mathit{\text{l}}}^{j}({\mathit{\text{x}}}_{\text{true}},\mathit{\text{R}})={[\mathit{\text{A}}\prime \mathit{\text{W A}}+\zeta \mathit{\text{R}}]}^{-1}\mathit{\text{A}}\prime \mathit{\text{W A}}{\mathit{\delta}}^{j},$$

(5)

where ** R** is the Hessian of the regularizer

As is evident from (5), the local impulse response depends on the regularizer through its Hessian, ** R**. We would like to design

$$\underset{\mathit{\text{R}}}{\text{arg}min}\sum _{\text{j}}{\Vert {\mathit{\text{S}}}^{j}{\mathit{\text{l}}}^{j}({\mathit{\text{x}}}_{\text{true}},\mathit{\text{R}})-{\mathit{\text{l}}}^{0}\Vert}^{2},$$

(6)

where *S** ^{j}* is a shift operator that recenters

This section first reviews the use of discrete Fourier transforms for resolution analysis, leading to a computationally intensive approach to regularization design. We then consider continuous-space analogs that lead to simplified designs.

An analytical expression for *A′W Aδ** ^{j}* is derived in Appendix §B, showing that

Let ** Q** denote an orthonormal discrete Fourier transform matrix centered at pixel

$${\mathit{\text{l}}}^{j}\approx {\left[\mathit{\text{Q}}\prime {\mathbf{\Lambda}}^{j}\mathit{\text{Q}}+\zeta \mathit{\text{Q}}\prime {\mathbf{\Gamma}}^{j}\mathit{\text{Q}}\right]}^{-1}\mathit{\text{Q}}\prime {\mathbf{\Lambda}}^{j}\mathit{\text{Q}}{\mathit{\delta}}^{j}$$

(7)

$$={\mathit{\text{Q}}}^{\prime}\left[\frac{{\mathbf{\Lambda}}^{j}}{{\mathbf{\Lambda}}^{j}+\zeta {\mathbf{\Gamma}}^{j}}\right]\mathit{Q}{\mathit{\delta}}^{j},$$

(8)

where the matrices in the bracketed term are diagonal and the division operation is element-wise on the diagonal. These approximations are accurate only for row and column indices that are “sufficiently close” to pixel *j*. The terms in brackets in (8) correspond to the local frequency response of our estimator. We would like to match these to the frequency response *L*^{0} of a target PSF (that will be discussed in §III-B) as closely as possible, i.e., we want

$${\mathit{\text{L}}}^{j}\triangleq \frac{{\mathbf{\Lambda}}^{j}}{{\mathbf{\Lambda}}^{j}+\zeta {\mathbf{\Gamma}}^{j}}\approx {\mathit{\text{L}}}^{0}.$$

(9)

Based on (9), one might consider a DFT formulation of the regularization design using the following minimization approach:

$${\widehat{\mathbf{\Gamma}}}^{j}=\underset{\Gamma \in \mathcal{T}}{\text{arg}min}\Vert {\mathit{L}}^{0}-\frac{{\mathbf{\Lambda}}^{j}}{{\mathbf{\Lambda}}^{j}+\zeta {\mathbf{\Gamma}}^{j}}\Vert ,$$

(10)

where *L*^{0} is the frequency response of *l*^{0} and denotes the set of possible frequency responses for the regularizer limited by its structure that will be enumerated in §III-B. Alternatively, as a preview of the methods in §IV-A, we can cross multiply the terms in (9) yielding the simpler design criterion:

$${\widehat{\mathbf{\Gamma}}}^{\text{j}}=\underset{\Gamma \in \mathcal{T}}{\text{arg}min}\Vert {\mathit{L}}^{0}({\mathbf{\Lambda}}^{\text{j}}+\zeta {\mathbf{\Gamma}}^{j})-{\mathbf{\Lambda}}^{j}\Vert .$$

(11)

This approach is similar to the formulation in [6]. Because **Λ*** ^{j}* is the DFT of

Prior to [10], regularization design methods were based on discrete Fourier transforms. As shown in Appendix B, the continuous-space analog of **Λ*** ^{j}* in polar coordinates (

$$\begin{array}{ll}{\mathcal{L}}^{\mathit{\text{j}}}(\rho ,\Phi )\hfill & \approx \frac{{w}^{\mathit{\text{j}}}(\Phi )\frac{1}{\left|\rho \right|}}{{w}^{\mathit{\text{j}}}(\Phi )\frac{1}{\left|\rho \right|}+\zeta {R}^{\mathit{\text{j}}}(\rho ,\Phi )}\hfill \\ \hfill & =\frac{{w}^{\mathit{\text{j}}}(\Phi )}{{w}^{\mathit{\text{j}}}(\Phi )+\zeta \left|\rho \right|{R}^{j}(\rho ,\Phi )},\hfill \end{array}$$

(12)

where *R ^{j}*(

Regularizers control roughness by penalizing differences between neighboring pixels. Indexing the image ** x** as a 2D function

$${c}_{l}(n,m)=\frac{1}{\sqrt{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}}(\delta (n,m)-\delta (n-{n}_{\mathit{\text{l}}},m-{m}_{\mathit{\text{l}}})),$$

(13)

where typically (*n _{l}*,

$$R(\mathit{\text{x}})=\sum _{n,m}\sum _{l=1}^{L}\frac{1}{2}{\left|({c}_{\mathit{\text{l}}}\ast \ast x)(n,m)\right|}^{2},$$

(14)

where ** denotes 2D convolution. This conventional regularizer assigns the same weight to the differences between each neighbor. In this paper, we make the regularizer spatially adaptive with the addition of weighting coefficients ${r}_{\mathit{\text{l}}}^{j}$ as follows:

$$R(\mathit{\text{x}})=\sum _{n,m}\sum _{l=1}^{L}{r}_{\mathit{\text{l}}}^{j}\frac{1}{2}{\left|({c}_{\mathit{\text{l}}}\ast \ast x)(n,m)\right|}^{2},$$

(15)

where *j* is the columnized pixel index which is a function of (*n*, *m*). The objective of this paper is to design coefficients
$\left\{{r}_{\mathit{\text{l}}}^{j}\right\}$. To this end, we must analyze the local frequency response *R ^{j}*(

Taking the Fourier transform of (13) yields:

$$\begin{array}{cc}|{C}_{\mathit{\text{l}}}({\omega}_{1},{\omega}_{2}){|}^{2}\hfill & =\frac{1}{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}|1-{e}^{-\mathit{\text{i}}({\omega}_{1}{n}_{\mathit{\text{l}}}+{\omega}_{2}{m}_{\mathit{\text{l}}})}{|}^{2}\\ & =\frac{1}{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}(2-2cos({\omega}_{1},{n}_{\mathit{\text{l}}}+{\omega}_{2},{m}_{\mathit{\text{l}}}).\hfill \end{array}$$

(16)

Combining (16) and (15) yields an accurate expression for the local frequency response of our regularizer that leads to the slower but more accurate Full Integral Iterative NNLS, *FIIN*, method, described in Appendix §A. Next we use an approximation to simplify (16) that will lead to a fast, Angular Integral Mostly Analytical, *AIMA* method.

Using the approximation 2 − 2 cos(*x*) ≈ *x*^{2}, which we will refer to as the *AIMA* approximation, (16) simplifies to

$${\left|{C}_{\mathit{\text{l}}}({\omega}_{1},{\omega}_{2})\right|}^{2}\approx \frac{1}{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}{({\omega}_{1}{n}_{\mathit{\text{l}}}+{\omega}_{2}{m}_{\mathit{\text{l}}})}^{2}.$$

(17)

One can think about isotropy intuitively in polar coordinates as eliminating angular dependence. Therefore we convert (17) to polar frequency coordinates to simplify the analysis. Using frequency and sampling relationships from the DFT, *ω*_{1} = 2*π*Δ* _{x}ρ* cos(Φ) and

$$\begin{array}{ll}{\left|{C}_{\mathit{\text{l}}}({\omega}_{1},{\omega}_{2})\right|}^{2}\hfill & \approx \hfill \\ \hfill & \frac{1}{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}{({n}_{\mathit{\text{l}}}2\pi {\Delta}_{x}\rho cos(\Phi )+{m}_{\mathit{\text{l}}}2\pi {\Delta}_{y}\rho sin(\Phi ))}^{2}\hfill \\ \hfill & =\frac{1}{{n}_{\mathit{\text{l}}}^{2}+{m}_{\mathit{\text{l}}}^{2}}{(2\pi \rho )}^{2}{({n}_{\mathit{\text{l}}}cos(\Phi )+{m}_{\mathit{\text{l}}}sin(\Phi ))}^{2}\hfill \\ \hfill & ={(2\pi \rho )}^{2}{cos}^{2}(\Phi -{\Phi}_{l}),\hfill \end{array}$$

(18)

where
${\Phi}_{l}\triangleq {tan}^{-1}\frac{{m}_{\mathit{\text{l}}}}{{n}_{\mathit{\text{l}}}}$. Taking the Fourier transform of (15) and combining with (18), our final expression for the local frequency response of the regularizer near the *j*th pixel is

$${R}^{\mathit{\text{j}}}(\rho ,\Phi )=\sum _{l=1}^{L}{r}_{\mathit{\text{l}}}^{j}{(2\pi \rho )}^{2}{cos}^{2}(\Phi -{\Phi}_{l}).$$

(19)

For the usual choice of *L* = 4 and for (*n _{l}*,

Substituting (19) into (12) yields a simple expression for the local frequency response of a PWLS estimator. We want to design each regularization coefficient vector ${\mathit{\text{r}}}^{j}=({r}_{1}^{j},\dots ,{r}_{\mathit{\text{L}}}^{j})$ such that our frequency response matches that of the target as closely as possible. We know that the local frequency response associated with PULS is isotropic at the center of the field of view so we select that to be our target.

Using (*n _{l}*,

$${R}_{0}(\rho ,\Phi )={(2\pi \rho )}^{2}$$

(20)

and without the approximation is

$${R}_{0}(\rho ,\Phi )=4-2cos(2\pi \rho cos\Phi )-2sin(2\pi \rho sin\Phi ).$$

(21)

For an unweighted cost function and a parallel beam geometry, the continuous-space frequency response that is analogous to *F*(*A′ Aδ** ^{j}*) is
$\frac{1}{|\mathit{\rho}|}$. As shown in (49) in Appendix §B, for uniform weights (

$${H}^{j}(\rho ,\Phi )=\frac{2}{J({s}^{j}(\Phi ))\left|\rho \right|},$$

(22)

where *J*(*s ^{j}*(Φ)) is the Jacobian for the change of coordinates from parallel-beam to fan-beam geometries as defined in (44), and

$${\mathcal{L}}^{0}\triangleq \frac{\frac{2}{J(0)\left|\rho \right|}}{\frac{2}{J(0)\left|\rho \right|}+\zeta {R}_{0}(\rho ,\Phi )}$$

(23)

$$=\frac{1}{1+\zeta \left|\rho \right|0.5J(0){R}_{0}(\rho ,\Phi )},$$

(24)

where ^{0} is the continuous-space analog of *L*^{0} in (9). Next we will use these continuous-space analogs to solve for the regularization coefficients *r** ^{j}*.

This section describes the Angular Integral Mostly Analytical, *AIMA*, approach. This approach is appropriately named because the continuous space approximations removes all dependence on *ρ* leaving a single integral over the angular coordinate Φ. This also results in a very simple problem that can be solved mostly analytically.

We try to design regularization coefficients *r** ^{j}* to match our designed local frequency response (12) to the target local frequency response (24) as follows:

$$\frac{{w}^{\mathit{\text{j}}}(\Phi )}{{w}^{\mathit{\text{j}}}(\Phi )+\zeta \left|\rho \right|{R}^{\mathit{\text{j}}}(\rho ,\Phi )}\approx \frac{1}{1+\zeta \left|\rho \right|0.5J(0){(2\pi \rho )}^{2}}.$$

(25)

Cross multiplying and simplifying yields

$$\begin{array}{c}{w}^{\mathit{\text{j}}}(\Phi )\zeta \left|\rho \right|0.5J(0){(2\pi \rho )}^{2}\approx \zeta \left|\rho \right|{R}^{\mathit{\text{j}}}(\rho ,\Phi )\hfill \\ \phantom{\rule{1.5em}{0ex}}{w}^{\mathit{\text{j}}}(\Phi )0.5J(0){(2\pi \rho )}^{2}\approx \sum _{l=1}^{L}{r}_{\mathit{\text{l}}}^{j}{(2\pi \rho )}^{2}{cos}^{2}(\Phi -{\Phi}_{l})\hfill \\ \phantom{\rule{0.5em}{0ex}}{\stackrel{\sim}{w}}^{j}(\Phi )\triangleq {w}^{\mathit{\text{j}}}(\Phi )0.5J(0)\approx \sum _{l=1}^{L}{r}_{\mathit{\text{l}}}^{j}{cos}^{2}(\Phi -{\Phi}_{l}).\hfill \end{array}$$

(26)

We design *r** ^{j}* by minimizing the difference between both sides of (26):

$${\mathit{r}}^{\mathit{\text{j}}}=\underset{\mathit{\text{r}}\ge 0}{\text{arg}min}\frac{1}{2\pi}{{\int}_{0}^{2\pi}\left({\stackrel{\sim}{w}}^{j}(\Phi )-\sum _{l=1}^{L}{r}_{l}{cos}^{2}(\Phi -{\Phi}_{l})\right)}^{2}d\Phi .$$

(27)

We solve for *L* coefficients using the above minimization for each pixel *j*. We constrain the *L* coefficients to be non-negative which ensures that the penalty function is convex. This expression also applies to parallel-beam geometries, where *J*(0) = 1.

We can think of the minimization in (27) as a projection of * ^{j}* = (Φ) onto the space spanned by {cos

$${cos}^{2}(\Phi -{\Phi}_{l})=\frac{1}{2}\cdot 1+\frac{cos(2{\Phi}_{l})}{2\sqrt{2}}\left[\sqrt{2}cos(2\Phi )\right]+\frac{sin(2{\Phi}_{l})}{2\sqrt{2}}\left[\sqrt{2}sin(2\Phi )\right].$$

(28)

The three orthonormal basis functions are

$$\begin{array}{l}{p}_{1}(\Phi )=1\\ {p}_{2}(\Phi )=\sqrt{2}cos(2\Phi )\\ {p}_{3}(\Phi )=\sqrt{2}sin(2\Phi ).\end{array}$$

Using (27), we write
${\sum}_{l=1}^{L}{r}_{\mathit{\text{l}}}^{j}{cos}^{2}(\Phi -{\Phi}_{l})={\mathit{\text{PTr}}}^{j}$, where ** P** is an operator whose columns are

Using (28), the minimization problem (27) simplifies to the following expression:

$${\mathit{\text{r}}}^{j}=\underset{\mathit{\text{r}}\ge 0}{\text{arg}min}{\Vert \mathit{\text{Tr}}-{\mathit{\text{b}}}^{j}\Vert}^{2},$$

(29)

where *b*^{j}*P*** ^{j}* (·) and

If there are too many zeros in *r** ^{j}*, there will be zeros in the Hessian, possibly degrading

Now we formulate our problem so that non-negative least squares (NNLS) algorithms will accommodate this new constraint. Let ^{j}*r** ^{j}* +

$$\begin{array}{ll}{\widehat{\mathit{\text{r}}}}^{j}\hfill & =\underset{\mathit{r}\ge 0}{\text{arg}min}\Vert \mathit{\text{Tr}}-({\mathit{\text{b}}}^{j}-\mathit{\text{T}}{\mathit{\epsilon}}^{j})\Vert \hfill \\ \hfill & =\underset{\mathit{r}\ge 0}{\text{arg}min}\Vert \mathit{\text{Tr}}-{\overline{\mathit{b}}}^{j}\Vert ,\hfill \end{array}$$

(30)

where ^{j}*b** ^{j}* −

Using a second-order neighborhood (*L* = 4) we select (*n _{l}*,

$$\begin{array}{c}\mathit{\text{T}}=\frac{1}{2}\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \\ \hfill 1/\sqrt{2}\hfill & \hfill -1/\sqrt{2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1/\sqrt{2}\hfill & \hfill -1/\sqrt{2}\hfill \end{array}\right]\\ {\overline{\mathit{\text{b}}}}^{\text{j}}=\left[\begin{array}{c}{\text{d}}_{1}^{\text{j}}\\ \sqrt{2}{d}_{2}^{j}\\ \sqrt{2}{d}_{3}^{j}\end{array}\right],{d}^{\mathit{\text{j}}}=\left[\begin{array}{c}\frac{1}{2\pi}(1-\alpha ){\int}_{0}^{2\pi}{\stackrel{\sim}{w}}^{j}(\Phi )d\Phi \\ \frac{1}{2\pi}{\int}_{0}^{2\pi}{\stackrel{\sim}{w}}^{j}(\Phi )cos(2\Phi )d\Phi \\ \frac{1}{2\pi}{\int}_{0}^{2\pi}{\stackrel{\sim}{w}}^{j}(\Phi )sin(2\Phi )d\Phi \end{array}\right].\end{array}$$

(31)

Observe that *Tε** ^{j}* = [

The system (31) is under-determined, which is somewhat intuitive since one can obtain approximately isotropic smoothing using only the horizontal and vertical neighbors, or only the diagonal neighbors. For the purposes of regularization design, an under-determined situation is desirable since it allows us to use the “extra” degrees of freedom to ensure non-negativity even when anisotropic regularization is needed.

We could solve the minimization (30) using an iterative *NNLS algorithm* [15, p. 158]. However, using the properties of ** T** and

First octant of quadratic penalty design space showing the four regions where different constraints are active.

- If ${d}_{2}\ge \frac{1}{2}{d}_{1}$ and ${d}_{3}\le \frac{2}{3}{d}_{2}-\frac{1}{3}{d}_{1}$, then ${r}_{1}=\frac{4}{3}\left({d}_{1}+{d}_{2}\right)$,
*r*_{2}=*r*_{3}=*r*_{4}= 0. - If ${d}_{3}\ge \frac{2}{3}{d}_{2}-\frac{1}{3}{d}_{1}$ and ${d}_{3}+{d}_{2}\ge \frac{1}{2}{d}_{1}$, then ${r}_{1}=\frac{8}{5}\left[\frac{1}{2}{d}_{1}+\frac{3}{2}{d}_{2}-{d}_{3}\right]$,
*r*_{2}=*r*_{4}= 0, ${r}_{3}=\frac{12}{5}\left[{d}_{3}-\left(\frac{2}{3}{d}_{2}-\frac{1}{3}{d}_{1}\right)\right]$. - If ${d}_{3}+{d}_{2}\le \frac{1}{2}{d}_{1}$ and ${d}_{2}\ge \frac{1}{4}{d}_{1}$, then there are multiple non-negative
that exactly solve Ψ(*r*) =*r***0**. The minimum-norm solution is*r*_{1}= 4*d*_{2},*r*_{2}= 0,*r*_{3}=*d*_{1}− 2*d*_{2}+ 2*d*_{3}, ${r}_{4}=2\left[\frac{1}{2}{d}_{1}-({d}_{2}+{d}_{3})\right]$. - If ${d}_{2}\le \frac{1}{4}{d}_{1}$, then there are multiple non-negative
that are exact solutions. The natural choice is the minimum-norm*r*given by the pseudo-inverse solution*r*=*r**T*^{†}, where ${r}_{1}=2\left(\frac{1}{4}{d}_{1}+{d}_{2}\right)$, ${r}_{2}=2\left(\frac{1}{4}{d}_{1}-{d}_{2}\right)$, ${r}_{3}=2\left(\frac{1}{4}{d}_{1}+{d}_{3}\right)$, ${r}_{4}=2\left(\frac{1}{4}{d}_{1}-{d}_{3}\right)$.*d*

The analytical solution presented above is for the usual first-order differences (13). For higher-order differences or neighborhoods, it would appear to become increasingly cumbersome to solve (27) analytically, so an iterative NNLS approach may be more appealing. This can still be practical since ** T** is quite small. The analytical solution above is a continuous function of

For practical implementation, we simply discretize the integrals in (31) [16]. This presents interpolation issues in extracting a discretized version of * ^{j}*(Φ), from

We first investigated imaging a phantom consisting of two uniform rings that highlight the effects of non-uniformities and anisotropy [7]. Afterwards, we studied real CT data.

We simulated a 2D 3rd-generation fan-beam CT system using distance-driven forward and backprojections [17]. The rotation center is 40.8cm from the detector, and the source is 94.9cm from the detector. The axis of rotation is at the center of the object. The simulated imaging system has 888 rays per view spaced 1mm apart, and 984 evenly spaced view angles over a 2*π* rotation. The reconstructed images consisted of a 512 × 512 grid of 1 mm pixels. We chose a *ζ* in (24) such that the target PSF has a full width half max (FWHM) of 3.18 mm.

We simulated a noiseless fan-beam sinogram without scatter using a phantom consisting of a background disk and two rings each of thickness 1mm shown in Fig. 2. We generated the noiseless sinogram by taking line-integrals through the analytical phantom with the same system geometry. We used the plug-in weighting *w _{i}* =

Regularization penalty coefficients used in reconstruction of ring phantom. The four images are
${r}_{\mathit{\text{l}}}^{j}$ for *l* = 1, …, 4

We reconstructed images ** $\widehat{x}$** using several methods. We first created an image uniformly blurred by the target PSF. We then reconstructed images using (i) conventional regularization,

Images of right-most ring, Upper-Left: uniformly blurred by target PSF. Upper-Right: reconstructed using conventional regularization. Mid-Left: reconstructed using certainty-based regularization. Mid-Right: reconstructed using *AIMA* regularization, with **...**

Images of left-most ring, Upper-Left: uniformly blurred by target PSF. Upper-Right: reconstructed using conventional regularization. Mid-Left: reconstructed using certainty-based regularization. Mid-Right: reconstructed using *AIMA* regularization, with **...**

Fig. 6 and Fig. 7 show profiles around the two rings of the reconstructed images using the various regularization methods relative to the mean intensity of the rings from our target, PULS reconstruction with conventional regularization. This verifies that *AIMA* and *FIIN* improve resolution uniformity. Here, 0 radians corresponds to the rightmost point of that ring and the angles are measured clockwise.

The second study used a similar imaging geometry using one slice of real CT data from a GE scanner described in [18] with weightings ** W** computed in the same way as [19].

Fig. 8 displays an image reconstructed by PWLS using conventional regularization. The impulse response locations are denoted with crosses and a region is marked which will be displayed in Fig. 9. Fig. 9 show windowed reconstructions using conventional regularization, certainty based regularization, *AIMA* with *α* = 0.1 and *α* = 0, and *FIIN* with *α* = 0.

Reconstructions windowed from 800 to 1200 HUs using, from top to bottom, conventional regularization, certainty based regularization, *AIMA* with *α* = 0.1, *AIMA* with *α* = 0, *FIIN* with *α* = 0.

Fig. 10--1313 show local impulse responses for the five regularization methods at several locations calculated analytically using (5). These figures show from left to right, the target impulse response, and local impulses responses for conventional regularization, certainty based regularization [3], *AIMA* with *α* = 0.1, *AIMA* with *α* = 0, and *FIIN* with *α* = 0. Contour plots of the LIR are displayed below at 0.9, 0.75, 0.5, 0.25, and 0.1 of the maximum value of the target PSF. The LIR becomes more anisotropic near the edge of the FOV. Our Fourier-based regularization scheme compensates for this anisotropy better than the certainty-based approach of [3].

Impulse Responses at (-100,-100). From left to right, target, conventional regularization, certainty based regularization, *AIMA* regularization with *α* = 0.1, *AIMA* regularization with *α* = 0, *FIIN* regularization with *α* = 0.

Impulse Responses at (-130,100). From left to right, target, conventional regularization, certainty based regularization, *AIMA* regularization with *α* = 0.1, *AIMA* regularization with *α* = 0, *FIIN* regularization with *α* = 0.

To quantify the performance of these regularizers on real CT data, we computed the PSF of the regularizers at every 10 pixels within the body. Then we calculated the FWHM of the PSFs at 181 evenly spaced angles, and computed the RMS error between the actual FWHM and the FWHM of the target. Histograms of these errors are displayed in Fig. 14. The mean of the RMS errors for conventional regularization, certainty based regularization, the *AIMA* method with *α* = 0.1, the *AIMA* method with *α* = 0, and the *FIIN* method with *α* = 0 are 2.7, 2.7, 2.3, 2.5, and 2.0 respectively.

Fig. 4 and Fig. 5 provide a qualitative understanding of the spatial resolution properties of various regularization methods for the ring phantom. The phantom used consists of rings of uniform intensity and uniform width, thus images with uniform spatial resolution should have rings with uniform width and uniform intensity. Conventional regularization creates rings with sharper spatial resolution near the edge of the field of view. Certainty based regularization improves ring uniformity but anisotropy remains. The last three reconstructions using *AIMA* with *α* = 0.1, *AIMA* with *α* = 0 and *FIIN* provide rings that look almost identical, and which have a more uniform ring width than conventional, and certainty based regularization.

Fig. 6 and Fig. 7 show the amplitude of the rings traced clockwise. This confirms our initial assessment, that certainty based regularization achieves more uniform spatial resolution than conventional regularization, and that the *AIMA* and *FIIN* methods are all very similar and outperform the previous approaches.

Fig. 9 display a quadrant of windowed reconstructions using real CT data to illustrate the images produced using these regularization design methodologies. However, since we do not know the “truth” for this data, these images provide only a qualitative illustration of the effect of regularization design on spatial resolution. The impulse responses in Fig. 10-Fig. 13 illustrate the effect of regularization design on spatial resolution at various locations. These figures confirm that *AIMA* and *FIIN* methods improved isotropy over conventional and certainty-based regularization. The histogram plot of Fig. 14 and the mean of the PSF errors mentioned previously confirm this. *AIMA* has lower FWHM RMS error both certainty based regularization and conventional regularization. *FIIN* has the lowest FWHM RMS error, however it is much slower than *AIMA*.

The resulting impulse responses from the *AIMA* and *FIIN* methods are not completely isotropic. This may seem to contradict the dramatic improvement these regularization design methods achieved with the ring phantom. However, recall that we are trying to approximate * ^{j}*(Φ) using 3 basis functions (see §IV-A) for the

The analysis of this paper focuses on the resolution nonuniformities cased by statistical weightings, not the resolution variation due to detector response and magnification. A more general regularization design with similar parameterization is discussed in [6]. Using the techniques in this paper to account for these effects is an open problem.

*AIMA* is quite efficient. Computing certainty based regularization takes the time of about 1 backprojection. In *AIMA*, we must first compute * ^{j}*(Φ) which takes the time of about 1 backprojection, and then solve the analytical solution which is very fast.

The results presented indicate that we have a functional and efficient regularization design structure that yields approximately uniform and isotropic spatial resolution for 2D fan-beam systems. The methods also apply to the simpler case of 2D (parallel-beam) PET as described in [10]. In the future, we will extend these methods to 3D PET [9] and cone-beam CT. We will also study the noise properties. Another open question is how to design non-quadratic regularization [20].

Impulse Responses at (150,-120). From left to right, target, conventional regularization, certainty based regularization, *AIMA* regularization with *α* = 0.1, *AIMA* regularization with *α* = 0, *FIIN* regularization with *α* = 0.

We would like to thank GE and J.B. Thibault for giving us the real CT data used in this paper. We would also like to thank the reviewers who suggested that we examine approximation in (17) which led to Appendix §A.

This work was supported in part by the National Institutes of Health/National Cancer Institute (NIH/NCI) under Grant P01 CA87634 and in part by GE Medical Systems.

This appendix describes the slower Double Integral Iterative NNLS, *FIIN*, method that does not use the *AIMA* approximation 2 − 2 cos(*x*) ≈ *x*^{2} in (17) in the expression of the local frequency response of ** R**. This design method is appropriately named because it uses the full integral over all frequency domain polar coordinates instead of just the angular coordinates. The matrices involved in this method are slightly more complex than those in

$${R}^{\mathit{\text{j}}}(\rho ,\Phi )=\sum _{l=1}^{L}\frac{{r}_{\mathit{\text{l}}}^{j}}{{n}_{l}^{2}+{m}_{l}^{2}}(2-2cos({u}_{\mathit{\text{l}}}(\rho ,\Phi ))),$$

where *u _{l}*(

Following (25), we try to match the LIR at each pixel *j* to the target:

$$\begin{array}{l}\frac{{w}^{j}(\Phi )}{{w}^{j}(\Phi )+\beta \left|\rho \right|{R}^{j}(\rho ,\Phi )}\approx \frac{1}{1+\beta \left|\rho \right|0.5J(0){R}_{0}(\rho ,\Phi )}\hfill \\ \phantom{\rule{0.1em}{0ex}}{w}^{j}(\Phi )+\beta \left|\rho \right|{R}^{j}(\rho ,\Phi )\approx {w}^{j}(\Phi )(1+\beta \left|\rho \right|0.5J(0){R}_{0}(\rho ,\Phi ))\hfill \\ \phantom{\rule{3.4em}{0ex}}\beta \left|\rho \right|{R}^{j}(\rho ,\Phi )\approx \beta \left|\rho \right|0.5J(0){w}^{j}(\Phi ){R}_{0}(\rho ,\Phi ))\hfill \\ \phantom{\rule{5em}{0ex}}{R}^{j}(\rho ,\Phi )\approx {\stackrel{\sim}{w}}^{j}(\Phi ){R}_{0}(\rho ,\Phi )),\hfill \end{array}$$

where * ^{j}*(Φ) = 0.5

$${\text{r}}^{j}=\underset{\mathit{r}\ge 0}{\text{arg}min}{\Vert \mathit{\text{Tr}}-{\mathit{\text{b}}}^{j}\Vert}^{2},$$

(32)

where *b** ^{j}* is a vector of inner products between

This appendix considers fan-beam geometries and uses continuous-space analysis to analyze the Fourier transform of the Grammian operator ** A′W A** to simplify the regularization design problem in (12). One can use polar coordinates Φ,

Fig. 15 illustrates the fan-beam geometry that we consider. Let *P* be the rotation isocenter. *D*_{0d} denotes the distance from the point *P* to the detector, *D*_{s0} denotes the distance from the X-ray source to *P*, and *D*_{fs} denotes the distance from the X-ray source to the focal point of the detector arc. Define *D*_{sd} *D*_{0d} + *D*_{s0} to be the total distance from the X-ray source to the center of the detector, and *D*_{fd} *D*_{sd} + *D*_{fs} to be the total distance from the focal point to the center of the detector. This formulation encompasses a variety of system configurations by allowing the detector focal point to differ from the X-ray source location. For flat detectors, *D*_{fs} = ∞. For third-generation X-ray CT systems, *D*_{fs} = 0. For fourth generation X-ray CT systems, *D*_{fs} = −*D*_{s0}.

Let *s* [−*s*_{max}, *s*_{max}] denote the (signed) *arc length* along the detector, where *s* = 0 corresponds to the detector center. Assuming detector elements are equally spaced along the detector, arc length is the natural parameterization. The various angles have the following relationships:

$$\begin{array}{l}\alpha (s)=\frac{s}{{D}_{\text{fd}}}\\ \gamma (s)={tan}^{-1}\left(\frac{\left({D}_{\text{fd}}\right)sin\alpha (s)}{\left({D}_{\text{fd}}\right)cos\alpha (s)-{D}_{\text{fs}}}\right)\end{array}$$

where the two most important cases are

$$\gamma (s)=\{\begin{array}{ll}\alpha (s),\hfill & {D}_{fs}=0\hfill \\ {tan}^{-1}s/{D}_{sd},\hfill & {D}_{fs}=\infty .\hfill \end{array}$$

(33)

The relationship between *s* and *γ* is:

$$s=\{\begin{array}{l}({D}_{\text{fd}})\phantom{\rule{0.3em}{0ex}}\left[\gamma -arcsin\phantom{\rule{0.1em}{0ex}}\left(\frac{{D}_{\text{fs}}}{{D}_{\text{fd}}}sin\gamma \right)\right],\hfill \\ {D}_{\text{sd}}tan\gamma ,{D}_{\text{fs}}=\infty .\hfill \end{array}0\le {D}_{\text{fs}}<\infty $$

(34)

The ray corresponding to detector element *s* and angle *β* is (*s, β*) = {(*x, y*) : *x* cos (*s, β*) + *y* sin (*s, β*) = *r*(*s*)}, where

$$\begin{array}{c}\varphi (s,\beta )\triangleq \beta +\gamma (s)\hfill \\ \phantom{\rule{1em}{0ex}}r(s)\triangleq {D}_{\text{s}0}sin\gamma (s).\hfill \end{array}$$

(35)

The range of *r* depends on the position of the X-ray source and the size of the detector:

$$\left|r(s)\right|\le {r}_{max}\triangleq {D}_{\text{s}0}sin{\gamma}_{max},$$

(36)

where *γ*_{max} *γ*(*s*_{max}) and *s*_{max} is half the total detector arc length. The radius *r*_{max} defines the circular *field of view* of the imaging system. The *fan angle* is 2 *γ*_{max}.

The line-integral projection *p*(*s, β*) of *f* along (*s, β*) is^{2}:

$$\begin{array}{ll}p(s,\beta )\hfill & ={\int}_{\mathcal{L}(s,\beta )}f(x,y)\text{d}\ell \hfill \\ \hfill & =\iint f(x,y)\delta (xcos(\varphi (s,\beta ))+ysin\varphi (s,\beta )-r(s))\text{d}x\phantom{\rule{0.2em}{0ex}}\text{d}y,\hfill \end{array}$$

(37)

for |*s*| ≤ *s*_{max} and 0 ≤ *β* < *β*_{max}. We require *β*_{max} ≥ *π*+2 *γ*_{max} to ensure complete sampling of the FOV (otherwise the impulse response would be highly anisotropic).

The usual inner product for fan-beam projection space is

$$\langle {p}_{1},{p}_{2}\rangle ={\int}_{-{s}_{max}}^{{s}_{max}}{\int}_{0}^{{\beta}_{max}}{p}_{1}\left(s,\beta \right){p}_{2}\left(s,\beta \right)\text{d}s\phantom{\rule{0.2em}{0ex}}\text{d}\beta .$$

Using this inner product in projection space, and the usual ^{2} inner product in image space, the adjoint of is given by

$$\begin{array}{l}(\mathcal{P}\ast p)(x,y)={\int}_{-{s}_{max}}^{{s}_{max}}{\int}_{0}^{{\beta}_{max}}p\left(s,\beta \right)\hfill \\ \delta (xcos\varphi \left(s,\beta \right)+ysin\varphi \left(s,\beta \right)-r(s))p(s,\beta )\phantom{\rule{0.2em}{0ex}}\text{d}s\phantom{\rule{0.2em}{0ex}}\text{d}\beta ,\hfill \end{array}$$

(38)

where *r*(*s*) and (*s, β*) were defined in (35). We will next extend a common derivation of backproject then filter (BPF) tomographic reconstruction to accommodate a user defined weighting, and then change the coordinates from parallel-beam to fan-beam space.

When we analyze the local impulse response, we typically consider recentered local impulse functions. In this derivation we will start with an uncentered local impulse response, and center it at the end by removing a phase term in the frequency domain. The un-centered local impulse response of the Grammian operator is,

$$\mathcal{\text{P}}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j}\left(x,y\right)={\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}\delta ({x}^{j}cos(\phi )+{y}^{j}sin(\phi )-r)w(r,\phi )\delta (xcos(\phi )+ysin(\phi )-r)\mathit{\text{drd}}\phi .$$

Using the sampling property with the first *δ*, define

$${w}^{j}(\phi )\triangleq w(r,\phi )=w({x}^{j}cos(\phi )+{y}^{j}sin(\phi ),\phi ).$$

(39)

We denote the Fourier transform of *δ ^{j}*(

$$\begin{array}{ll}\mathcal{P}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j}(x,y)\hfill & ={\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}{\int}_{-\infty}^{\infty}{F}^{j}(\rho ,\phi ){e}^{\mathit{\text{i}}2\pi \rho r}d\rho \phantom{\rule{0.2em}{0ex}}{w}^{j}(\phi )\delta (xcos\phi +ysin\phi -r)\mathit{\text{drd}}\phi \hfill \\ \hfill & ={\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}{F}^{j}(\rho ,\Phi ){w}^{j}(\phi ){e}^{\mathit{\text{i}}2\pi \rho (xcos\phi +ysin\phi )}d\rho d\phi .\hfill \end{array}$$

(40)

This is nearly the inverse Fourier transform in signed polar coordinates except for a *ρ* scale factor. Dividing by *ρ*,

$$\mathcal{P}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j}(x,y)={\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}\frac{{w}^{j}(\phi )F({\delta}^{j}(x,y))}{\rho}{e}^{\mathit{\text{i}}2\pi \rho (xcos\phi +ysin\phi )}d\rho d\phi .$$

The local impulse response is recentered, *h ^{j}*(

$${H}^{j}\left(\rho ,\Phi \right)=\frac{{w}^{j}(\Phi )}{\rho}.$$

(41)

The natural indexing for fan-beam data is arc-length *s* and angle *β*. The analogs to *r, ϕ* for parallel-beam systems are *r*(*s*), *ϕ*(*s, β*) as defined in (35). The weighting expression *w*(*r, ϕ*) is indexed as *w*(*s, β*) in fan-beam coordinates. We start by looking at the fan-beam projection in terms of the analogs to parallel-beam coordinates,

$$\left(\mathcal{P}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j}\right)\left(x,y\right)={\int}_{-{\beta}_{max}}^{{\beta}_{max}}{\int}_{-{s}_{max}}^{{s}_{max}}w\left(s,\beta \right)\delta ({x}^{j}cos(\varphi (s,\beta ))+{y}^{j}sin(\varphi (s,\beta ))-r(s))\delta (xcos(\varphi (s,\beta ))+ysin(\varphi (s,\beta ))-r(s))\mathit{\text{dsd}}\beta .$$

We use the change of variables

$$r\prime =r(s)={D}_{\text{s}0}sin\gamma (s)$$

(42)

$$\varphi \prime =\varphi (s,\beta )=\beta +\gamma (s)$$

(43)

as defined in (35) which has the corresponding Jacobian determinant

$$J(s)=\left|{D}_{\text{s}0}cos\gamma (s)\left|\right|\dot{\gamma}(s)\right|.$$

(44)

Then,

$$\begin{array}{l}(\mathcal{P}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j})(x,y)={\int}_{0}^{2\pi}{\int}_{-{r}_{max}}^{{r}_{max}}w(r\prime ,\varphi \prime )\delta ({x}^{\mathit{\text{j}}}cos\varphi \prime +{y}^{\mathit{\text{j}}}sin\varphi \prime -r\prime )\delta (xcos\varphi \prime +ysin\varphi \prime -r\prime )\frac{1}{J(s(r\prime ))}\mathit{\text{drd}}\varphi .\end{array}$$

In this expression,

$$\begin{array}{c}w(r,\varphi )\triangleq w(s(r),\beta (r,\varphi ))\hfill \\ \phantom{\rule{1.2em}{0ex}}s(r)={\gamma}^{-1}(arcsin(r/{D}_{\text{s}0}))\hfill \\ \phantom{\rule{0.2em}{0ex}}\beta (r,\varphi )=\varphi -arcsin(r/{D}_{\text{s}0}).\hfill \end{array}$$

Using the sampling property of the first *δ* as in the parallel beam case,

$${w}^{j}(\varphi \prime )\triangleq w({x}^{j}cos\varphi \prime +{y}^{j}sin\varphi \prime ,\varphi \prime )$$

(45)

$${s}^{j}(\varphi \prime )\triangleq s(r\prime ){\mid}_{r\prime ={x}^{j}cos\varphi \prime +{y}^{j}sin\varphi \prime}.$$

(46)

Again, let *F ^{j}*(

$$\begin{array}{l}(\mathcal{P}\ast \mathcal{\text{WP}}{\mathit{\delta}}^{j})(x,y)\\ ={\int}_{-{r}_{max}}^{{r}_{max}}{\int}_{0}^{2\pi}{\int}_{-\infty}^{\infty}\frac{{F}^{j}(\rho ,\Phi ){w}^{j}(\varphi \prime )}{J({s}^{j}(\varphi \prime ))}{e}^{\mathit{\text{i}}2\pi pr\prime}\delta (xcos\varphi \prime +ysin\varphi \prime -r\prime )d\rho d\varphi \prime dr\prime \\ ={\int}_{0}^{2\pi}{\int}_{-\infty}^{\infty}\frac{{F}^{j}(\rho ,\Phi ){w}^{j}(\varphi )}{J({s}^{j}(\varphi \prime ))}{e}^{\mathit{\text{i}}2\pi \rho (xcos\varphi \prime +ysin\varphi \prime )}d\rho d\varphi \prime \\ ={\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}\frac{{F}^{j}(\rho ,\Phi ){w}^{j}(\varphi \prime )}{J({s}^{j}(\varphi \prime ))}{e}^{\mathit{\text{i}}2\pi \rho (xcos\varphi \prime +ysin\varphi \prime )}d\rho d\varphi \prime \\ +{\int}_{0}^{\pi}{\int}_{-\infty}^{\infty}\frac{{F}^{j}(\rho ,\Phi ){w}^{j}(\varphi \u2033)}{J({s}^{j}(\varphi \u2033))}{e}^{\mathit{\text{i}}2\pi \rho (xcos\varphi \u2033+ysin\varphi \u2033)}d\rho d\varphi \u2033,\end{array}$$

where ″ = *′* + *π*. This is similar to the parallel-beam derivation except that we have two integrals. We can convert each integral into the inverse Fourier transform as we did in the parallel-beam case and strip out the phase term *F ^{j}*(

$${H}^{j}(\rho ,\Phi )=\frac{{\stackrel{\sim}{w}}^{j}(\Phi )}{\left|\rho \right|},$$

(47)

where,

$${\stackrel{\sim}{w}}^{j}(\Phi )=\frac{1}{J({s}^{j}(\Phi ))}\left[{w}^{j}(\varphi ){|}_{\varphi =\Phi}+{w}^{j}(\varphi ){|}_{\varphi =\Phi +\pi}\right].$$

(48)

Because of the absolute value function in (44), *J*(*s*) is invariant to the *π* phase shift. For the case where we have uniform weighting, *w*(*s, β*) = 1 and therefore ** W** =

$${H}^{j}(\rho ,\Phi )=\frac{2}{J({s}^{j}(\Phi ))\left|\rho \right|}.$$

(49)

We use this equation in the calculation of a target local impulse response (24).

In this appendix, we show that $\sqrt{{d}_{2}^{2}+{d}_{3}^{2}}\le {d}_{1}$, as claimed below (31) in §IV-C. Squaring the integrals in (31), we have:

$$\begin{array}{l}{d}_{1}^{2}=\frac{1}{\pi}{\int}_{0}^{\pi}{\int}_{0}^{\pi}{w}^{\mathit{\text{j}}}(X){w}^{\mathit{\text{j}}}(Y)\mathit{\text{dX dY}}\\ {d}_{2}^{2}=\frac{1}{\pi}{\int}_{0}^{\pi}{\int}_{0}^{\pi}cos(2X){w}^{\mathit{\text{j}}}(X)cos(2Y){w}^{\mathit{\text{j}}}(Y)\mathit{\text{dX dY}}\\ {d}_{3}^{2}=\frac{1}{\pi}{\int}_{0}^{\pi}{\int}_{0}^{\pi}sin(2X){w}^{\mathit{\text{j}}}(X)sin(2Y){w}^{\mathit{\text{j}}}(Y)\mathit{\text{dX dY}}.\end{array}$$

In particular,

$$\begin{array}{ll}{d}_{2}^{2}+{d}_{3}^{2}\hfill & =\frac{1}{\pi}{\int}_{0}^{\pi}{\int}_{0}^{\pi}{w}^{\mathit{\text{j}}}(X){w}^{\mathit{\text{j}}}(Y)\left[cos(2X)cos(2Y)+sin(2X)sin(2Y)\right]\mathit{\text{dX dY}}\hfill \\ \hfill & =\frac{1}{\pi}{\int}_{0}^{\pi}{\int}_{0}^{\pi}{w}^{\mathit{\text{j}}}(X){w}^{\mathit{\text{j}}}(Y)cos(2X-2Y)\mathit{\text{dX dY}}.\hfill \end{array}$$

Thus,
${d}_{2}^{2}+{d}_{3}^{2}\le {d}_{1}^{2}$ since *w ^{j}*(Φ) ≥ 0 for all Φ and cos(2

^{1}In this paper we write PWLS and penalized unweighted least squares (PULS) because the techniques used in this paper can be generalized for non-quadratic regularization. However, the scope of this paper is quadratic regularization and all regularizers used here are quadratic.

^{2}Practically speaking, the integral should be restricted to the field of view:
$\sqrt{{x}^{2}+{y}^{2}}\le {r}_{max}$, but this restriction would complicate analysis by introducing a shift variance into the problem, so we ignore it.

1. Fessler JA. Statistical image reconstruction methods for transmission tomography. In: Sonka M, Fitzpatrick J Michael, editors. Handbook of Medical Imaging, Volume 2. Medical Image Processing and Analysis. SPIE; Bellingham: 2000. pp. 1–70.

2. Fessler JA. Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Trans Med Imag. 1994 Jun;13(2):290–300. [PubMed]

3. Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs. IEEE Trans Im Proc. 1996 Sep;5(9):1346–58. [PubMed]

4. Joseph PM. The influence of modulation transfer function shape on CT image quality. Radiology. 1982 Oct;145(1):179–85. [PubMed]

5. Stayman JW, Fessler JA. Regularization for uniform spatial resolution properties in penalized-likelihood image reconstruction. IEEE Trans Med Imag. 2000 Jun;19(6):601–15. [PubMed]

6. Stayman JW, Fessler JA. Compensation for nonuniform resolution using penalized-likelihood reconstruction in space-variant imaging systems. IEEE Trans Med Imag. 2004 Mar;23(3):269–84. [PubMed]

7. Nuyts J, Fessler JA. A penalized-likelihood image reconstruction method for emission tomography, compared to post-smoothed maximum-likelihood with matched spatial resolution. IEEE Trans Med Imag. 2003 Sep;22(9):1042–52. [PubMed]

8. Qi J, Leahy RM. Resolution and noise properties of MAP reconstruction for fully 3D PET. IEEE Trans Med Imag. 2000 May;19(5):493–506. [PubMed]

9. Shi H, Fessler JA. Quadratic regularization design for 3d cylindrical PET. Proc IEEE Nuc Sci Symp Med Im Conf. 2005;4:2301–5.

10. Fessler JA. Analytical approach to regularization design for isotropic spatial resolution. Proc IEEE Nuc Sci Symp Med Im Conf. 2003;3:2022–6.

11. Shi H, Fessler JA. Quadratic regularization design for fan beam transmission tomography. Proc SPIE 5747, Medical Imaging 2005: Image Proc. 2005:2023–33.

12. Fessler JA. Technical Report 297. Comm and Sign Proc Lab, Dept of EECS, Univ of Michigan; Ann Arbor, MI: Aug, 1995. Resolution properties of regularized image reconstruction methods; pp. 48109–2122.

13. Li Q, Asma E, Qi J, Bading JR, Leahy RM. Accurate estimation of the Fisher information matrix for the PET image reconstruction problem. IEEE Trans Med Imag. 2004 Sep;23(9):1057–64. [PubMed]

14. Clinthorne NH, Pan TS, Chiao PC, Rogers WL, Stamos JA. Preconditioning methods for improved convergence rates in iterative reconstructions. IEEE Trans Med Imag. 1993 Mar;12(1):78–83. [PubMed]

15. Lawson CL, Hanson RJ. Solving least squares problems. Prentice-Hall; 1974.

16. Fessler JA. Matlab tomography toolbox. 2004. Available from http://www.eecs.umich.edu/~fessler.

17. De Man B, Basu S. Distance-driven projection and backprojection. Proc IEEE Nuc Sci Symp Med Im Conf. 2002;3:1477–80.

18. Thibault JB, Sauer K, Bouman C, Hsieh J. A three-dimensional statistical approach to improved image quality for multi-slice helical CT. Med Phys. 2007 Nov;34(11):4526–44. [PubMed]

19. Thibault JB, Bouman CA, Sauer KD, Hsieh J. A recursive filter for noise reduction in statistical iterative tomographic imaging. Proc SPIE 6065, Computational Imaging IV. 2006:60650X.

20. Ahn S, Leahy RM. Spatial resolution properties of nonquadratically regularized image reconstruction for PET. Proc IEEE Intl Symp Biomed Imag. 2006:287–90.

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