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

**|**HHS Author Manuscripts**|**PMC2923589

Formats

Article sections

- Abstract
- I. Introduction
- II. Local shift-Invariance approximations
- III. The Analytical Variance Prediction
- IV. Fan-beam Geometry
- V. Quadratic Regularization: R0(ρ, Φ)
- VI. Variance Prediction Implementation
- VII. Simulation Results
- VIII. Conclusion and Discussion
- References

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 August 18.

Published in final edited form as:

PMCID: PMC2923589

NIHMSID: NIHMS217507

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

Yingying Zhang-O'Connor: ude.hcimu.scee@zyy; 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.

Accurate predictions of image variances can be useful for reconstruction algorithm analysis and for the design of regularization methods. Computing the predicted variance at every pixel using matrix-based approximations [1] is impractical. Even most recently adopted methods that are based on local discrete Fourier approximations are impractical since they would require a forward and backprojection and two FFT calculations for every pixel, particularly for shift-variant systems like fan-beam tomography. This paper describes new “analytical” approaches to predicting the approximate variance maps of 2D images that are reconstructed by penalized-likelihood estimation with quadratic regularization in fan-beam geometries. The simplest of the proposed analytical approaches requires computation equivalent to one backprojection and some summations, so it is computationally practical even for the data sizes in X-ray CT. Simulation results show that it gives accurate predictions of the variance maps. The parallel-beam geometry is a simple special case of the fan-beam analysis. The analysis is also applicable to 2D PET.

Statistical methods have obtained increasing attention in tomographic image reconstruction due to improved noise and resolution properties. These methods are usually nonlinear and shift-variant. To analyze the statistical characteristics of the reconstructed images, one would like to be able to predict the variances and covariances of estimated pixel values. The variance information provides an uncertainty measure of the reconstructed image and may aid regularization parameter selection.

The existing noise analysis methods can be divided into two categories: iteration based and estimator based. The iteration-based variance predictions are studied in [2], [3] as a function of the iteration number for the maximum-likelihood expectation maximization algorithm based on the “stopping rule” to terminate the iterations before convergence. The estimator-based variance predictions are independent of the particular algorithm and iterations, [1], [4], [5]. Our proposed method falls in the estimator-based category. We will give a brief overview on the existing estimator-based methods and our proposed method.

The estimator-based analysis for the mean and variance proposed in [1] uses the partial derivatives of the cost function and Taylor approximations. The approximations are matrix form and give accurate results. However, the predictions involve the inversion of the Hessian matrices and therefore are computationally expensive. Based on this work, a great deal of effort has been given to simplify these matrix methods [4], [5]. All these methods, that we refer to the DFT approximations, are based on a factorization of the system matrix and circulant approximations to the Hessian matrices to precompute and store a great portion of the calculations. The factorization of the system matrix into geometric and object-dependent portions is specially useful for the shift-varying imaging systems. However, these DFT approximations still require in precomputation one forward and backprojection and two FFT calculations, one for likelihood Hessian and one for penalty Hessian, for each location of interest. Moreover, the expressions are still in matrix form and provide little direct insight into the noise properties.

Our proposed approximations in this paper are still based on the results given in [1] but turn to a very different analysis approach. Instead of working in the discrete space, we use the discrete space Fourier transform (DSFT) and Parseval's theorem to bridge from *the discrete space* to *the continuous space*. Using local shift-invariance approximations and local Fourier analysis, we derive “analytical” closed-form expressions for the local impulse response and local frequency response of the Gram operator and the regularization operator. The final approximations eliminate the need of FFTs for variance predictions, greatly reducing computation for cases where the variance is to be predicted at numerous pixel locations. Furthermore, these approximations provide insight into the resolution and noise properties of the reconstructed images.

Because our analysis is built on the previous work in [1], we briefly repeat its main results here. The goal of transmission image reconstruction is to estimate an attenuation image *μ*[*$\stackrel{\u20d7}{n}$*] from projection data ** Y**, where

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

where ** μ** = (

$$\Phi (\mathit{\text{\mu}},\mathit{\text{Y}})=-L(\mathit{\text{\mu}},\mathit{\text{Y}})+\alpha R(\mathit{\text{\mu}}).$$

(1)

As a concrete example, for transmission tomography under the Poisson noise model, the log-likelihood is

$$L(\mathit{\text{\mu}},\mathit{\text{Y}})=\text{\u2211}_{i}{Y}_{i}log\phantom{\rule{0.2em}{0ex}}({\overline{Y}}_{i}(\mathit{\text{\mu}}))-{\overline{Y}}_{i}(\mathit{\text{\mu}}).$$

(2)

For mono-energetic transmission scans, the measurement means are modeled by

$${\overline{Y}}_{i}(\mathit{\text{\mu}})={b}_{i}{e}^{-{[\mathit{\text{A\mu}}]}_{i}}+{r}_{i},$$

(3)

where ** A** is the system matrix,

We focus on regularization terms of the following form:

$$R(\mathit{\text{\mu}})=\text{\u2211}_{l=1}^{4}{R}_{l}(\mathit{\text{\mu}})$$

(4)

$${R}_{l}(\mathit{\text{\mu}})=\text{\u2211}_{\overrightarrow{n},\overrightarrow{n}-{\overrightarrow{m}}_{l}\in S}{r}_{l}[\overrightarrow{n}]\frac{1}{2}{(\mu [\overrightarrow{n}]-\mu [\overrightarrow{n}-{\overrightarrow{m}}_{l}])}^{2},$$

(5)

where *S* {*$\stackrel{\u20d7}{n}$ _{j}* :

A *N* × *M* lattice with approximately circular FOV. Only the pixels with indices are estimated. In this example, *p*= |*S*| = 12.

The goal of this work is to approximate the covariance matrix Cov{**} efficiently yet accurately, motivated by the problem of designing the regularizer ***R*(** μ**). The proposed prediction methods can be generalized to other log-likelihood terms including 2D emission tomography by modifying

The following approximation to the *p* × *p* covariance matrix of ** was derived in [1]:**

$$\mathit{\text{K}}={({\mathit{\text{A}}}^{\prime}\mathit{\text{WA}}+\alpha \mathbf{\text{P}})}^{-1}{\mathit{\text{A}}}^{\prime}\mathit{\text{WA}}{({\mathit{\text{A}}}^{\prime}\mathit{\text{WA}}+\alpha \mathbf{\text{P}})}^{-1},$$

(6)

where **P** is the Hessian matrix of the roughness penalty. For transmission tomography with the models (1) and (2), ** W** = diag{

$$\text{Cov}\{\widehat{\mu}[{\overrightarrow{n}}_{j}],\widehat{\mu}[{\overrightarrow{n}}_{k}]\}\approx {\mathit{\text{e}}}_{j}^{\text{'}}\mathit{\text{K}}{\mathit{\text{e}}}_{k},$$

(7)

where e* _{j}* denotes the

The matrix method described in (6) and (7) has been used in various applications [5], [9]. Simulation and experimental results have confirmed the accuracy of this covariance approximation in image regions where the non-negativity constraint is usually inactive. However, evaluating (7) is relatively expensive. In this paper, we introduce “continuous space analysis” and use “local stationarity” to develop fast approximations for the variance and covariance of the reconstructed image [*$\stackrel{\u20d7}{n}$*].

The paper is organized as follows. Section II briefly reviews the matrix method and the local shift-invariance approximations. Section III proposes the general analytical approach for the variance approximation. Section IV and V apply this method to fan-beam geometry and quadratic regularization. Section VI and VII analyze the single integral approach used and give simulation results for two types of quadratic regularization, including a comparison of the predicted, DFT-based and empirical standard deviation images. Finally, discussion and conclusions are in Section VIII.

The matrix method described in (6) and (7) is very expensive to compute, even for the variance at a single pixel. To accelerate computation, local shift-invariance approximations are usually used in practice, *e.g.*, [4], [5], [9]–[11].

Let ** M** denote one of the

$$\begin{array}{ll}y[\overrightarrow{n}]\hfill & ={\delta}_{\mathcal{S}}[\overrightarrow{n}]\sum _{{\overrightarrow{n}}^{\prime}\in \mathcal{S}}h(\overrightarrow{n},{\overrightarrow{n}}^{\prime})x[{\overrightarrow{n}}^{\prime}]\hfill \\ & ={\delta}_{\mathcal{S}}[\overrightarrow{n}]\sum _{\overrightarrow{n}\text{'}}h(\overrightarrow{n},{\overrightarrow{n}}^{\prime})x[{\overrightarrow{n}}^{\prime}]{\delta}_{\mathcal{S}}[{\overrightarrow{n}}^{\prime}],\hfill \end{array}$$

(8)

where is an indicator function of *$\stackrel{\u20d7}{n}$* defined as follows:

$${\delta}_{\mathcal{S}}[\overrightarrow{n}]\triangleq \{\begin{array}{cc}1,& \overrightarrow{n}\in \mathcal{S}\\ 0,& \text{otherwise}\end{array}.$$

(9)

In other words, the elements of ** M** correspond to

Near a given location *$\stackrel{\u20d7}{n}$*_{0} of interest, we define a *local impulse response* of ** M** as follows

$${h}_{0}(\overrightarrow{m})\triangleq h({\overrightarrow{n}}_{0}+\lambda \overrightarrow{m},{\overrightarrow{n}}_{0}-(1-\lambda )\overrightarrow{m}){\delta}_{\mathcal{S}}[{\overrightarrow{n}}_{0}+\lambda \overrightarrow{m}]{\delta}_{\mathcal{S}}({\overrightarrow{n}}_{0}-(1-\lambda )\overrightarrow{m}],$$

(10)

where ^{2}, denotes the set of integers. Usually we choose *λ* = 1. However, sometimes we can approximate *h* even for non-integer arguments, in which case *λ* = 1/2 may also be useful [12, p. 870].

We say that *h*(*$\stackrel{\u20d7}{n}$*, *$\stackrel{\u20d7}{n}$*′) is locally shift invariant near *$\stackrel{\u20d7}{n}$*_{0} if *h*(*$\stackrel{\u20d7}{n}$*, *$\stackrel{\u20d7}{n}$*′) ≈ *h*_{0}(*$\stackrel{\u20d7}{n}$* − *$\stackrel{\u20d7}{n}$*′) for *$\stackrel{\u20d7}{n}$* and *$\stackrel{\u20d7}{n}$*′ close to *$\stackrel{\u20d7}{n}$*_{0}. The approximation should be accurate provided *$\stackrel{\u20d7}{n}$* and *$\stackrel{\u20d7}{n}$*′ are “sufficiently close” to *$\stackrel{\u20d7}{n}$*_{0} relative to the width of *h*_{0}. Thus, if the operator ** M** is approximately locally shift invariant near

$$y[\overrightarrow{n}]\approx {\delta}_{\mathcal{S}}[\overrightarrow{n}]\sum _{\overrightarrow{n}\text{'}}{h}_{0}(\overrightarrow{n}-{\overrightarrow{n}}^{\prime})x[{\overrightarrow{n}}^{\prime}]{\delta}_{\mathcal{S}}[{\overrightarrow{n}}^{\prime}],$$

(11)

or equivalently ** y** ≈

Let **T** be the *NM* × *p* matrix such that

$${\text{T}}_{1+n+\mathit{\text{mN}},j}=\{\begin{array}{cc}1,& {\overrightarrow{n}}_{j}=(n,m)\\ 0,& \text{otherwise}\end{array},$$

for *n* = 0, … *N* − 1 and *m* = 0, … *M* − 1. The purpose of **T** is to embed the *p* elements of ** μ** (as shown in Fig 1) back to the 2D

In the spirit of the local shift-invariance approximations presented in Section II, we approximate the covariance matrix in (6) near a given location *$\stackrel{\u20d7}{n}$*_{0} by

$$\begin{array}{l}\mathit{\text{K}}\approx {\mathit{\text{K}}}_{0}\triangleq {\mathit{\text{T}}}^{\prime}{\stackrel{\u2323}{\mathit{\text{K}}}}_{0}\mathit{\text{T}}\\ {\stackrel{\u2323}{\mathit{\text{K}}}}_{0}\triangleq {({\text{F}}_{0}+\alpha {\mathbf{\text{P}}}_{0})}^{-1}{\text{F}}_{0}{({\text{F}}_{0}+\alpha {\mathbf{\text{P}}}_{0})}^{-1},\end{array}$$

where F_{0} and **P**_{0} are the *NM* × *NM* BTTB approximations corresponding to ** A′WA** and

$$\text{Cov}\{\widehat{\mu}[\overrightarrow{n}],\widehat{\mu}[{\overrightarrow{n}}^{\prime}]\}\approx \langle {\stackrel{\u2323}{\mathit{\text{K}}}}_{0}{\mathbf{\text{e}}}_{{\overrightarrow{n}}^{\prime}},{\mathbf{\text{e}}}_{\overrightarrow{n}}\rangle ,$$

(12)

where e* _{$\stackrel{\u20d7}{n}$}* is

Two useful approximations to (12) follow from Parseval's theorem. One option is to interpret the arguments in (11) with a suitable modulo *N* or *M*. In this case, the inner product defined in (12) is in the form of circulant convolution and can be approximated by FFTs:

$$\text{Cov}\left\{\widehat{\mu}[\overrightarrow{n}],\widehat{\mu}[{\overrightarrow{n}}^{\prime}]\right\}\approx \frac{1}{NM}\sum _{\overrightarrow{k}=\overrightarrow{0}}^{\overrightarrow{N}-1}{P}_{d0}[\overrightarrow{k}]{e}^{i{\overrightarrow{\omega}}_{\overrightarrow{k}}\cdot (\overrightarrow{n}-\overrightarrow{n})},$$

(13)

for *$\stackrel{\u20d7}{n}$, $\stackrel{\u20d7}{n}$*′ ≈ *$\stackrel{\u20d7}{n}$*_{0}, where *$\stackrel{\u20d7}{N}$* = (*N, M*), * _{$\stackrel{\u20d7}{k}$}*, = (2

$${P}_{d0}[\overrightarrow{k}]\triangleq \frac{{\Gamma}_{0}[\overrightarrow{k}]}{{({\Gamma}_{0}[\overrightarrow{k}]+\alpha {\Omega}_{0}[\overrightarrow{k}])}^{2}},$$

with

$$\begin{array}{l}{\text{F}}_{0}\approx \mathit{\text{Q}}{\mathbf{\Gamma}}_{0}{\mathit{Q}}^{\prime}\\ {\mathbf{\text{P}}}_{0}\approx \mathit{\text{Q}}{\Omega}_{0}{\mathit{\text{Q}}}^{\prime},\end{array}$$

where ** Q** is the 2D (

$$\begin{array}{ll}\text{Var}\left\{\widehat{\mu}[\overrightarrow{n}]\right\}& \approx \langle {\stackrel{\u2323}{\mathit{\text{K}}}}_{0}{\mathbf{\text{e}}}_{\overrightarrow{n}},{\mathbf{\text{e}}}_{\overrightarrow{n}}\rangle \hfill \\ & \approx \frac{1}{NM}\sum _{\overrightarrow{k}=0}^{\overrightarrow{N}-1}\frac{{\Gamma}_{0}[\overrightarrow{k}]}{{({\Gamma}_{0}[\overrightarrow{k}]+\alpha {\Omega}_{0}[\overrightarrow{k}])}^{2}}.\hfill \end{array}$$

(14)

Generally, evaluating this expression for a single pixel requires a forward and backprojection and two FFTs Computation of this DFT approximation is still expensive for realistic image sizes when the variance must be computed for many or all pixels, particularly for shift-variant systems like fan-beam tomography.

An alternative option is to consider *μ*[*$\stackrel{\u20d7}{n}$*] to be defined over all of ^{2} (two-dimensional integer space) in which case (12) is in the form of ordinary convolution that can be expressed using the discrete-space Fourier transform (DSFT) as follows:

$$\text{Cov}\left\{\widehat{\mu}[\overrightarrow{n}],\widehat{\mu}[{\overrightarrow{n}}^{\prime}]\right\}\approx {\int}_{-\pi}^{\pi}{\int}_{-\pi}^{\pi}{P}_{d0}(\overrightarrow{\omega}){\text{e}}^{i\overrightarrow{\omega}\cdot (\overrightarrow{n}-{\overrightarrow{n}}^{\prime})}\frac{d\overrightarrow{\omega}}{{(2\pi )}^{2}},$$

(15)

where *P _{d}*

$${\text{P}}_{\text{d}0}(\overrightarrow{\omega})\triangleq \frac{{H}_{\text{d}0}(\overrightarrow{\omega})}{{[{H}_{\text{d}0}(\overrightarrow{\omega})+\alpha {R}_{\text{d}0}(\overrightarrow{\omega})]}^{2}},$$

(16)

where *H _{d}*

For regularizer design the standard deviation map of the reconstructed image is one quantity of interest, and our numerical investigation will focus on variance prediction. However, the methodology applies readily to approximate covariances.

Using the DSFT approximation (15), we approximate the variance at pixel *$\stackrel{\u20d7}{n}$*_{0} as follows:

$$\text{Var}\left\{\stackrel{\u02c6}{\mu}[{\overrightarrow{n}}_{0}]\right\}\approx {\int}_{-\pi}^{\pi}{\int}_{-\pi}^{\pi}{\text{P}}_{\text{d}0}(\overrightarrow{\omega})\frac{d\overrightarrow{\omega}}{{(2\pi )}^{2}}.$$

(17)

Let Δ denote the sample spacing in the reconstructed image. By making the change of variable, = (2*πρ*Δ)*$\stackrel{\u20d7}{e}$*_{Φ} where *$\stackrel{\u20d7}{e}$*_{Φ} (cos Φ, sinΦ), we rewrite (17) in terms of polar frequency coordinates (*ρ*, Φ) as follows:

$$\text{Var}\left\{\widehat{\mu}[{\overrightarrow{n}}_{0}]\right\}\approx {\Delta}^{2}{\int}_{0}^{2\pi}{\int}_{0}^{{\rho}_{max}}{P}_{0}(\rho ,\Phi )\rho \phantom{\rule{0.2em}{0ex}}\text{d}\rho \phantom{\rule{0.2em}{0ex}}\text{d}\Phi ,$$

(18)

where ${\mathit{\text{\rho}}}_{max}=\frac{1}{2\Delta}$, and we define

$${P}_{0}(\rho ,\Phi )\triangleq {P}_{\text{d}0}(2\pi \rho \Delta {\overrightarrow{e}}_{\Phi})=\frac{{H}_{0}(\rho ,\Phi )}{{[{H}_{0}(\rho ,\Phi )+\alpha {R}_{0}(\rho ,\Phi )]}^{2}}.$$

(19)

We defined *H*_{0} and *R*_{0} similarly in terms of *H*_{d0} and *R*_{d0}. The variance prediction (18) applies to any 2D geometry. The next section specializes (18) by finding analytical approximations to the local frequency response *H*_{0}(*ρ*,Φ) for the fan-beam geometry.

The following analysis is focused on equiangular fan-beam transmission tomography with an arc detector. However, the method generalizes readily to flat detectors, *i.e.*, equidistant sampling and to parallel-beam geometries. As illustrated in Fig. 2, fan-beam rays are indexed by coordinates (*s*, *β*), where *β* is the angle of the source relative to the *y* axis and *s* is the arc length along the detector. For the case where the detector focal point is at the source position, *γ*(*s*) = *s*/*D*_{sd}, where *γ* is the angle of the ray relative to the source and *D*_{sd} is the source to detector distance. The relation between parallel-beam and fan-beam coordinates is [16]:

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

(20)

$$\phi (s,\beta )=\beta +\gamma (s),$$

(21)

where *D*_{s0} is the source-to-rotation center distance.

To predict variance images in fan-beam transmission tomography using (18), we need to determine the local frequency response *H*_{0}(*ρ*, Φ), or equivalently *H*_{d0}(). We first find the local impulse response.

Consider the 2D object model based on a common basis function *χ*(*$\stackrel{\u20d7}{x}$*) superimposed ona*N* × *M* Cartesian grid as follows:

$$\mu (\overrightarrow{x})=\sum _{\overrightarrow{n}\in \mathcal{S}}\mu [\overrightarrow{n}]\chi \left(\frac{\overrightarrow{x}-{\overrightarrow{x}}_{c}[\overrightarrow{n}]}{\Delta}\right),$$

(22)

where *$\stackrel{\u20d7}{x}$* ^{2} denotes the 2D coordinates of the continuous image space, and *$\stackrel{\u20d7}{x}$ _{c}*[

$$\begin{array}{ll}{\overrightarrow{x}}_{c}[\overrightarrow{n}]\hfill & =(\overrightarrow{n}-{\overrightarrow{\omega}}_{\overrightarrow{x}})\Delta ,\phantom{\rule{0.4em}{0ex}}\overrightarrow{n}\in \mathcal{S}\hfill \\ {\overrightarrow{\omega}}_{\overrightarrow{x}}\hfill & =(\overrightarrow{N}-\overrightarrow{1})/2+{\overrightarrow{c}}_{\overrightarrow{x}},\hfill \end{array}$$

where *$\stackrel{\u20d7}{N}$* = (*N*, *M*) and the user-selectable parameter *$\stackrel{\u20d7}{c}$ _{$\stackrel{\u20d7}{x}$}* denotes an optional spatial offset for the image center.

For simplicity, we assume here that the detector blur *b*(*s*) is locally shift invariant, independent of source position *β*, and acts only along the *s* coordinate. Then we model the mean projections as follows:

$${\overline{y}}_{\beta}[{s}_{k}]=\int b({s}_{k}-{s}^{\prime}){p}_{\phi ({s}^{\prime},\beta )}(r({s}^{\prime}))\phantom{\rule{0.2em}{0ex}}\text{d}{s}^{\prime}$$

(23)

for *s _{k}* = (

$${p}_{\phi}(r)=\int \mu (rcos\phi -\ell sin\phi ,rsin\phi +\ell cos\phi )\phantom{\rule{0.2em}{0ex}}\text{d}\ell .$$

Substituting the basis expansion model in (22) for the object into the measurement model (23) and simplifying leads to the linear model

$${\overline{y}}_{\beta}[{s}_{k}]=\sum _{\overrightarrow{n}\in \mathcal{S}}a({s}_{k},\beta ;\overrightarrow{n})\mu [\overrightarrow{n}],$$

where the fan-beam system matrix elements are samples of the following fan-beam projection of a single basis function centered at *$\stackrel{\u20d7}{x}$ _{c}*[

$$a(s,\beta ;\overrightarrow{n})=\int b(s-{s}^{\prime})\Delta g\left(\frac{r({s}^{\prime})-{r}_{\phi ({s}^{\prime},\beta )}[\overrightarrow{n}]}{\Delta},\phi ({s}^{\prime},\beta )\right)\phantom{\rule{0.2em}{0ex}}\text{d}{s}^{\prime},$$

(24)

where *g*(·, *)* is the Radon transform of *χ*(*$\stackrel{\u20d7}{x}$*) at angle and

$${r}_{\phi}[\overrightarrow{n}]\triangleq {\overrightarrow{x}}_{c}[\overrightarrow{n}]\cdot {\overrightarrow{e}}_{\phi},$$

with $\stackrel{\u20d7}{e}$* _{}* (cos, sin).

Then the elements of the *Gram matrix* are given exactly by

$$\begin{array}{cc}{h}_{\text{d}}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]& =\{\begin{array}{cc}{[{\mathit{\text{A}}}^{\prime}\mathit{\text{WA}}]}_{j{j}^{\prime}},& \overrightarrow{n}={\overrightarrow{n}}_{j}\in \mathcal{S},{\overrightarrow{n}}^{\prime}={\overrightarrow{n}}_{{j}^{\prime}}\in \mathcal{S}\\ 0,& \text{otherwise}\end{array}\\ & ={\stackrel{\u2323}{h}}_{d}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\eta ({\overrightarrow{x}}_{c}[\overrightarrow{n}])\eta ({\overrightarrow{x}}_{c}[\overrightarrow{n}])\hfill \end{array}$$

(25)

where *η*(*$\stackrel{\u20d7}{x}$ _{c}*[

$${\stackrel{\u2323}{h}}_{\text{d}}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\triangleq \sum _{l=1}^{{n}_{\text{A}}}\sum _{k=1}^{{n}_{s}}\omega ({s}_{k},{\beta}_{l})\phantom{\rule{0.2em}{0ex}}a({s}_{k},{\beta}_{l};\overrightarrow{n})\phantom{\rule{0.2em}{0ex}}a({s}_{k},{\beta}_{l};{\overrightarrow{n}}^{\prime})$$

(26)

and *ω*(*s*, *β*) denotes the weighting associated with ** W** and

$${\stackrel{\u2323}{h}}_{d}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\approx \frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{\text{S}}}{\int}_{0}^{2\pi}{\int}_{-\infty}^{\infty}w(s,\beta )\phantom{\rule{0.2em}{0ex}}a(s,\beta ;\overrightarrow{n})\phantom{\rule{0.2em}{0ex}}a(s,\beta ;{\overrightarrow{n}}^{\prime})\text{d}s\phantom{\rule{0.2em}{0ex}}\text{d}\beta ,$$

(27)

where Δ* _{β}* is the source angular sampling interval. Notice that

We develop locally shift-invariant approximations to _{d}[*$\stackrel{\u20d7}{n}$; $\stackrel{\u20d7}{n}$*′] in (27) by reparameterizing variables *s*, *β* using analogs of fan-to-parallel beam rebinning. The following locally shift-invariant approximation to *h*_{d}[*$\stackrel{\u20d7}{n}$; $\stackrel{\u20d7}{n}$*′] is derived in detail in Appendix I:

$${\stackrel{\u2323}{h}}_{\text{d}}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\approx \frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{S}}{\int}_{0}^{2\pi}{w}_{0}(\phi ){\stackrel{\u2323}{h}}_{0}(\Delta (\overrightarrow{n}-{\overrightarrow{n}}^{\prime})\cdot {\overrightarrow{e}}_{\phi},\phi )\phantom{\rule{0.2em}{0ex}}\text{d}\phi ,$$

(28)

where the following 1D autocorrelation is with respect to *r*:

$${\stackrel{\u2323}{h}}_{0}(r,\phi )\triangleq {a}_{0}(r,\phi )\ast {a}_{0}(r,\phi ),$$

and *a*_{0}(*r*, ) is a locally parallel-beam version of the system model defined in (51) (see Appendix I). The angle-dependent weighting *w*_{0}() is associated with pixel *$\stackrel{\u20d7}{n}$*_{0}, accounting for the position-dependent magnification as follows:

$${w}_{0}(\phi )\triangleq \left|{m}_{0}(\phi )\right|w(s({r}_{0}(\phi )),\beta ({r}_{0}(\phi ),(\phi ))$$

(29)

$$\begin{array}{l}{r}_{0}(\phi )\triangleq {r}_{\phi}[{\overrightarrow{n}}_{0}]\\ {m}_{0}(\phi )\triangleq \frac{\partial}{\partial r}s(r){|}_{r={r}_{0}(\phi )}=\frac{{D}_{\text{sd}}/{D}_{\text{s0}}}{\sqrt{1-{({r}_{0}(\phi )/{D}_{\text{s0}})}^{2}}},\end{array}$$

(30)

where *s*(*r*) and *β*(*r*, ) are the inverse of (20) and (21). The shape of the local impulse response (28) is a modification of 1/*r* (cf [17]) with statistically modulated angular weighting. The key property of (28) is that it is locally shift invariant, except for edge effects. This approximation should be reasonably accurate provided that *$\stackrel{\u20d7}{n}$* and *$\stackrel{\u20d7}{n}$*′ are “sufficiently close” to *$\stackrel{\u20d7}{n}$*_{0}, the coordinates of the pixel of interest.

Having found the local impulse response approximation (28), the next step is to find the local frequency response. This requires consideration of the edge effects in (25).

The following local frequency response near a point *$\stackrel{\u20d7}{n}$*_{0} is derived in detail in Appendix II:

$${H}_{0}(\rho ,\Phi )\approx \frac{1}{{\Delta}_{\beta}{\Delta}_{\text{S}}{\Delta}^{2}}{\int}_{0}^{2\pi}{w}_{0}(\phi ){S}_{\phi}(\rho ,\Phi )\phantom{\rule{0.2em}{0ex}}\text{d}\phi ,$$

(31)

where the following function captures both detector response effects and edge effects:

$${S}_{\phi}(\rho ,\Phi )={\left|{A}_{0}(\rho cos(\Phi -\phi ),\phi )\right|}^{2}\cdot {d}_{0}(\phi )sin{\text{c}}^{2}({d}_{0}(\phi )\rho sin(\Phi -\phi )),$$

(32)

*d*_{0}() denotes the length of the chord through *$\stackrel{\u20d7}{n}$*_{0} through the FOV at angle (+*π*/2), and *A*_{0}(*ν*, ) is the 1D FT of *a*_{0}(*r*, ) with respect to *r*.

The local frequency response of the Gram operator in (31) is very accurate However direct implementation of (31) is still computationally demanding. We present here two types of further approximations to simplify (31).

As *d*_{0}() → ∞, one can show that for large |*ρ*|,

$${d}_{0}(\phi )sin{c}^{2}({d}_{0}(\phi )\rho sin(\Phi -\phi ))\to \delta (\rho sin(\Phi -\phi )).$$

Therefore the sinc^{2} term is sharply peaked near Φ = and Φ = ± *π*, so we consider the further simplifying approximation

$${\int}_{0}^{2\pi}{w}_{0}(\phi ){S}_{\phi}(\rho ,\Phi )\phantom{\rule{0.2em}{0ex}}d\phi \approx {w}_{0}(\Phi ){\left|{A}_{0}(\rho ,\Phi )\right|}^{2}{G}_{0}(\rho ,\Phi ),$$

(33)

where

$${G}_{0}(\rho ,\Phi )={\int}_{0}^{2\pi}{d}_{0}(\phi )sin{c}^{2}({d}_{0}(\phi )\rho sin(\Phi -\phi ))\phantom{\rule{0.2em}{0ex}}\text{d}\phi .$$

(34)

Substituting into (31) leads to the “Type I” approximation:

$${H}_{0}(\rho ,\Phi )\approx {H}_{01}(\rho ,\Phi )\triangleq \frac{{w}_{0}(\Phi )}{{\Delta}_{\beta}{\Delta}_{\text{S}}{\Delta}^{2}}{\left|{A}_{0}(\rho ,\Phi )\right|}^{2}{G}_{0}(\rho ,\Phi ).$$

(35)

Although *H*_{01}(*ρ*, Φ) is not separable, we can precompute *w*_{0}(Φ) and tabulate *G*_{0}(*ρ*, Φ) once for all pixels for coarsely sampled Φ. Accurately computing *G*_{0}(*ρ*, Φ) is crucial, therefore finely sampled is necessary in (33).

We can simplify further by using the sifting property of the Dirac impulse:

$${\int}_{0}^{2\pi}{w}_{0}(\phi ){S}_{\phi}\left(\rho ,\Phi \right)\text{d}\phi \approx \frac{2}{\left|\rho \right|}{w}_{0}(\Phi ){\left|{A}_{0}\left(\rho ,\Phi \right)\right|}^{2}.$$

Because typically *A*_{0}(*ν*, ) varies slowly, we also consider the following further approximation:

$${A}_{0}(v,\phi )\approx {A}_{0}(0,\phi ).$$

Combining all the above approximations yields the following separable approximation to the local frequency response:

$${H}_{0}(\rho ,\Phi )\approx {H}_{02}(\rho ,\Phi )\triangleq \frac{2{\left|{A}_{0}(0,\phi )\right|}^{2}}{{\Delta}_{\beta}{\Delta}_{\text{S}}{\Delta}^{2}}\frac{{w}_{0}(\Phi )}{\left|\rho \right|}.$$

(36)

This “Type II” separable form agrees with the familiar FT of
$\frac{1}{\mathit{\text{r}}}$. Figure 3 shows the profiles of two types of local frequency responses for *$\stackrel{\u20d7}{n}$*_{0} at image center in unweighted case. We can see that two profiles agrees with each other closely. The discrepancy is mainly at low frequencies for both the unweighted and weighted cases. The discrepancy is more obvious when *$\stackrel{\u20d7}{n}$*_{0} is off image center.

To evaluate the variance using (18) and (19), we also need the local frequency response of quadratic regularization, *R*_{0}(*ρ*, Φ), [7], [8], [18], [19].

Practical regularization methods are based on the differences between neighboring pixel values. For a discrete-space 2D object *μ*[*$\stackrel{\u20d7}{n}$*], a typical quadratic roughness penalty is given in (4) and (5) for 1st-order differences. The *r _{l}*[

$$R(\mathit{\text{\mu}})=\sum _{\overrightarrow{n}}\sum _{l=1}^{L}{r}_{l,0}\frac{1}{2}{\left(({c}_{l}\ast \ast \mu )[\overrightarrow{n}]\right)}^{2}.$$

The *r _{l}*

$${c}_{l}[\overrightarrow{n}]={\xi}_{l}({\delta}_{2}[\overrightarrow{n}]-{\delta}_{2}[\overrightarrow{n}-{\overrightarrow{m}}_{l}]),$$

or for a 2nd-order difference:

$${c}_{l}[\overrightarrow{n}]={\xi}_{l}({\delta}_{2}[\overrightarrow{n}]-{\delta}_{2}[\overrightarrow{n}-{\overrightarrow{m}}_{l}])\ast \ast {\xi}_{l}({\delta}_{2}[\overrightarrow{n}]-{\delta}_{2}[\overrightarrow{n}-{\overrightarrow{m}}_{l}]),$$

where *ξ _{l}* = ‖

Applying Parseval's theorem, we can rewrite *R*(** μ**) as follows:

$$R(\mathit{\text{\mu}})=\sum _{l=1}^{L}{\int}_{-\pi}^{\pi}{\int}_{-\pi}^{\pi}\frac{1}{2}{r}_{l,0}{\left|{C}_{l}(\overrightarrow{\omega})U(\overrightarrow{\omega})\right|}^{2}\frac{d\overrightarrow{\omega}}{{(2\pi )}^{2}},$$

(37)

where $\mu [\overrightarrow{n}]\phantom{\rule{0.2em}{0ex}}\stackrel{\text{FT}}{\u27f7}\phantom{\rule{0.2em}{0ex}}U(\overrightarrow{\omega})$ and the DSFT of a Λ-order (where Λ ) difference has the following magnitude:

$$\begin{array}{ll}\left|{C}_{l}(\overrightarrow{\omega})\right|\hfill & ={\xi}_{l}^{\wedge}{\left|1-{e}^{-\iota (\overrightarrow{\omega}\cdot {\overrightarrow{m}}_{l})}\right|}^{\wedge}\hfill \\ & ={\xi}_{l}^{\wedge}{2}^{\wedge}{sin}^{\wedge}\left(\frac{1}{2}(\overrightarrow{\omega}\cdot {\overrightarrow{m}}_{l})\right).\hfill \end{array}$$

In the polar coordinates of (19):

$${\left|{C}_{l}(\rho ,\Phi )\right|}^{2}={\left|{C}_{l}(2\pi \rho \Delta {\overrightarrow{e}}_{\Phi})\right|}^{2}={\xi}_{l}^{2\wedge}{4}^{\wedge}{sin}^{2\wedge}(\pi \Delta \rho {\overrightarrow{e}}_{\Phi}\cdot {\overrightarrow{m}}_{l}).$$

(38)

Thus, the *Type I* local frequency response for the regularization operator is

$$\begin{array}{cc}{R}_{0}\left(\rho ,\Phi \right)& ={R}_{01}\left(\rho ,\Phi \right)=\sum _{l=1}^{L}{r}_{l,0}{\left|{C}_{l}\left(\rho ,\Phi \right)\right|}^{2}\\ & =\sum _{l=1}^{L}{r}_{l,0}{\xi}_{l}^{2\wedge}{4}^{\wedge}{sin}^{2\wedge}\left(\pi \Delta \rho {\overrightarrow{e}}_{\Phi}\cdot {\overrightarrow{m}}_{l}\right).\end{array}$$

(39)

Applying the approximation sin(*x*) ≈ *x* to (38) yields:

$$\begin{array}{ll}{\left|{C}_{l}\left(\rho ,\Phi \right)\right|}^{2}& \approx {\xi}_{l}^{2\wedge}{\left({\overrightarrow{m}}_{l}\cdot {\overrightarrow{e}}_{\Phi}\right)}^{2\wedge}{\left(2\pi \Delta \rho \right)}^{2\wedge}\\ & ={\left(2\pi \Delta \rho \right)}^{2\wedge}{\xi}_{l}^{\left(1-2/\upsilon \right)2\wedge}{cos}^{2\wedge}\left(\Phi -{\phi}_{l}\right),\end{array}$$

where the angle between the *l*th neighbors is

$${\phi}_{l}\triangleq {tan}^{-1}\frac{{m}_{l}}{{n}_{l}}.$$

With this simplification, the *Type II* local frequency response of the regularizer is approximately separable in (*ρ*, Φ):

$${R}_{0}\left(\rho ,\Phi \right)\approx {R}_{02}\left(\rho ,\Phi \right)={\left(2\pi \rho \Delta \right)}^{2\wedge}{\stackrel{\sim}{R}}_{0}(\Phi ),$$

(40)

where

$${\stackrel{\sim}{R}}_{0}(\Phi )\triangleq \sum _{l=1}^{L}{\xi}_{l}^{\left(1-2/\upsilon \right)2\wedge}{r}_{l,0}{cos}^{2\wedge}\left(\Phi -{\phi}_{l}\right).$$

This separable form agrees with the familiar FT of the differentiation operation.

Having obtained the approximations to *H*_{0}(*ρ*, Φ), the local frequency response of the Gram operator given in (35) and (36), and to *R*_{0}(*ρ*, Φ), the local frequency response of the regularizer given in (39) and (40), we can discretize the integral (18) again to compute the variance image. There are two variance prediction expressions for fan-beam transmission tomography based on the Type I *H*_{01}(*ρ*, Φ) given in (35) and *R*_{01}(*ρ*, Φ) given in (39), and the Type II *H*_{02}(*ρ*, Φ) given in (36) and *R*_{02}(*ρ*, Φ) given in (40).

The variance prediction using *H*_{01}(*ρ*, Φ) in (35) and *R*_{01}(*ρ*, Φ) in (39) involves a double integral and can be implemented by a double summation. We call this prediction the double integral (DI) approach. The location-dependent weighting function *w*_{0}(Φ) can be precomputed, with the computation equivalent to one back-projection. We can coarsely sample Φ because *P*_{0}(*ρ*, Φ) is fairly smooth along Φ.

The separability properties of *H*_{02}(*ρ*, Φ) in (36) and *R*_{02}(*ρ*, Φ) in (40) enable us to carry out the inner integral over *ρ* analytically. Therefore the double-integral in (18) is reduced to one single integral over Φ. Substituting (36) and (40) into (18) yields the remarkably simple expression:

$$\text{Var}\left\{\widehat{\mu}\left[\overrightarrow{n}\right]\right\}\approx {\int}_{0}^{2\pi}\frac{\zeta /3}{2{\left|{A}_{0}\left(0,\Phi \right)\right|}^{2}{w}_{0}(\Phi )+\alpha 4{\pi}^{2}\zeta {\stackrel{\sim}{R}}_{0}(\Phi )}\phantom{\rule{0.2em}{0ex}}\text{d}\Phi ,$$

(41)

where
$\zeta ={\rho}_{max}^{3}{\Delta}_{\beta}{\Delta}_{\text{S}}{\Delta}^{4}$ is a constant. We call this prediction the single integral (SI) approach. We evaluate this integral using a finite summation, with *w*_{0}(Φ) and _{0}(Φ) precomputed. Therefore, computing (41) is equivalent to one back-projection.

We found empirically that the SI approach (41) gave predictions that could be improved by including a single global scale factor, presumably because the SI approach (41) uses many approximations to achieve its simple form. In particular, we found that the SI method under estimates the variance, presumably because the “Fisher information” implied by Type II local frequency response in (36) is too large for low spatial frequencies. To determine the scale factor, we assumed that the DFT-based approach and the analytical approach should produce equivalent results at the image center Specifically, we used the predicted variance for unweighted least squares estimator with standard quadratic penalty (QPULS) for unit variance data.

For QPULS estimator for unit variance data, the statistical weighting, *w*(*s*, *β*) becomes 1. Consider the standard quadratic penalty with first-order (Λ = 1) differences and second-order neighborhood (L = 4), for which _{1,2,3,4} = 0, *π*/2, *π*/4, and −*π*/4 and *r _{l}*

$${\stackrel{\sim}{R}}_{0}(\Phi )=\sum _{l=1}^{L}\Vert {\overrightarrow{m}}_{l}\Vert =1+\sqrt{2}.$$

(42)

Finally, to determine the scale factor, we computed the ratio of the variance predicted by the DFT approach over that predicted by (41). For the parameters used in our simulations, this factor was (1.13)^{2}. This value would need to recomputed for different system geometries or regularization parameters, but is otherwise patient independent.

To evaluate the performance of the proposed methods, we implemented the variance predictions for fan-beam tomographic images reconstructed by quadratically penalized likelihood algorithm. We simulated 1000 realizations of fan-beam transmission scans using a 256×256 Zubal phantom [22] and a blank scan of 10^{6} counts/detector. The corresponding sinogram size was 444 samples in *s*, spaced by Δ_{s} ≈ 2 mm and 492 source positions over 360°. We simulated the geometry of the GE LightSpeed Pro CT scanner with the source-to-detector distance around 949 mm, the isocenter-to-detector distance 408 mm and Δ = 500/256 mm.

An ellipse support was used for *S*, with *p* = 43892. For simplicity in (34) we used the width of the central profile through the FOV:

$${d}_{0}\left(\phi \right)\approx d\left(\phi \right)\triangleq \frac{2{z}_{1}{z}_{2}}{\sqrt{{z}_{1}^{2}{sin}^{2}\phi +{z}_{2}^{2}{cos}^{2}\phi}},$$

(43)

where *z*_{1} = 244.1 mm and *z*_{2} = 220.7 mm are the semi-major and semi-minor axes of the ellipse. This approximation is exact when *$\stackrel{\u20d7}{n}$*_{0} is at the ellipse center. The approximation to *d*() becomes less accurate as *$\stackrel{\u20d7}{x}$ _{c}*[

For simplicity, we model the detector response^{2} in (23) by a shift-invariant rectangle of width Δ_{s}:

$$b\left(s\right)=\frac{1}{{\Delta}_{\text{S}}}\text{rect}\left(\frac{s}{{\Delta}_{\text{S}}}\right).$$

In the case of a square pixel basis *χ*(*$\stackrel{\u20d7}{x}$*) = rect_{2}(*$\stackrel{\u20d7}{x}$*)^{3}, we have from (51) (see Appendix I)

$${A}_{0}\left(v,\phi \right)=sin\text{c}\left(\frac{{\Delta}_{\text{S}}v}{{m}_{0}\left(\phi \right)}\right){\Delta}^{2}sinc\left(v\Delta cos\phi \right)sinc\left(v\Delta sin\phi \right),$$

(44)

which we substitute into (32). In our simulation, we make the following simplification:

$${m}_{0}\left(\phi \right)\approx {m}_{c}\left(\phi \right)={m}_{c}={D}_{\text{sd}}/{D}_{\text{s}0},$$

where *m _{c}* is the value of

We chose the regularization parameter *α* = 2^{11} to give FWHM = 1.72 pixels, *i.e.*, 3.4 mm, at the center of the image. For each realization, ** was reconstructed using 70 iterations of the convergent incremental optimization transfer algorithm (PL-IOT) [23] with 41 subsets and no nonnegativity constraint. The initial images were the filtered back projection (FBP) images with equivalent spatial resolution, obtained by post-filtering ramp-filtered FBP images with the designed PSF. We computed the sample standard deviation pixel by pixel within the finite support ***S* used in reconstruction. All images and profiles are shown in Hounsfield unit (HU).

Two prediction approaches are investigated here: the DI approximation (18) with Type I *H*_{01} (*ρ*, Φ) in (35) and *R*_{01} (*ρ*, Φ) in (39), and the SI approximation (41) with _{0}(Φ) in (40). The former formula was derived with fewer approximations while the latter formula involves more approximations. The accuracy and computation time are compared below.

We considered two types of regularization: standard and certainty-based [24]. In both cases, we implemented (39) and (40) for regularization with first-order (Λ = 1) differences and second-order neighborhood (*L* = 4). In both cases, the standard deviation images predicted by the DI approach are displayed while both DI and SI predictions are compared in the profile plots.

We first considered a standard quadratic penalty where

$${r}_{l,0}={\kappa}_{c}^{2},$$

and
${\kappa}_{c}^{2}$ is the value of *κ*^{2} [*$\stackrel{\u20d7}{n}$*_{0}] at the image center in (45) below. This choice assures that the resolution of the two studies is matched at the image center.
${\stackrel{\sim}{R}}_{0}(\Phi )=\left(1+\sqrt{2}\right){\kappa}_{c}^{2}$ is a constant for all pixels.

Figure 4 shows the reconstructed images and empirical standard deviation images. The empirical standard deviation image for FBP is also shown. The average of FBP standard deviations is around 2.2 HU, approximately 1.8 times higher than that of PL-IOT, 1.2 HU, illustrating the noise advantage of the statistical reconstruction methods at matched resolution.

Predicted and empirical standard deviation images (in HU) and central profiles for Zabul phantom for PL fan-beam transmission image reconstruction using the standard quadratic penalty:
${\stackrel{\sim}{R}}_{0}=\left(1+\sqrt{2}\right){\kappa}_{c}^{2}$.

Figure 4 also shows the central horizontal and vertical profiles of the standard deviation maps. The analytical, the FFT-based and the empirical standard deviations agree with one another very closely within the object. The largest discrepancy within the object was about 10% in the left lung for unknown reasons.

We next investigate a more complicated regularizer that was designed to achieve nearly uniform spatial resolution [24]. In this case, we used space-varying regularizer:

$${r}_{l,0}={\kappa}^{2}\left[{\overrightarrow{n}}_{0}\right],$$

where

$${\kappa}^{2}\left[{\overrightarrow{n}}_{0}\right]\triangleq \frac{1}{2\pi}{\int}_{0}^{2\pi}w\left(s({r}_{0}\left(\phi \right)\right),\beta \left({r}_{0}\left(\phi \right),\phi )\right)\text{d}\phi .$$

(45)

Here, _{0}(Φ) is still independent of Φ, but varies spatially. Computing the “certainty map” (45) requires a simple back-projection with fan-to-parallel beam rebinning.

Figure 5 shows the reconstructed images, standard deviation images and central horizontal and vertical profiles. In this case the average of FBP standard deviations is around 2.2 HU, approximately 1.8 times higher than that of PL-IOT, 0.7 HU. The analytical, the FFT-based and the empirical standard deviations agree with one another very closely within the object.

Predicted and empirical standard deviation images (in HU) and central profiles for Zubal phantom for PL fan-beam transmission image reconstruction using the certainty-based quadratic penalty.

In both cases, the analytical and the FFT-based predictions are somewhat less accurate near the edge of the finite support used in image reconstruction. This is probably due to the fact that the “local stationarity” approximation is less accurate in this area where the statistical weights *w*(*s*, *β*) can vary rapidly. The approximation (43) may also vary rapidly in our study, so it may be possible to improve accuracy near the edges of the FOV by using *d*_{0}().

In our calculations, we used 123 samples in Φ and 128 samples in *ρ* in (18) to predict a 256×256 standard deviation image. Both DI and SI predictions precompute *w*_{0}(Φ) and *G*_{0}(*ρ*, Φ). The precomputation time for *w*_{0}(Φ) was about 19 seconds on dual Intel Xeon 3.40GHz CPU. The precomputation time for *G*_{0}(*ρ*, Φ) with 492 samples in was 2.3 seconds. The DI prediction requires no scale factor precomputation and the computation time was about 135 seconds. The SI prediction requires the scale factor precomputation that is (1.13)^{2} in our case, and the computation time for prediction was about 0.6 second. In contrast, the FFT-based prediction needed about 374 seconds to compute only one single central profile. As expected, the DI prediction is slightly more accurate than the SI prediction, particularly near edges. The SI prediction matches a bit better with the FFT-based prediction because the scale factor calibration was based on FFT-predicted values. For the purposes of regularization design or noise exploration, we believe that the very fast SI approach is adequate.

Because we only compute two central profiles of the FFT-based prediction in each case, we compute the normalized rms (NRMS) percent errors only for these two central profiles. For ${\kappa}_{c}^{2}$ case, the NRMS percent errors for FFT, DI and SI are 6.6%, 6.8% and 6.6%; for ${\mathit{\text{\kappa}}}_{0}^{2}$ case, the NRMS percent errors for FFT, DI and SI are 6.5%, 6.0% and 8.3%, respectively.

This paper has developed analytical variance approximations for 2D tomography. The double integral (18) with (35) and (39), and the single integral (41) provide fast and accurate variance predictions for a quadratically penalized likelihood estimator in fan-beam tomography. The simplest of the proposed approaches (41) requires one backprojection with some additional summations, which is much less computation than previous FFT-based methods. In fact, using the proposed methods, we can predict the variance map in much less time than it takes to reconstruct a single image. The proposed approximations are especially useful in the case that the variance information is needed for many or all pixels, such as when choosing space-varying regularization parameters [6]. The empirical results from the simulated fan beam CT transmission scans demonstrate that the proposed variance approximations are very accurate. Future work will explore using these predictions for regularization design.

Although we focused on variance prediction, by using (15) we could also easily predict covariances and thus predict the covariance of an ROI whose size is small enough that the local approximations are sufficiently accurate. However, if only a single local autocorrelation function is needed, then the FFT approach is probably easier to use. For analysis of detectability of lesions with unknown locations, autocorrelations at many spatial positions may be needed [10], [25]–[28], in which case the proposed approach based on (15) can save computation. The matrix method described in (6) and (7) is also applicable to other imaging modalities, such as PET and SPECT [1]. Therefore the proposed methods are also readily extended to those imaging modalities, with different considerations of the weighting function.

The proposed analytical variance approximations are only investigated in the case of the shift-invariant detector blur. We can also generalize the analysis to shift-variant detector blur where the local shift-invariance approximation is applicable, *e.g.*, for varifocal collimators in SPECT. In this case, *b*(*s* − *s*′) is replaced by *b*(*s*, *s*′) in (24) and *b*_{0}(*r*, ) in (52) becomes

$${b}_{0}\left(r,\phi \right)\triangleq {m}_{0}\left(\phi \right)b({s}_{0}\left(\phi \right)+{m}_{0}\left(\phi \right)r,{s}_{0}\left(\phi \right)),$$

where

$${\text{s}}_{0}\left(\phi \right)\triangleq {D}_{\text{sd}}arcsin\frac{{r}_{0}\left(\phi \right)}{{D}_{\text{s}0}}.$$

This paper has focused on 2D fan-beam geometry. 3D generalization of these methods can be done by applying the same principles [29]. This paper has also focused on analytical variance approximations for the case of quadratic regularization. An interesting challenge for future work is to generalize the analysis to the case of edge-preserving non-quadratic regularization. The analysis in [30] may be a useful starting point.

This work is supported in part by NIH/NCI grant P01 CA87634 and GE Medical Systems.

Reparameterize variables *s* and *β* in (27) according to the inversion of (20) and (21):

$$\begin{array}{cc}\hfill s\to s(r)& ={D}_{\text{sd}}arcsin(r/{D}_{\text{s}0})\hfill \\ \hfill \beta \to \beta (r,\phi )& =\phi -arcsin(r/{D}_{\text{s}0}).\hfill \end{array}$$

Then the fan-to-parallel beam rebinning of *a*(*s*, *β; $\stackrel{\u20d7}{n}$*) is

$$\begin{array}{ll}a(s\left(r\right),\beta \left(r,\phi \right);\overrightarrow{n})& \approx \int b(s\left(r\right)-s\left({r}^{\prime}\right))\Delta \cdot g\left(\frac{{r}^{\prime}-{r}_{\phi}\left[\overrightarrow{n}\right]}{\Delta},\phi \right)|\dot{s}\left({r}^{\prime}\right)|d{r}^{\prime},\\ & \triangleq a\left(r,\phi ;\overrightarrow{n}\right)\end{array}$$

(46)

because *r*(*s*(*r*′)) = *r*′ and (*s*(*r*′), *β*(*r*, )) ≈ for *r* ≈ *r*′.

A first-order Taylor expansion of *s*(*r*) around *r*′ yields

$$s(r)-s({r}^{\prime})\approx \dot{s}({r}^{\prime})(r-{r}^{\prime}).$$

Substituting into (46), the system matrix elements become

$$a\left(r,\phi ;\overrightarrow{n}\right)\approx \int b\left(\dot{s}\left({r}^{\prime}\right)\left(r-{r}^{\prime}\right)\right)\Delta \cdot g\left(\frac{{r}^{\prime}-{r}_{\phi}\left[\overrightarrow{n}\right]}{\Delta},\phi \right)|\dot{s}\left({r}^{\prime}\right)|d{r}^{\prime}.$$

(47)

Substituting (47) into (27) and changing variables from (*s*, *β*) to (*r*, ) using (20) and (21) yields the local impulse approximation,

$${\stackrel{\u2323}{h}}_{\text{d}}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\approx \frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{\text{S}}}{\int}_{0}^{2\pi}{\int}_{-\infty}^{\infty}\overline{w}(r,\phi )a(r,\phi ;\overrightarrow{n})a(r,\phi ;{\overrightarrow{n}}^{\prime})\cdot \left|J(r)\right|\text{d}r\text{d}\phi =\frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{\text{S}}}{\int}_{0}^{2\pi}\stackrel{\u2323}{w}(\phi ;\overrightarrow{n};{\overrightarrow{n}}^{\prime}){\stackrel{\u2323}{h}}_{\phi}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\text{d}\phi ,$$

(48)

where |*J*(*r*)| = |(*r*)| is the determinant of Jacobian matrix, and

$${\stackrel{\u2323}{h}}_{\phi}\left[\overrightarrow{n};{\overrightarrow{n}}^{\prime}\right]\triangleq {\int}_{-\infty}^{\infty}a\left(r-{r}_{\phi}\left[\overrightarrow{n}\right],\phi ;\overrightarrow{n}\right)a\left(r-{r}_{\phi}\left[{\overrightarrow{n}}^{\prime}\right],\phi ;{\overrightarrow{n}}^{\prime}\right)\text{d}r$$

(49)

$$\begin{array}{cc}\hfill \stackrel{\u2323}{w}\left(\phi ;\overrightarrow{n};{\overrightarrow{n}}^{\prime}\right)& \triangleq \frac{{\int}_{-\infty}^{\infty}\overline{w}\left(r,\phi \right)|J\left(r\right)|a(r-{r}_{\phi}\left[\overrightarrow{n}\right],\phi )a(r-{r}_{\phi}\left[{\overrightarrow{n}}^{\prime}\right],\phi )\text{d}r}{{\int}_{-\infty}^{\infty}a(r-{r}_{\phi}\left[\overrightarrow{n}\right],\phi )a(r-{r}_{\phi}\left[{\overrightarrow{n}}^{\prime}\right],\phi )\text{d}r}\hfill \\ \hfill \overline{w}\left(r,\phi \right)& \triangleq w(s\left(r\right),\beta \left(r,\phi \right)).\hfill \end{array}$$

Let *r*_{0}() *r _{}*[

$$\dot{s}\left({r}^{\prime}\right)\approx \dot{s}\left({r}_{0}\left(\phi \right)\right)\triangleq {m}_{0}\left(\phi \right).$$

(50)

Substituting (50) into (47) and simplifying yields

$$\begin{array}{cc}a\left(r,\phi ;\overrightarrow{n}\right)\hfill & \approx a\left(r,\phi ;{\overrightarrow{n}}_{0}\right)\triangleq {a}_{0}\left(r,\phi \right)\hfill \\ & =\int {b}_{0}\left(r-{r}_{\phi}\left[\overrightarrow{n}\right]-{r}^{\u2033},\phi \right)\Delta g\left(\frac{{r}^{\u2033}}{\Delta},\phi \right)d{r}^{\u2033},\hfill \end{array}$$

(51)

with *r*′ = *r*′ − *r _{}*[

$${b}_{0}\left(r,\phi \right)\triangleq {m}_{0}\left(\phi \right)b\left({m}_{0}\left(\phi \right)r\right).$$

(52)

Therefore, we further simply (48) as follows

$${\stackrel{\u2323}{h}}_{\text{d}}[\overrightarrow{n};{\overrightarrow{n}}^{\prime}]\approx \frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{\text{S}}}{\int}_{0}^{2\pi}{w}_{0}\left(\phi \right){\stackrel{\u2323}{h}}_{0}\left(\Delta \left(\overrightarrow{n}-{\overrightarrow{n}}^{\prime}\right)\cdot {\overrightarrow{e}}_{\phi},\phi \right)d\phi ,$$

(53)

where because (*r*, ) often varies slowly in *r* relative to the typically sharp peak of *a*_{0}(*r*, ) at *r* = 0,

$$\begin{array}{ll}\hfill {\stackrel{\u2323}{h}}_{0}(r,\phi )& \triangleq {a}_{0}(r,\phi )\u2605{a}_{0}(r,\phi ),\hfill \\ \hfill \stackrel{\u2323}{w}(\phi ;\overrightarrow{n};{\overrightarrow{n}}^{\prime})& \approx \frac{{\int}_{-\infty}^{\infty}\overline{w}(r,\phi )|{m}_{0}(\phi )|{a}_{0}^{2}(r-{r}_{\phi}\left[{\overrightarrow{n}}_{0}\right],\phi )\text{d}r}{{\int}_{-\infty}^{\infty}{a}_{0}^{2}(r-{r}_{\phi}\left[{\overrightarrow{n}}_{0}\right],\phi )\text{d}r}\hfill \\ & \approx |{m}_{0}(\phi )|\overline{w}({r}_{0}(\phi ),\phi )\triangleq {w}_{0}(\phi ),\hfill \end{array}$$

(54)

where denotes a 1D autocorrelation with respect to *r*.

The simplest approach to finding the local frequency response would be to take the 2D Fourier transform of the local impulse response in (28), while ignoring the “edge effects” caused by the support functions in (25). We found this approach to yield suboptimal accuracy. One way to consider the edge effects is to use a triangular function with the angular-dependent width:

$${\stackrel{\sim}{h}}_{\phi}\left[\overrightarrow{n};{\overrightarrow{n}}_{0}\right]\triangleq {\stackrel{\sim}{h}}_{0}\left(\Delta \overrightarrow{n}\cdot {\overrightarrow{e}}_{\phi},\phi \right){\eta}_{1}\left(\Delta \overrightarrow{n}\right)$$

where

$${\eta}_{1}\left(\overrightarrow{x}\right)\triangleq \text{tri}\left(\frac{\overrightarrow{x}\cdot {\overrightarrow{e}}_{\phi}^{\perp}}{{d}_{0}\left(\phi \right)}\right),$$

(55)

and *d*_{0}() denotes the length of the chord through *$\stackrel{\u20d7}{n}$*_{0} through the FOV at angle ( + *π*/2). This approach is inspired by circulant approximations for Toeplitz matrices [13], [31], [32], preserving the nonnegative definite property of ** A′WA**. This choice might not be optimal in our application, and further investigation may be beneficial.

One alternative way to consider edge effects is to use the coordinate transformation recommended for analyzing quasi-stationary noise in [12, p. 870] as follows:

$$\begin{array}{ll}{\stackrel{\sim}{h}}_{\phi}\left[\overrightarrow{n};{\overrightarrow{n}}_{0}\right]\hfill & \triangleq {h}_{\phi}\left[{\overrightarrow{n}}_{0}+\overrightarrow{n}/2;{\overrightarrow{n}}_{0}-\overrightarrow{n}/2\right]\\ & ={\stackrel{\u2323}{h}}_{0}(\Delta \overrightarrow{n}\cdot {\overrightarrow{e}}_{\phi},\phi ){\eta}_{2}\left(\Delta \overrightarrow{n}\right),\end{array}$$

where *$\stackrel{\u20d7}{x}$*_{0} *$\stackrel{\u20d7}{x}$ _{c}*[

$${\eta}_{2}(\overrightarrow{x})\triangleq \eta ({\overrightarrow{x}}_{0}+\overrightarrow{x}/2)\eta ({\overrightarrow{x}}_{0}-\overrightarrow{x}/2).$$

(56)

This approach yields a local impulse response that is symmetric in *$\stackrel{\u20d7}{n}$*, thus ensuring that its spectrum is real.

Another alternative is to refer all displacements relative to the point *$\stackrel{\u20d7}{n}$*_{0} as follows:

$$\begin{array}{ll}{\stackrel{\sim}{h}}_{\phi}\left[\overrightarrow{n};{\overrightarrow{n}}_{0}\right]\hfill & \triangleq {h}_{\phi}\left[{\overrightarrow{n}}_{0}+\overrightarrow{n};{\overrightarrow{n}}_{0}\right]\\ & ={\stackrel{\u2323}{h}}_{0}(\Delta \overrightarrow{n}\cdot {\overrightarrow{e}}_{\phi},\phi ){\eta}_{3}\left(\Delta \overrightarrow{n}\right),\end{array}$$

where

$${\eta}_{3}(\overrightarrow{x})\triangleq \eta ({\overrightarrow{x}}_{0}+\overrightarrow{x})\eta ({\overrightarrow{x}}_{0}).$$

(57)

This choice is not symmetric in *$\stackrel{\u20d7}{n}$* but it better corresponds to the local Fourier analysis based on the DFT of *A′WA*e* _{j}*. For simplicity, we could also approximate (57) as follows:

$${\eta}_{4}(\overrightarrow{x})\approx {\eta}_{0}(\overrightarrow{x})=\eta (\overrightarrow{x})\eta ({\overrightarrow{x}}_{0}).$$

(58)

This choice also yields a local impulse response that is symmetric in *$\stackrel{\u20d7}{n}$* provided *η*(*$\stackrel{\u20d7}{x}$*) is symmetric itself.

We focus on *η*_{1} (*$\stackrel{\u20d7}{x}$*) in (55) hereafter because it preserves the property of nonnegative definiteness [33]. Define the following “strip like” function:

$${s}_{\phi}(\overrightarrow{x})\triangleq {\stackrel{\u2323}{h}}_{0}(\overrightarrow{x}\cdot {\overrightarrow{e}}_{\phi},\phi ){\eta}_{1}(\overrightarrow{x}),$$

and ${s}_{\phi}(\overrightarrow{x})\phantom{\rule{0.2em}{0ex}}\underleftrightarrow{2\text{D FT}}\phantom{\rule{0.2em}{0ex}}{S}_{\phi}(u,v)$. Taking the DSFT of (53) yields the following result:

$${H}_{\text{d}0}\left(\overrightarrow{\omega}\right)=\frac{1}{{\Delta}_{\beta}}\frac{1}{{\Delta}_{\text{S}}}{\int}_{0}^{2\pi}{w}_{0}\left(\phi \right){H}_{\phi}\left(\overrightarrow{\omega}\right)d\phi ,$$

(59)

where *H _{}*() is the spectrum of

$$\begin{array}{ll}{H}_{\phi}\left(\overrightarrow{\omega}\right)\hfill & =\sum _{\overrightarrow{n}}{s}_{\phi}\left(\Delta \overrightarrow{n}\right){e}^{-\iota \left(\overrightarrow{\omega}\cdot \overrightarrow{n}\right)}\\ & \approx \frac{1}{{\Delta}^{2}}\int \int {s}_{\phi}\left(\overrightarrow{x}\right){e}^{-\iota \frac{1}{\Delta}\left(\overrightarrow{\omega}\cdot \overrightarrow{x}\right)}\text{d}\overrightarrow{x}\\ & =\frac{1}{{\Delta}^{2}}{S}_{\phi}\left(\frac{\overrightarrow{\omega}}{2\pi \Delta}\right).\end{array}$$

(60)

The 2D FT of *s _{}*(

$${s}_{0}(x,y)={\stackrel{\u2323}{h}}_{0}(x,0)\text{tri}\left(\frac{y}{{d}_{0}(0)}\right)\phantom{\rule{0.2em}{0ex}}\phantom{\rule{0.2em}{0ex}}\underleftrightarrow{2\text{D FT}}\phantom{\rule{0.2em}{0ex}}{S}_{0}(u,\upsilon )={\left|{A}_{0}(u,0)\right|}^{2}{d}_{0}(0)sin{\text{c}}^{2}({d}_{0}(0)\upsilon ),$$

where *A*_{0}(*ν*, ) is associated with the detector response and basis effect, given in (44). By the rotation property of the 2D FT:

$${S}_{\phi}\left(\rho ,\Phi \right)\approx {\left|{A}_{0}\left(\rho cos\left(\Phi -\phi \right),\phi \right)\right|}^{2}\cdot {d}_{0}\left(\phi \right)sin{c}^{2}\left({d}_{0}\left(\phi \right)\rho sin\left(\Phi -\phi \right)\right).$$

Therefore, using (54) and (60), the local frequency response *H*_{0}(*ρ*,Φ) around a point *$\stackrel{\u20d7}{n}$*_{0} is

$${H}_{0}\left(\rho ,\Phi \right)\approx \frac{1}{{\Delta}_{\beta}{\Delta}_{\text{s}}{\Delta}^{2}}{\int}_{0}^{2\pi}{w}_{0}\left(\phi \right){S}_{\phi}\left(\rho ,\Phi \right)\text{d}\phi .$$

(61)

^{1}Throughout the paper we use the subscript “0” to indicate dependence on a given pixel location *$\stackrel{\u20d7}{n}$*_{0}.

^{2}A more accurate model could include detector deadspace and crosstalk effects.

^{3}rect_{2}(*$\stackrel{\u20d7}{x}$*) rect(*x*)rect(*y*) is a 2D square function.

1. Fessler JA. Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography. IEEE Trans Im Proc. 1996 Mar;5(3):493–506. [PubMed]

2. Barrett HH, Wilson DW, Tsui BMW. Noise properties of the EM algorithm: I. Theory. Phys Med Biol. 1994 May;39(5):833–46. [PubMed]

3. Wilson DW, Tsui BMW, Barrett HH. Noise properties of the EM algorithm: II. Monte Carlo simulations. Phys Med Biol. 1994 May;39(5):847–72. [PubMed]

4. Qi J, Leahy RM. A theoretical study of the contrast recovery and variance of MAP reconstructions from PET data. IEEE Trans Med Imag. 1999 Apr;18(4):293–305. [PubMed]

5. Stayman JW, Fessler JA. Efficient calculation of resolution and covariance for fully-3D SPECT. IEEE Trans Med Imag. 2004 Dec;23(12):1543–56. [PubMed]

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

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

8. Fessler JA. Tech Rep 308. Comm and Sign Proc Lab., Dept of EECS, Univ of Michigan; Ann Arbor, MI: Mar, 1997. Spatial resolution properties of penalized weighted least-squares image reconstruction with model mismatch. 48109-2122. Available from http://www.eecs.umich.edu/~fessler.

9. Qi J, Huesman RH. Theoretical study of lesion detectability of MAP reconstruction using computer observers. IEEE Trans Med Imag. 2001 Aug;20(8):815–22. [PubMed]

10. Khurd P, Gindi G. Rapid computation of LROC figures of merit using numerical observers (for SPECT/PET reconstruction) IEEE Trans Nuc Sci. 2005 June;52(3):618–26. [PMC free article] [PubMed]

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

12. Barrett HH, Myers KJ. Foundations of image science. Wiley; New York: 2003.

13. Chan RH, Ng MK. Conjugate gradient methods for Toeplitz systems. SIAM Review. 1996 Sept;38(3):427–82.

14. Fessler JA, Booth SD. Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction. IEEE Trans Im Proc. 1999 May;8(5):688–99. [PubMed]

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

16. Peng H, Stark H. Direct Fourier reconstruction in fan-beam tomography. IEEE Trans Med Imag. 1987 Sept;6(3):209–19. [PubMed]

17. Older JK, Johns PC. Matrix formulation of computed tomogram reconstruction. Phys Med Biol. 1993 Aug;38(8):1051–64.

18. Unser M, Aldroubi A, Eden M. Recursive regularization filters: design, properties, and applications. IEEE Trans Patt Anal Mach Int. 1991 Mar;13(3):272–7.

19. Stark H, Olsen ET. Projection-based image restoration. J Opt Soc Am A. 1992 Nov;9(11):1914–9.

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

21. 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 Sept;22(9):1042–52. [PubMed]

22. Zubal G, Gindi G, Lee M, Harrell C, Smith E. High resolution anthropomorphic phantom for Monte Carlo analysis of internal radiation sources. IEEE Symposium on Computer-Based Medical Systems. 1990:540–7.

23. Ahn S, Fessler JA, Blatt D, Hero AO. Convergent incremental optimization transfer algorithms: Application to tomography. IEEE Trans Med Imag. 2006 Mar;25(3):283–96. [PubMed]

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

25. Khurd P, Gindi G. Fast LROC analysis of Bayesian reconstructed emission tomographic images using model observers. Phys Med Biol. 2005 Apr;50(7):1519–32. [PMC free article] [PubMed]

26. Qi J, Huesman RH. Fast approach to evaluate MAP reconstruction for lesion detection and localization. Proc SPIE 5372, Medical Imaging 2004: Image Perception, Observer Performance, and Technology Assessment. 2004:273–82.

27. Khurd PK, Gindi GR. LROC model observers for emission to-mographic reconstruction. Proc SPIE 5372, Medical Imaging 2004: Image Perception, Observer Performance, and Technology Assessment. 2004:509–20.

28. Yendiki A, Fessler JA. Analysis of unknown-location signal detectability for regularized tomographic image reconstruction. Proc IEEE Intl Symp Biomed Imag. 2006:279–83.

29. Zhang-O'Connor Y, Fessler JA. Fast and accurate variance predictions for 3D cone-beam CT with quadratic regularization. spie. 2007 Submitted.

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

31. Chan TF. An optimal circulant preconditioner for Toeplitz systems. SIAM J Sci Stat Comp. 1988 July;9(4):766–71.

32. Chan TF, Olkin JA. Circulant preconditioners for Toeplitz-block matrices. Numer Algorithms. 1994;6:89–101.

33. Tyrtyshnikov E. Optimal and superoptimal circulant preconditioners. SIAM J Matrix Anal Appl. 1992;13(2):459–73.

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