PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2010 August 18.
Published in final edited form as:
PMCID: PMC2923589
NIHMSID: NIHMS217507

Fast predictions of variance images for fan-beam transmission tomography with quadratic regularization

Yingying Zhang-O'Connor, Student member, IEEE and Jeffrey A. Fessler, Fellow, IEEE

Abstract

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.

Index Terms: variance approximation, local discrete Fourier analysis, fan-beam tomography, penalized-likelihood image reconstruction

I. Introduction

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 μ[n] from projection data Y, where n is a vector denoting the 2D image pixel location. We focus here on penalized-likelihood estimators obtained by minimizing a cost function as follows:

μ^=argminμΦ(μ,Y),

where μ = (μ[n1], …, μ[np]) [set membership] Rp (p-dimensional real space). The cost function includes a negative log-likelihood term and a regulation term:

Φ(μ,Y)=L(μ,Y)+αR(μ).
(1)

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

L(μ,Y)=iYilog(Y¯i(μ))Y¯i(μ).
(2)

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

Y¯i(μ)=bie[]i+ri,
(3)

where A is the system matrix, bi denotes the blank scan, and ri denotes the additive contribution of scatter to the ith ray.

We focus on regularization terms of the following form:

R(μ)=l=14Rl(μ)
(4)
Rl(μ)=n,nmlSrl[n]12(μ[n]μ[nml])2,
(5)

where S [triangle, equals] {nj : j = 1,…,p} denotes the subset of the N × M lattice that is estimated and ml [set membership] {(1, 0), (0, 1), (1, 1), (−1, 1)}. The roughness penalty (4) involves the horizontal, vertical and diagonal neighbors and allows for the possibility of using regularization coefficients {rl[n]} that vary both with spatial location and direction [6], [7]. In general 1 ≤ pNM and p < NM because the physical field of view (FOV) is a subset of the lattice, see Fig. 1.

Fig. 1
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{[mu]} 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 W in (??) below.

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

K=(AWA+αP)1AWA(AWA+αP)1,
(6)

where P is the Hessian matrix of the roughness penalty. For transmission tomography with the models (1) and (2), W = diag{Yi}. In practice Yi is unknown, so we plug in Yi as an approximation [8]. The covariance between pixels [mu][nk] and [mu][nj] can be approximated using (6) as follows:

Cov{μ^[nj],μ^[nk]}ej'Kek,
(7)

where ej denotes the jth unit column vector of length p.

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 [mu][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.

II. Local shift-Invariance approximations

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 p × p matrices in (6), such as A′WA or P, or inverses or sums thereof. Then a matrix-vector operation y = Mx can be expressed equivalently as

y[n]=δS[n]nSh(n,n)x[n]=δS[n]n'h(n,n)x[n]δS[n],
(8)

where An external file that holds a picture, illustration, etc.
Object name is nihms217507ig1.jpg is an indicator function of n defined as follows:

δS[n]{1,nS0,otherwise.
(9)

In other words, the elements of M correspond to Mkj = h(nk, nj).

Near a given location n0 of interest, we define a local impulse response of M as follows1:

h0(m)h(n0+λm,n0(1λ)m)δS[n0+λm]δS(n0(1λ)m],
(10)

where m [set membership] Z2, Z 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(n, n′) is locally shift invariant near n0 if h(n, n′) ≈ h0(nn′) for n and n′ close to n0. The approximation should be accurate provided n and n′ are “sufficiently close” to n0 relative to the width of h0. Thus, if the operator M is approximately locally shift invariant near n0, then we can approximate the superposition sum (8) by (almost) a convolution sum:

y[n]δS[n]n'h0(nn)x[n]δS[n],
(11)

or equivalently yM0x, where the p × p matrix M0 is defined by [M0]kj = h0(nknj). The expression (11) is almost a convolution sum, except for the “edge conditions” of the indicator functions. If the point n0 is not “too close” to the boundaries of the support mask An external file that holds a picture, illustration, etc.
Object name is nihms217507ig2.jpg, then we may able to disregard the indicator functions and treat the expression as a convolution.

Let T be the NM × p matrix such that

T1+n+mN,j={1,nj=(n,m)0,otherwise,

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 N × M lattice. Then M0 = T[M with breve]0T, where ([M with breve]0)n, n = h0(nn′) is an NM × NM matrix that is block Toeplitz with Toeplitz blocks (BTTB). Thus we can make a circulant approximation to [M with breve]0, [13]. Such approximations are often reasonably accurate except near the edges of the FOV, where the differences between “Toeplitz” and “circulant” end conditions are largest. The local impulse response (10) and the corresponding circulant approximation are two key tools for analysis.

III. The Analytical Variance Prediction

In the spirit of the local shift-invariance approximations presented in Section II, we approximate the covariance matrix in (6) near a given location n0 by

KK0TK0TK0(F0+αP0)1F0(F0+αP0)1,

where F0 and P0 are the NM × NM BTTB approximations corresponding to A′WA and P, respectively. Then we approximate the covariance between pixels [mu][n] and [mu][n′] in (7) by the following inner product:

Cov{μ^[n],μ^[n]}K0en,en,
(12)

where en is nth unit vector of length NM.

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:

Cov{μ^[n],μ^[n]}1NMk=0N1Pd0[k]eiωk(nn),
(13)

for n, n′ ≈ n0, where N = (N, M), [omega with right arrow above]k, = (2πk1/N, 2πk2/M) and

Pd0[k]Γ0[k](Γ0[k]+αΩ0[k])2,

with

F0QΓ0QP0QΩ0Q,

where Q is the 2D (N, M)-point orthonormal DFT matrix. The diagonal matrices Γ0 and Ω0 have diagonal elements Γ0 [k] and Ω0 [k] that are the 2D DFT coefficients of the local impulse response of A′WA and P near n0, respectively. This DFT/FFT approximation has been used in [4], [14], [15] to predict variance at a single pixel:

Var{μ^[n]}K0en,en1NMk=0N1Γ0[k](Γ0[k]+αΩ0[k])2.
(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 μ[n] to be defined over all of Z2 (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:

Cov{μ^[n],μ^[n]}ππππPd0(ω)eiω(nn)dω(2π)2,
(15)

where Pd0([omega with right arrow above]) is the local spectrum of K0, given as follows:

Pd0(ω)Hd0(ω)[Hd0(ω)+αRd0(ω)]2,
(16)

where Hd0([omega with right arrow above]) is the local frequency response of the Gram matrix A′WA and Rd0([omega with right arrow above]) is the local frequency response of P near n0. To our knowledge, this paper is the first to use (15) to develop analytical variance approximations as a faster alternative to the DFT approach (14).

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 n0 as follows:

Var{μˆ[n0]}ππππPd0(ω)dω(2π)2.
(17)

Let Δ denote the sample spacing in the reconstructed image. By making the change of variable, [omega with right arrow above] = (2πρΔ)eΦ where eΦ [triangle, equals] (cos Φ, sinΦ), we rewrite (17) in terms of polar frequency coordinates (ρ, Φ) as follows:

Var{μ^[n0]}Δ202π0ρmaxP0(ρ,Φ)ρdρdΦ,
(18)

where ρmax=12Δ, and we define

P0(ρ,Φ)Pd0(2πρΔeΦ)=H0(ρ,Φ)[H0(ρ,Φ)+αR0(ρ,Φ)]2.
(19)

We defined H0 and R0 similarly in terms of Hd0 and Rd0. The variance prediction (18) applies to any 2D geometry. The next section specializes (18) by finding analytical approximations to the local frequency response H0(ρ,Φ) for the fan-beam geometry.

IV. 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/Dsd, where γ is the angle of the ray relative to the source and Dsd is the source to detector distance. The relation between parallel-beam and fan-beam coordinates is [16]:

Fig. 2
Angular coordinates in fan-beam geometry.
r(s)=Ds0sinγ(s)
(20)
φ(s,β)=β+γ(s),
(21)

where Ds0 is the source-to-rotation center distance.

A. Local Impulse Response

To predict variance images in fan-beam transmission tomography using (18), we need to determine the local frequency response H0(ρ, Φ), or equivalently Hd0([omega with right arrow above]). We first find the local impulse response.

Consider the 2D object model based on a common basis function χ(x) superimposed onaN × M Cartesian grid as follows:

μ(x)=nSμ[n]χ(xxc[n]Δ),
(22)

where x [set membership] R2 denotes the 2D coordinates of the continuous image space, and xc[n] denotes the center of the basis function. Typically

xc[n]=(nωx)Δ,nSωx=(N1)/2+cx,

where N = (N, M) and the user-selectable parameter cx 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:

y¯β[sk]=b(sks)pφ(s,β)(r(s))ds
(23)

for sk = (kwSS and k = 1,…, ns, where ΔS is the sample spacing in s, wS is defined akin to wx, and p[var phi](r) is the 2D Radon transform of μ(x):

pφ(r)=μ(rcosφsinφ,rsinφ+cosφ)d.

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

y¯β[sk]=nSa(sk,β;n)μ[n],

where the fan-beam system matrix elements are samples of the following fan-beam projection of a single basis function centered at xc[n]:

a(s,β;n)=b(ss)Δg(r(s)rφ(s,β)[n]Δ,φ(s,β))ds,
(24)

where g(·, [var phi]) is the Radon transform of χ(x) at angle [var phi] and

rφ[n]xc[n]eφ,

with e[var phi] [triangle, equals] (cos[var phi], sin[var phi]).

Then the elements of the Gram matrix are given exactly by

hd[n;n]={[AWA]jj,n=njS,n=njS0,otherwise=hd[n;n]η(xc[n])η(xc[n])
(25)

where η(xc[n]) [triangle, equals] 1{n[set membership] S},

hd[n;n]l=1nAk=1nsω(sk,βl)a(sk,βl;n)a(sk,βl;n)
(26)

and ω(s, β) denotes the weighting associated with W and nA denotes the number of samples of the source position β. To simplify (25), we first use an integral to approximate the summation in (26) as follows:

hd[n;n]1Δβ1ΔS02πw(s,β)a(s,β;n)a(s,β;n)dsdβ,
(27)

where Δβ is the source angular sampling interval. Notice that hd[n; n′] in (27) is not shift invariant.

We develop locally shift-invariant approximations to hd[n; n′] in (27) by reparameterizing variables s, β using analogs of fan-to-parallel beam rebinning. The following locally shift-invariant approximation to hd[n; n′] is derived in detail in Appendix I:

hd[n;n]1Δβ1ΔS02πw0(φ)h0(Δ(nn)eφ,φ)dφ,
(28)

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

h0(r,φ)a0(r,φ)a0(r,φ),

and a0(r, [var phi]) is a locally parallel-beam version of the system model defined in (51) (see Appendix I). The angle-dependent weighting w0([var phi]) is associated with pixel n0, accounting for the position-dependent magnification as follows:

w0(φ)|m0(φ)|w(s(r0(φ)),β(r0(φ),(φ))
(29)
r0(φ)rφ[n0]m0(φ)rs(r)|r=r0(φ)=Dsd/Ds01(r0(φ)/Ds0)2,
(30)

where s(r) and β(r, [var phi]) 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 n and n′ are “sufficiently close” to n0, the coordinates of the pixel of interest.

B. Local frequency response

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 n0 is derived in detail in Appendix II:

H0(ρ,Φ)1ΔβΔSΔ202πw0(φ)Sφ(ρ,Φ)dφ,
(31)

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

Sφ(ρ,Φ)=|A0(ρcos(Φφ),φ)|2d0(φ)sinc2(d0(φ)ρsin(Φφ)),
(32)

d0([var phi]) denotes the length of the chord through n0 through the FOV at angle ([var phi]+π/2), and A0(ν, [var phi]) is the 1D FT of a0(r, [var phi]) with respect to r.

C. Further approximations of local frequency response

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

1) Type I non-separable form

As d0([var phi]) → ∞, one can show that for large |ρ|,

d0(φ)sinc2(d0(φ)ρsin(Φφ))δ(ρsin(Φφ)).

Therefore the sinc2 term is sharply peaked near Φ = [var phi] and Φ = [var phi] ± π, so we consider the further simplifying approximation

02πw0(φ)Sφ(ρ,Φ)dφw0(Φ)|A0(ρ,Φ)|2G0(ρ,Φ),
(33)

where

G0(ρ,Φ)=02πd0(φ)sinc2(d0(φ)ρsin(Φφ))dφ.
(34)

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

H0(ρ,Φ)H01(ρ,Φ)w0(Φ)ΔβΔSΔ2|A0(ρ,Φ)|2G0(ρ,Φ).
(35)

Although H01(ρ, Φ) is not separable, we can precompute w0(Φ) and tabulate G0(ρ, Φ) once for all pixels for coarsely sampled Φ. Accurately computing G0(ρ, Φ) is crucial, therefore finely sampled [var phi] is necessary in (33).

2) Type II separable form

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

02πw0(φ)Sφ(ρ,Φ)dφ2|ρ|w0(Φ)|A0(ρ,Φ)|2.

Because typically A0(ν, [var phi]) varies slowly, we also consider the following further approximation:

A0(v,φ)A0(0,φ).

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

H0(ρ,Φ)H02(ρ,Φ)2|A0(0,φ)|2ΔβΔSΔ2w0(Φ)|ρ|.
(36)

This “Type II” separable form agrees with the familiar FT of 1r. Figure 3 shows the profiles of two types of local frequency responses for n0 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 n0 is off image center.

Fig. 3
Type I and Type II local frequency responses H01(ρ, 0) and H02(ρ, 0) for n0 at image center in unweighted case: w(s, β) = 1. H01(0, 0) and H02(0, 0) are not shown because the value of H02(0, 0) blows up.

V. Quadratic Regularization: R0(ρ, Φ)

To evaluate the variance using (18) and (19), we also need the local frequency response of quadratic regularization, R0(ρ, Φ), [7], [8], [18], [19].

Practical regularization methods are based on the differences between neighboring pixel values. For a discrete-space 2D object μ[n], a typical quadratic roughness penalty is given in (4) and (5) for 1st-order differences. The rl[n] values are possibly space variant. For the purpose of local frequency response analysis, we examine the characteristics of R(μ) near a pixel n0 of interest, so we define rl,0 [triangle, equals] rl [n0] assuming rl[n] values vary smoothly. Then, the quadratic roughness penalty near a pixel n0 has the following form:

R(μ)=nl=1Lrl,012((clμ)[n])2.

The rl,0 values are design parameters that affect the directionality of the regularization and hence the shape of the PSF Each cl[n] is a (typically) high-pass filter. For a first-order difference:

cl[n]=ξl(δ2[n]δ2[nml]),

or for a 2nd-order difference:

cl[n]=ξl(δ2[n]δ2[nml])ξl(δ2[n]δ2[nml]),

where ξl = ‖mlυ/2, m = (nl, ml) denotes the spatial offsets to the neighboring pixels, and υ is the power of weights for diagonal neighbors that can be chosen by the user. For example, common practice chooses υ = 1 [20], [21].

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

R(μ)=l=1Lππππ12rl,0|Cl(ω)U(ω)|2dω(2π)2,
(37)

where μ[n]FTU(ω) and the DSFT of a Λ-order (where Λ [set membership] N) difference has the following magnitude:

|Cl(ω)|=ξl|1eι(ωml)|=ξl2sin(12(ωml)).

In the polar coordinates of (19):

|Cl(ρ,Φ)|2=|Cl(2πρΔeΦ)|2=ξl24sin2(πΔρeΦml).
(38)

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

R0(ρ,Φ)=R01(ρ,Φ)=l=1Lrl,0|Cl(ρ,Φ)|2=l=1Lrl,0ξl24sin2(πΔρeΦml).
(39)

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

|Cl(ρ,Φ)|2ξl2(mleΦ)2(2πΔρ)2=(2πΔρ)2ξl(12/υ)2cos2(Φφl),

where the angle between the lth neighbors is

φltan1mlnl.

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

R0(ρ,Φ)R02(ρ,Φ)=(2πρΔ)2R0(Φ),
(40)

where

R0(Φ)l=1Lξl(12/υ)2rl,0cos2(Φφl).

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

VI. Variance Prediction Implementation

Having obtained the approximations to H0(ρ, Φ), the local frequency response of the Gram operator given in (35) and (36), and to R0(ρ, Φ), 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 H01(ρ, Φ) given in (35) and R01(ρ, Φ) given in (39), and the Type II H02(ρ, Φ) given in (36) and R02(ρ, Φ) given in (40).

A. Double integral approach

The variance prediction using H01(ρ, Φ) in (35) and R01(ρ, Φ) 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 w0(Φ) can be precomputed, with the computation equivalent to one back-projection. We can coarsely sample Φ because P0(ρ, Φ) is fairly smooth along Φ.

B. Single integral approach

The separability properties of H02(ρ, Φ) in (36) and R02(ρ, Φ) 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:

Var{μ^[n]}02πζ/32|A0(0,Φ)|2w0(Φ)+α4π2ζR0(Φ)dΦ,
(41)

where ζ=ρmax3ΔβΔSΔ4 is a constant. We call this prediction the single integral (SI) approach. We evaluate this integral using a finite summation, with w0(Φ) and R0(Φ) precomputed. Therefore, computing (41) is equivalent to one back-projection.

C. Implementation of the single integral prediction

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 [var phi]1,2,3,4 = 0, π/2, π/4, and −π/4 and rl,0 = (1,1,1,1). We used υ = 1 both in calibration and reconstruction as is the common practice in quadratic regularization. For these choices, the Type II R02 (ρ, Φ) in (40) becomes independent of Φ:

R0(Φ)=l=1Lml=1+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.

VII. Simulation Results

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

d0(φ)d(φ)2z1z2z12sin2φ+z22cos2φ,
(43)

where z1 = 244.1 mm and z2 = 220.7 mm are the semi-major and semi-minor axes of the ellipse. This approximation is exact when n0 is at the ellipse center. The approximation to d([var phi]) becomes less accurate as xc[n0] approaches the edge of the ellipse support.

For simplicity, we model the detector response2 in (23) by a shift-invariant rectangle of width Δs:

b(s)=1ΔSrect(sΔS).

In the case of a square pixel basis χ(x) = rect2(x)3, we have from (51) (see Appendix I)

A0(v,φ)=sinc(ΔSvm0(φ))Δ2sinc(vΔcosφ)sinc(vΔsinφ),
(44)

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

m0(φ)mc(φ)=mc=Dsd/Ds0,

where mc is the value of m0([var phi]) at the image center.

We chose the regularization parameter α = 211 to give FWHM = 1.72 pixels, i.e., 3.4 mm, at the center of the image. For each realization, [mu] 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 H01 (ρ, Φ) in (35) and R01 (ρ, Φ) in (39), and the SI approximation (41) with R0(Φ) 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.

A. Standard Quadratic Penalty Function

We first considered a standard quadratic penalty where

rl,0=κc2,

and κc2 is the value of κ2 [n0] at the image center in (45) below. This choice assures that the resolution of the two studies is matched at the image center. R0(Φ)=(1+2)κc2 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.

Fig. 4
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: R0=(1+2)κc2.

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.

B. Certainty-based Quadratic Penalty Function

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:

rl,0=κ2[n0],

where

κ2[n0]12π02πw(s(r0(φ)),β(r0(φ),φ))dφ.
(45)

Here, R0(Φ) 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.

Fig. 5
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 d0([var phi]).

C. Computation Time and Accuracy

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 w0(Φ) and G0(ρ, Φ). The precomputation time for w0(Φ) was about 19 seconds on dual Intel Xeon 3.40GHz CPU. The precomputation time for G0(ρ, Φ) with 492 samples in [var phi] 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 κc2 case, the NRMS percent errors for FFT, DI and SI are 6.6%, 6.8% and 6.6%; for κ02 case, the NRMS percent errors for FFT, DI and SI are 6.5%, 6.0% and 8.3%, respectively.

VIII. Conclusion and Discussion

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(ss′) is replaced by b(s, s′) in (24) and b0(r, [var phi]) in (52) becomes

b0(r,φ)m0(φ)b(s0(φ)+m0(φ)r,s0(φ)),

where

s0(φ)Dsdarcsinr0(φ)Ds0.

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.

Acknowledgments

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

Appendix I: Derivation of Local Impulse Response

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

ss(r)=Dsdarcsin(r/Ds0)ββ(r,φ)=φarcsin(r/Ds0).

Then the fan-to-parallel beam rebinning of a(s, β; n) is

a(s(r),β(r,φ);n)b(s(r)s(r))Δg(rrφ[n]Δ,φ)|s˙(r)|dr,a(r,φ;n)
(46)

because r(s(r′)) = r′ and [var phi](s(r′), β(r, [var phi])) ≈ [var phi] for rr′.

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

s(r)s(r)s˙(r)(rr).

Substituting into (46), the system matrix elements become

a(r,φ;n)b(s˙(r)(rr))Δg(rrφ[n]Δ,φ)|s˙(r)|dr.
(47)

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

hd[n;n]1Δβ1ΔS02πw¯(r,φ)a(r,φ;n)a(r,φ;n)|J(r)|drdφ=1Δβ1ΔS02πw(φ;n;n)hφ[n;n]dφ,
(48)

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

hφ[n;n]a(rrφ[n],φ;n)a(rrφ[n],φ;n)dr
(49)
w(φ;n;n)w¯(r,φ)|J(r)|a(rrφ[n],φ)a(rrφ[n],φ)dra(rrφ[n],φ)a(rrφ[n],φ)drw¯(r,φ)w(s(r),β(r,φ)).

Let r0([var phi]) [triangle, equals] r[var phi][n0]. Because s(r) is fairly smooth, we make the following approximation for r′ ≈ r0([var phi]):

s˙(r)s˙(r0(φ))m0(φ).
(50)

Substituting (50) into (47) and simplifying yields

a(r,φ;n)a(r,φ;n0)a0(r,φ)=b0(rrφ[n]r,φ)Δg(rΔ,φ)dr,
(51)

with r′ = r′ − r[var phi][n] and

b0(r,φ)m0(φ)b(m0(φ)r).
(52)

Therefore, we further simply (48) as follows

hd[n;n]1Δβ1ΔS02πw0(φ)h0(Δ(nn)eφ,φ)dφ,
(53)

where because w(r, [var phi]) often varies slowly in r relative to the typically sharp peak of a0(r, [var phi]) at r = 0,

h0(r,φ)a0(r,φ)a0(r,φ),w(φ;n;n)w¯(r,φ)|m0(φ)|a02(rrφ[n0],φ)dra02(rrφ[n0],φ)dr|m0(φ)|w¯(r0(φ),φ)w0(φ),
(54)

where [large star] denotes a 1D autocorrelation with respect to r.

Appendix II: Derivation of Local Frequency Response

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:

hφ[n;n0]h0(Δneφ,φ)η1(Δn)

where

η1(x)tri(xeφd0(φ)),
(55)

and d0([var phi]) denotes the length of the chord through n0 through the FOV at angle ([var phi] + π/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:

hφ[n;n0]hφ[n0+n/2;n0n/2]=h0(Δneφ,φ)η2(Δn),

where x0 [triangle, equals] xc[n0], and η2 denotes the support of the image,

η2(x)η(x0+x/2)η(x0x/2).
(56)

This approach yields a local impulse response that is symmetric in n, thus ensuring that its spectrum is real.

Another alternative is to refer all displacements relative to the point n0 as follows:

hφ[n;n0]hφ[n0+n;n0]=h0(Δneφ,φ)η3(Δn),

where

η3(x)η(x0+x)η(x0).
(57)

This choice is not symmetric in n but it better corresponds to the local Fourier analysis based on the DFT of A′WAej. For simplicity, we could also approximate (57) as follows:

η4(x)η0(x)=η(x)η(x0).
(58)

This choice also yields a local impulse response that is symmetric in n provided η(x) is symmetric itself.

We focus on η1 (x) in (55) hereafter because it preserves the property of nonnegative definiteness [33]. Define the following “strip like” function:

sφ(x)h0(xeφ,φ)η1(x),

and sφ(x)2D FTSφ(u,v). Taking the DSFT of (53) yields the following result:

Hd0(ω)=1Δβ1ΔS02πw0(φ)Hφ(ω)dφ,
(59)

where H[var phi]([omega with right arrow above]) is the spectrum of h[var phi][n;n0] = s[var phi]n), as follows:

Hφ(ω)=nsφ(Δn)eι(ωn)1Δ2sφ(x)eι1Δ(ωx)dx=1Δ2Sφ(ω2πΔ).
(60)

The 2D FT of s[var phi](x) is easiest to see for the case [var phi] = 0:

s0(x,y)=h0(x,0)tri(yd0(0))2D FTS0(u,υ)=|A0(u,0)|2d0(0)sinc2(d0(0)υ),

where A0(ν, [var phi]) is associated with the detector response and basis effect, given in (44). By the rotation property of the 2D FT:

Sφ(ρ,Φ)|A0(ρcos(Φφ),φ)|2d0(φ)sinc2(d0(φ)ρsin(Φφ)).

Therefore, using (54) and (60), the local frequency response H0(ρ,Φ) around a point n0 is

H0(ρ,Φ)1ΔβΔsΔ202πw0(φ)Sφ(ρ,Φ)dφ.
(61)

Footnotes

1Throughout the paper we use the subscript “0” to indicate dependence on a given pixel location n0.

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

3rect2(x) [triangle, equals] rect(x)rect(y) is a 2D square function.

References

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.