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 July 28.
Published in final edited form as:
PMCID: PMC2911484
NIHMSID: NIHMS217495

Quadratic Regularization Design for 2D CT

Hugo R. Shi, Student Member, IEEEcorresponding author and Jeffrey A. Fessler, Fellow, IEEE

Abstract

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

Index Terms: Fan-Beam CT, Iterative Reconstruction, Regularization, Local Impulse Response, Spatial Resolution

I. Introduction

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

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

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

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

II. Local Impulse Reponse

Let y = (y1, …, yM) denote the vector of noisy transmission sinogram measurements recorded in a CT system. For simplicity, we consider the following mono-energetic formulation for the ensemble mean of the measured data:

y¯i(x)=E[yi]=bie[Ax]i+ri,
(1)

where A is the system matrix, x = (x1, x2, …, xN) is a discretized version of the object being imaged, bi denotes the blank scan, ri denotes the additive contribution of scatter, and [Ax]i=j=1Naijxj. The spatial resolution and therefore regularization design techniques are applicable to PL methods when appropriate weightings are chosen [3]. We estimate x by minimizing the following cost function:

Ψ(y,x)(y)AxW2+ζR(x)
(2)
x^(y)=argminxΦ(y,x),
(3)

where (y)ln(yrb) using element-wise division, R(x) is a regularizer that controls noise, ζ is a scalar that determines the resolution-noise tradeoff, and W = diag{wi} is a statistical weighting matrix. For transmission tomography with the model in (1), usually we have the plug-in weighting wi = yi corresponding to the PL estimator for (1) [12]. For emission tomography, the plug-in weighting is wi = 1/yi; for low count levels alternative weightings are preferable [13]. The regularization design methods in this paper apply to other statistical models than (1) by simply changing the weighting matrix W.

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

ljlimε0x^(y¯(xtrue+εδj))x^(y¯(xtrue))ε=yx^(y)|y=y¯(xtrue)xy¯(x)|x=xtrueδj,
(4)

where δj is a Kronecker impulse at the jth pixel and the gradient operations are matrices such that [yx^(y)]ji=yix^j(y) and [xy¯(x)]ij=xjy¯i(x) (this definition depends on the mean of the data, rather than the mean of x̂ and differs slightly from that used in [3], but yields the same final expression without requiring further approximations). Evaluating (4) [3], the LIR for the PWLS estimator (3) is

lj=lj(xtrue,R)=[AW A+ζR]1AW Aδj,
(5)

where R is the Hessian of the regularizer R(x).

As is evident from (5), the local impulse response depends on the regularizer through its Hessian, R. We would like to design R such that the local impulse response lj best matches a target response l0 at every pixel j. We could phrase this matrix optimization problem as:

argminRjSjlj(xtrue,R)l02,
(6)

where Sj is a shift operator that recenters lj to the center pixel. This matrix formulation of the design problem seems intractable, so we turn to the frequency domain to simplify the problem.

III. Frequency Domain Analysis

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

A. Discrete Fourier Analysis

An analytical expression for A′W Aδj is derived in Appendix §B, showing that A′W Aδj changes slowly with j. Thus A′W A is approximately locally circulant. The R matrix in (5) can be designed to be locally circulant as well. This allows us to factor out a Fourier transform matrix in (5).

Let Q denote an orthonormal discrete Fourier transform matrix centered at pixel j (we recenter the Fourier matrix to avoid the complex exponentials caused by non-centered impulse responses). Then A′W AQ′ Λj Q, where Λj=diag{λkj} and λj [triangle, equals] F(A′W Aδj), and RQ′ΓjQ, where Γj=diag{γkj} and γj [triangle, equals] F(j), where F(·) denotes the 2D DFT [8], [14]. Then we can approximate the LIR in (5) as:

lj[QΛjQ+ζQΓjQ]1QΛjQδj
(7)
=Q[ΛjΛj+ζΓj]Qδj,
(8)

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

LjΛjΛj+ζΓjL0.
(9)

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

Γ^j=argminΓTL0ΛjΛj+ζΓj,
(10)

where L0 is the frequency response of l0 and An external file that holds a picture, illustration, etc.
Object name is nihms217495ig1.jpg denotes the set of possible frequency responses for the regularizer limited by its structure that will be enumerated in §III-B. Alternatively, as a preview of the methods in §IV-A, we can cross multiply the terms in (9) yielding the simpler design criterion:

Γ^j=argminΓTL0(Λj+ζΓj)Λj.
(11)

This approach is similar to the formulation in [6]. Because Λj is the DFT of A′W Aδj, calculating Λj requires one forward projection and backward projection per pixel. Thus, regularization design based on (10) or (11) would be very slow.

Prior to [10], regularization design methods were based on discrete Fourier transforms. As shown in Appendix B, the continuous-space analog of Λj in polar coordinates (ρ, Φ) is the frequency response wj(Φ)|ρ|, where wj(Φ) is an expression that incorporates the Jacobian from the change of coordinates from parallel-beam to fan-beam geometry, and weights from W that correspond to rays that intersect pixel j at angle Φ. Substituting this into (9) and using continuous space analogs of Γj yields the following expression for the continuous space analog of Lj:

j(ρ,Φ)wj(Φ)1|ρ|wj(Φ)1|ρ|+ζRj(ρ,Φ)=wj(Φ)wj(Φ)+ζ|ρ|Rj(ρ,Φ),
(12)

where Rj(ρ, Φ) is the local frequency response of the regularizer near pixel j. The continuous space analogs of Λj and Γj simplify (12) to provide a more efficient approach to regularization design than (11), as detailed below in (27).

B. Regularization Structure

Regularizers control roughness by penalizing differences between neighboring pixels. Indexing the image x as a 2D function x(n, m), we define a first-order differencing function for a regularizer that penalizes the lth direction as

cl(n,m)=1nl2+ml2(δ(n,m)δ(nnl,mml)),
(13)

where typically (nl, ml) [set membership] {(0, 1), (1, 0), (1, 1), (−1, 1)}, corresponding to horizontal, vertical, and diagonal differences. A conventional quadratic regularizer can then be expressed as

R(x)=n,ml=1L12|(clx)(n,m)|2,
(14)

where ** denotes 2D convolution. This conventional regularizer assigns the same weight to the differences between each neighbor. In this paper, we make the regularizer spatially adaptive with the addition of weighting coefficients rlj as follows:

R(x)=n,ml=1Lrlj12|(clx)(n,m)|2,
(15)

where j is the columnized pixel index which is a function of (n, m). The objective of this paper is to design coefficients {rlj}. To this end, we must analyze the local frequency response Rj(ρ, Φ) of the Hessian R of the space-variant regularizer (15).

Taking the Fourier transform of (13) yields:

|Cl(ω1,ω2)|2=1nl2+ml2|1ei(ω1nl+ω2ml)|2=1nl2+ml2(22cos(ω1,nl+ω2,ml).
(16)

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

C. Efficient Approximation for R

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

|Cl(ω1,ω2)|21nl2+ml2(ω1nl+ω2ml)2.
(17)

One can think about isotropy intuitively in polar coordinates as eliminating angular dependence. Therefore we convert (17) to polar frequency coordinates to simplify the analysis. Using frequency and sampling relationships from the DFT, ω1 = 2πΔxρ cos(Φ) and ω2 = 2πΔyρ sin(Φ), where Δx, Δy denotes pixel size and (ρ, Φ) are polar frequency coordinates. For simplicity, we assume Δx = Δy = 1 for square pixels. Then:

|Cl(ω1,ω2)|21nl2+ml2(nl2πΔxρcos(Φ)+ml2πΔyρsin(Φ))2=1nl2+ml2(2πρ)2(nlcos(Φ)+mlsin(Φ))2=(2πρ)2cos2(ΦΦl),
(18)

where Φltan1mlnl. Taking the Fourier transform of (15) and combining with (18), our final expression for the local frequency response of the regularizer near the jth pixel is

Rj(ρ,Φ)=l=1Lrlj(2πρ)2cos2(ΦΦl).
(19)

For the usual choice of L = 4 and for (nl, ml) described below (13), Φl [set membership] {0, π/2, π/4, 3π/4}.

D. Target Local Frequency Response

Substituting (19) into (12) yields a simple expression for the local frequency response of a PWLS estimator. We want to design each regularization coefficient vector rj=(r1j,,rLj) such that our frequency response matches that of the target as closely as possible. We know that the local frequency response associated with PULS is isotropic at the center of the field of view so we select that to be our target.

Using (nl, ml) [set membership] {(0, 1), (1, 0)}, the squared magnitude response of a conventional regularizer in (14) using the AIMA approximation is,

R0(ρ,Φ)=(2πρ)2
(20)

and without the approximation is

R0(ρ,Φ)=42cos(2πρcosΦ)2sin(2πρsinΦ).
(21)

For an unweighted cost function and a parallel beam geometry, the continuous-space frequency response that is analogous to F(A′ Aδj) is 1|ρ|. As shown in (49) in Appendix §B, for uniform weights (wi = 1) we have the following local frequency response for fan-beam geometries

Hj(ρ,Φ)=2J(sj(Φ))|ρ|,
(22)

where J(sj(Φ)) is the Jacobian for the change of coordinates from parallel-beam to fan-beam geometries as defined in (44), and sj(Φ) is an index into the sinogram based on pixel j and angle Φ. Setting sj(Φ) = 0 to correspond with the center pixel in (49), the target local frequency response is:

02J(0)|ρ|2J(0)|ρ|+ζR0(ρ,Φ)
(23)
=11+ζ|ρ|0.5J(0)R0(ρ,Φ),
(24)

where An external file that holds a picture, illustration, etc.
Object name is nihms217495ig2.jpg0 is the continuous-space analog of L0 in (9). Next we will use these continuous-space analogs to solve for the regularization coefficients rj.

IV. The AIMA Approach

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

A. Solving for regularization coefficients

We try to design regularization coefficients rj to match our designed local frequency response (12) to the target local frequency response (24) as follows:

wj(Φ)wj(Φ)+ζ|ρ|Rj(ρ,Φ)11+ζ|ρ|0.5J(0)(2πρ)2.
(25)

Cross multiplying and simplifying yields

wj(Φ)ζ|ρ|0.5J(0)(2πρ)2ζ|ρ|Rj(ρ,Φ)wj(Φ)0.5J(0)(2πρ)2l=1Lrlj(2πρ)2cos2(ΦΦl)wj(Φ)wj(Φ)0.5J(0)l=1Lrljcos2(ΦΦl).
(26)

We design rj by minimizing the difference between both sides of (26):

rj=argminr012π02π(wj(Φ)l=1Lrlcos2(ΦΦl))2dΦ.
(27)

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

We can think of the minimization in (27) as a projection of wj = (Φ) onto the space spanned by {cos2(− Φl)}, which allows us to greatly simplify the problem and derive a computationally efficient and mostly analytical solution to the regularization design problem. Expanding these cosines in terms of a three term basis that is orthonormal with respect to the inner product f,g=12π02πf(Φ)g(Φ)dΦ yields

cos2(ΦΦl)=121+cos(2Φl)22[2cos(2Φ)]+sin(2Φl)22[2sin(2Φ)].
(28)

The three orthonormal basis functions are

p1(Φ)=1p2(Φ)=2cos(2Φ)p3(Φ)=2sin(2Φ).

Using (27), we write l=1Lrljcos2(ΦΦl)=PTrj, where P is an operator whose columns are p1, p2, and p3, and T is a 3 × L matrix of linear combination coefficients whose lth column is [12cos(2Φl)22sin(2Φl)22]T.. T is computed by taking the inner products of cos2(Φ − Φl) and p1, p2, and p3, i.e., Tmn = ∫ cos2(Φ − Φm) pn(Φ)dΦ.

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

rj=argminr0Trbj2,
(29)

where bj [triangle, equals] P*wj (·) and P* denotes the adjoint of P, i.e., bkj=pk(Φ)wj(Φ)dΦ, k = 1, 2, 3.

B. Zeros in the Hessian

If there are too many zeros in rj, there will be zeros in the Hessian, possibly degrading x̂. This can cause elongated impulse responses that may contribute to streak artifacts in the reconstructed image. This phenomenon occurred in [10] however we did not notice the artifacts due to the coarser spatial resolution in PET. It is present when using AIMA with the ring phantom in §V. The problem improves when using FIIN method of §A, thus we believe this phenomenon is caused by the approximation error in (17) as well as the non-negativity constraint. For AIMA, and to be safe using FIIN, we modify (29) to ensure that an adequate number of rlj values are greater than some εlj>0. Requiring the penalty coefficients for the vertical and horizontal directions be nonzero is sufficient to eliminate zeros in the Hessian. (A similar constraint could be created using the 2 diagonals instead of the vertical and horizontal neighbors; however the approximation in (17) is worse for diagonal neighbors). We turn to previous work to select εlj. In [3], we derived a certainty based weighting using a spatially variant (κj)2=iaij2wiiaij2 that seeks to provide uniform spatial resolution. In terms of the continuous space analogs used in this paper, (κj)2=1π0πwj(Φ)dΦ. This regularization design method can be implemented using the regularization structure presented in this paper by setting rlj=(κj)2 for (ml, nl) = (0, 1), (1, 0). This approach provides a convenient nominal value for the regularization strength at each pixel. We define the lower constraint vector εj such that εlj=α(κj)2 for vertical and horizontal neighbors, (ml, nl) = (0, 1), (1, 0), and εlj=0 for all other neighbors. Using a nonzero coefficient α improves noise at the expense of isotropy. For the results presented in this paper, we used α = 0.1.

Now we formulate our problem so that non-negative least squares (NNLS) algorithms will accommodate this new constraint. Let [r with macron]j [triangle, equals] rj + εj. Solving with the constraint of rj ≥ 0 ensures that [r with macron]jεj. Substituting [r with macron]j into (29) yields T[r with macron]jbj = T(rj + εj) − bj = Trj − (bjj). So our final cost function for regularization design is

r^j=argminr0Tr(bjTεj)=argminr0Trb¯j,
(30)

where bj [triangle, equals] bjj. We use coefficients [r with macron]j = [r with circumflex]j + εj in the regularizer (15). We next solve (30) analytically. It is this analytical solution that makes the AIMA more efficient than the FIIN method of Appendix §A.

C. A Mostly Analytical Solution

Using a second-order neighborhood (L = 4) we select (nl, ml) to be (1, 0), (0, 1), (1, 1), (1, −1) leading to the following Φl: Φ1 = 0, Φ2 = π/2, Φ3 = π/4, Φ4 = −π/4. So the terms in (30) are

T=12[11111/21/200001/21/2]b¯j=[d1j2d2j2d3j],dj=[12π(1α)02πwj(Φ)dΦ12π02πwj(Φ)cos(2Φ)dΦ12π02πwj(Φ)sin(2Φ)dΦ].
(31)

Observe that j = [α(κj)2 0 0]T, so εj affects only d1. Ignoring α, d1 is the continuous space analog of the certainty-based regularization weighting in [3], d2 is related to the horizontal and vertical directions, and d3 is related to the diagonal directions.

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

We could solve the minimization (30) using an iterative NNLS algorithm [15, p. 158]. However, using the properties of T and d, we can avoid iterations almost entirely by solving (30) analytically using the KKT conditions. When α = 0, d22+d32d1. as outlined in Appendix C. This inequality is usually true for the values of α used (typically around 0.1), however for pixels where it is not, the problem would have to be solved using a NNLS algorithm rather than this analytical solution. To use NNLS, one must add Tikhonov regularization to help (30) converge to a minimum norm solution which ensures that rj is a continuous function of dj which is a continuous function of wj(Φ). For more details on using NNLS to solve the regularization problem, see §IV-A in [9]. For the results presented in this paper, all pixels of both simulations satisfie d22+d32d1 for all pixels and NNLS was not used. The structure of T leads to eight-fold symmetry that simplifies analysis. If d2 < 0 we can solve for r using |d2| and then swap r1 with r2. If d3 < 0 we can solve for r using |d3| and then swap r3 with r4. If d3 > d2 we can solve for r with d2 and d3 interchanged, and then swap r1 with r3 and r2 with r4. Therefore, hereafter we focus on cases where 0 ≤ d3d2d1. Fig. 1 shows these first octant cases, numbered according to the number of nonzero elements of r.

Fig. 1
First octant of quadratic penalty design space showing the four regions where different constraints are active.
  1. If d212d1 and d323d213d1, then r1=43(d1+d2), r2 = r3 = r4 = 0.
  2. If d323d213d1 and d3+d212d1, then r1=85[12d1+32d2d3], r2 = r4 = 0, r3=125[d3(23d213d1)].
  3. If d3+d212d1 and d214d1, then there are multiple non-negative r that exactly solve [nabla]Ψ(r) = 0. The minimum-norm solution is r1 = 4d2, r2 = 0,
    r3 = d1 − 2d2 + 2d3, r4=2[12d1(d2+d3)].
  4. If d214d1, then there are multiple non-negative r that are exact solutions. The natural choice is the minimum-norm r given by the pseudo-inverse solution r = Td, where r1=2(14d1+d2), r2=2(14d1d2), r3=2(14d1+d3), r4=2(14d1d3).

The analytical solution presented above is for the usual first-order differences (13). For higher-order differences or neighborhoods, it would appear to become increasingly cumbersome to solve (27) analytically, so an iterative NNLS approach may be more appealing. This can still be practical since T is quite small. The analytical solution above is a continuous function of d, which in turn is a continuous function of wj(Φ). This continuity property would seem to be desirable for avoiding artifacts in the reconstructed images.

For practical implementation, we simply discretize the integrals in (31) [16]. This presents interpolation issues in extracting a discretized version of wj(Φ), from W. For the fan beam case, we use the analytical formula for wj(Φ) presented in Appendix §B (48). Those equations are presented with continuous space coordinates s that spans the length of the detector, and β, the angle of rotation of the system. In discrete implementation, we round these off to the nearest neighbor to get indices into W. The analytical solution presented in this section is very efficient over all values of j, and the bulk of the computation time is spent computing wj(Φ) which has the compute time of approximately one back-projection. For further speed optimizations, down-sampled back-projections could be explored. In practice, wj(Φ) is fairly smooth and the basis functions that we use to approximate wj(Φ) are even smoother. One can achieve good results by calculating wj(Φ) for fewer angles and interpolating the rest. The performance and efficacy of such a scheme still needs to be explored.

V. Results

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

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

We simulated a noiseless fan-beam sinogram without scatter using a phantom consisting of a background disk and two rings each of thickness 1mm shown in Fig. 2. We generated the noiseless sinogram by taking line-integrals through the analytical phantom with the same system geometry. We used the plug-in weighting wi = yi for this experiment. Fig. 3 shows penalty coefficients {rj} designed using the method of section IV-C that yielded the following results. The top images show coefficients in the horizontal and vertical neighbors, and the bottom images show coefficients for diagonal neighbors. There are substantial spatial variations in these coefficients.

Fig. 2
Ring phantom used for simulation study
Fig. 3
Regularization penalty coefficients used in reconstruction of ring phantom. The four images are rlj for l = 1, …, 4

We reconstructed images x̂ using several methods. We first created an image uniformly blurred by the target PSF. We then reconstructed images using (i) conventional regularization, rj = 1, (ml, nl) = (0, 1), (1, 0), (ii) certainty based regularization [3], (iii) AIMA method with α = 0.1, (iv) AIMA method with α = 0, and (v) FIIN method with α = 0. Fig. 4 is a closeup of the right-most ring reconstructed with the various methods listed above and Fig. 5 is a closeup of the left-most ring.

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

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

Fig. 6
Profiles through the right-most ring from each reconstructed image.
Fig. 7
Profiles through the left-most ring from each reconstructed image.

The second study used a similar imaging geometry using one slice of real CT data from a GE scanner described in [18] with weightings W computed in the same way as [19]. ζ was selected such that the target PSF had a FWHM of 1.51 mm.

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

Fig. 8
Reconstruction with conventional regularization without windowing with impulse responses marked
Fig. 9
Reconstructions windowed from 800 to 1200 HUs using, from top to bottom, conventional regularization, certainty based regularization, AIMA with α = 0.1, AIMA with α = 0, FIIN with α = 0.

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

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

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

Fig. 14
Plots of the error histogram for different impulse responses.

VI. Discussion

A. Spatial Resolution Properties

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

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

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

The resulting impulse responses from the AIMA and FIIN methods are not completely isotropic. This may seem to contradict the dramatic improvement these regularization design methods achieved with the ring phantom. However, recall that we are trying to approximate wj(Φ) using 3 basis functions (see §IV-A) for the AIMA, and 4 basis functions (see §A) for FIIN. With real data, wj(Φ) is a complicated function that cannot be parameterized using 3 or 4 basis functions. This aspect, along with the non-negativity constraint limits the performance of any regularization design technique with a finite number of parameters. wj(Φ) is much simpler for simple phantoms like the ring phantom, so AIMA yields better results there. Extensions of this regularization design to higher order penalties have the potential for more basis functions, and better performance.

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

B. Computation Time

AIMA is quite efficient. Computing certainty based regularization takes the time of about 1 backprojection. In AIMA, we must first compute wj(Φ) which takes the time of about 1 backprojection, and then solve the analytical solution which is very fast. FIIN also requires 1 backprojection to compute wj(Φ), however it then has to run a non-negative least squares problem for every pixel. Though this is a fairly small NNLS problem, T is 4 × 4, it adds much compute time since it must be calculated for each pixel. In general, due to the faster compute times, we recommend AIMA with α = 0.1. If accuracy is more important than compute time, FIIN can be used instead.

VII. Conclusion and Future Work

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

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

Acknowledgments

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

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

Appendix A: Slower Regularization Design

This appendix describes the slower Double Integral Iterative NNLS, FIIN, method that does not use the AIMA approximation 2 − 2 cos(x) ≈ x2 in (17) in the expression of the local frequency response of R. This design method is appropriately named because it uses the full integral over all frequency domain polar coordinates instead of just the angular coordinates. The matrices involved in this method are slightly more complex than those in AIMA so an iterative technique is needed. Converting (16) in §III-C into polar coordinates and combining with (15) yields

Rj(ρ,Φ)=l=1Lrljnl2+ml2(22cos(ul(ρ,Φ))),

where ul(ρ, Φ) = nl2πΔxρ cos(Φ) +ml2πΔyρ sin(Φ). The target we use in regularization is the LIR associated with PULS at the center of the image, expressed in (23).

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

wj(Φ)wj(Φ)+β|ρ|Rj(ρ,Φ)11+β|ρ|0.5J(0)R0(ρ,Φ)wj(Φ)+β|ρ|Rj(ρ,Φ)wj(Φ)(1+β|ρ|0.5J(0)R0(ρ,Φ))β|ρ|Rj(ρ,Φ)β|ρ|0.5J(0)wj(Φ)R0(ρ,Φ))Rj(ρ,Φ)wj(Φ)R0(ρ,Φ)),

where wj(Φ) = 0.5J(0)wj(Φ). To design R, we use a method similar to that of §IV-A and project wj(Φ)R0(ρ, Φ)) onto the space spanned by {1 − cos(ul(ρ, Φ))} that can be orthonormalized into L basis functions, {pl(ρ, Φ)} using Gramm-Schmidt. Then our regularization design problem simplifies to

rj=argminr0Trbj2,
(32)

where bj is a vector of inner products between wj(Φ)R0(ρ, Φ)) and the L orthonormal basis functions {pl(ρ, Φ)}, bkj=pk(ρ,Φ)wj(Φ)R0(ρ,Φ)dρdΦ, and T is a L × L matrix whose elements are the inner products between 2 − 2 cos(ul(ρ, Φ)) and {pl(ρ, Φ)}, Tmn = ∫ ∫ 2 − 2cos(um(ρ, Φ))pn(ρ, Φ)dρdΦ. We define P to be an operator whose columns are {p1, p2, …pL} so R(ρ, Φ) = PTrj. This design problem is then solved with an NNLS algorithm. This more accurate version is slower than that presented in §IV-A because there is no apparent analytical solution similar to the one presented in §IV-C.

Appendix B: Analysis of Grammian Operator

This appendix considers fan-beam geometries and uses continuous-space analysis to analyze the Fourier transform of the Grammian operator A′W A to simplify the regularization design problem in (12). One can use polar coordinates Φ, ρ, and continuous-space analysis to separate the angular and radial components of A′W A. Using this framework, isotropy can be thought of as eliminating dependence on the angular component Φ.

A. Fan-Beam Geometry

Fig. 15 illustrates the fan-beam geometry that we consider. Let P be the rotation isocenter. D0d denotes the distance from the point P to the detector, Ds0 denotes the distance from the X-ray source to P, and Dfs denotes the distance from the X-ray source to the focal point of the detector arc. Define Dsd [triangle, equals] D0d + Ds0 to be the total distance from the X-ray source to the center of the detector, and Dfd [triangle, equals] Dsd + Dfs to be the total distance from the focal point to the center of the detector. This formulation encompasses a variety of system configurations by allowing the detector focal point to differ from the X-ray source location. For flat detectors, Dfs = ∞. For third-generation X-ray CT systems, Dfs = 0. For fourth generation X-ray CT systems, Dfs = −Ds0.

Fig. 15
Illustration of fan beam geometry.

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

α(s)=sDfdγ(s)=tan1((Dfd)sinα(s)(Dfd)cosα(s)Dfs)

where the two most important cases are

γ(s)={α(s),Dfs=0tan1s/Dsd,Dfs=.
(33)

The relationship between s and γ is:

s={(Dfd)[γarcsin(DfsDfdsinγ)],Dsdtanγ,Dfs=.0Dfs<
(34)

The ray corresponding to detector element s and angle β is An external file that holds a picture, illustration, etc.
Object name is nihms217495ig2.jpg(s, β) = {(x, y) : x cos [var phi](s, β) + y sin [var phi](s, β) = r(s)}, where

ϕ(s,β)β+γ(s)r(s)Ds0sinγ(s).
(35)

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

|r(s)|rmaxDs0sinγmax,
(36)

where γmax [triangle, equals] γ(smax) and smax is half the total detector arc length. The radius rmax defines the circular field of view of the imaging system. The fan angle is 2 γmax.

The line-integral projection p(s, β) of f along An external file that holds a picture, illustration, etc.
Object name is nihms217495ig2.jpg(s, β) is2:

p(s,β)=(s,β)f(x,y)d=f(x,y)δ(xcos(ϕ(s,β))+ysinϕ(s,β)r(s))dxdy,
(37)

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

The usual inner product for fan-beam projection space is

p1,p2=smaxsmax0βmaxp1(s,β)p2(s,β)dsdβ.

Using this inner product in projection space, and the usual An external file that holds a picture, illustration, etc.
Object name is nihms217495ig2.jpg2 inner product in image space, the adjoint of An external file that holds a picture, illustration, etc.
Object name is nihms217495ig3.jpg is given by

(Pp)(x,y)=smaxsmax0βmaxp(s,β)δ(xcosϕ(s,β)+ysinϕ(s,β)r(s))p(s,β)dsdβ,
(38)

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

B. Parallel-Beam Grammian Analysis

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

PWPδj(x,y)=0πδ(xjcos(φ)+yjsin(φ)r)w(r,φ)δ(xcos(φ)+ysin(φ)r)drdφ.

Using the sampling property with the first δ, define

wj(φ)w(r,φ)=w(xjcos(φ)+yjsin(φ),φ).
(39)

We denote the Fourier transform of δj(x, y) as Fj(ρ, Φ). Then, using the Fourier slice theorem,

PWPδj(x,y)=0πFj(ρ,φ)ei2πρrdρwj(φ)δ(xcosφ+ysinφr)drdφ=0πFj(ρ,Φ)wj(φ)ei2πρ(xcosφ+ysinφ)dρdφ.
(40)

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

PWPδj(x,y)=0πwj(φ)F(δj(x,y))ρei2πρ(xcosφ+ysinφ)dρdφ.

The local impulse response is recentered, hj(x, y) = An external file that holds a picture, illustration, etc.
Object name is nihms217495ig3.jpg*An external file that holds a picture, illustration, etc.
Object name is nihms217495ig4.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms217495ig3.jpgδj(x+xj, y+yj) which eliminates the complex phase term Fj(ρ, [var phi]). Then, the local frequency response is

Hj(ρ,Φ)=wj(Φ)ρ.
(41)

C. Fan Beam Grammian Analysis

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

(PWPδj)(x,y)=βmaxβmaxsmaxsmaxw(s,β)δ(xjcos(ϕ(s,β))+yjsin(ϕ(s,β))r(s))δ(xcos(ϕ(s,β))+ysin(ϕ(s,β))r(s))dsdβ.

We use the change of variables

r=r(s)=Ds0sinγ(s)
(42)
ϕ=ϕ(s,β)=β+γ(s)
(43)

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

J(s)=|Ds0cosγ(s)||γ˙(s)|.
(44)

Then,

(PWPδj)(x,y)=02πrmaxrmaxw(r,ϕ)δ(xjcosϕ+yjsinϕr)δ(xcosϕ+ysinϕr)1J(s(r))drdϕ.

In this expression,

w(r,ϕ)w(s(r),β(r,ϕ))s(r)=γ1(arcsin(r/Ds0))β(r,ϕ)=ϕarcsin(r/Ds0).

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

wj(ϕ)w(xjcosϕ+yjsinϕ,ϕ)
(45)
sj(ϕ)s(r)r=xjcosϕ+yjsinϕ.
(46)

Again, let Fj(ρ, Φ) denote the Fourier transform of δj(x, y). Then, using the Fourier slice theorem,

(PWPδj)(x,y)=rmaxrmax02πFj(ρ,Φ)wj(ϕ)J(sj(ϕ))ei2πprδ(xcosϕ+ysinϕr)dρdϕdr=02πFj(ρ,Φ)wj(ϕ)J(sj(ϕ))ei2πρ(xcosϕ+ysinϕ)dρdϕ=0πFj(ρ,Φ)wj(ϕ)J(sj(ϕ))ei2πρ(xcosϕ+ysinϕ)dρdϕ+0πFj(ρ,Φ)wj(ϕ)J(sj(ϕ))ei2πρ(xcosϕ+ysinϕ)dρdϕ,

where [var phi]″ = [var phi] + π. This is similar to the parallel-beam derivation except that we have two integrals. We can convert each integral into the inverse Fourier transform as we did in the parallel-beam case and strip out the phase term Fj(ρ, Φ) by recentering the local impulse response. The local frequency responses of the Grammian operator is

Hj(ρ,Φ)=wj(Φ)|ρ|,
(47)

where,

wj(Φ)=1J(sj(Φ))[wj(ϕ)|ϕ=Φ+wj(ϕ)|ϕ=Φ+π].
(48)

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

Hj(ρ,Φ)=2J(sj(Φ))|ρ|.
(49)

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

Appendix C: Derivations of the analytical solution

In this appendix, we show that d22+d32d1, as claimed below (31) in §IV-C. Squaring the integrals in (31), we have:

d12=1π0π0πwj(X)wj(Y)dX dYd22=1π0π0πcos(2X)wj(X)cos(2Y)wj(Y)dX dYd32=1π0π0πsin(2X)wj(X)sin(2Y)wj(Y)dX dY.

In particular,

d22+d32=1π0π0πwj(X)wj(Y)[cos(2X)cos(2Y)+sin(2X)sin(2Y)]dX dY=1π0π0πwj(X)wj(Y)cos(2X2Y)dX dY.

Thus, d22+d32d12 since wj(Φ) ≥ 0 for all Φ and cos(2X − 2Y) ≤ 1 for all X and Y.

Footnotes

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

2Practically speaking, the integral should be restricted to the field of view: x2+y2rmax, but this restriction would complicate analysis by introducing a shift variance into the problem, so we ignore it.

References

1. Fessler JA. Statistical image reconstruction methods for transmission tomography. In: Sonka M, Fitzpatrick J Michael, editors. Handbook of Medical Imaging, Volume 2. Medical Image Processing and Analysis. SPIE; Bellingham: 2000. pp. 1–70.
2. Fessler JA. Penalized weighted least-squares image reconstruction for positron emission tomography. IEEE Trans Med Imag. 1994 Jun;13(2):290–300. [PubMed]
3. Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image reconstruction methods: Space-invariant tomographs. IEEE Trans Im Proc. 1996 Sep;5(9):1346–58. [PubMed]
4. Joseph PM. The influence of modulation transfer function shape on CT image quality. Radiology. 1982 Oct;145(1):179–85. [PubMed]
5. Stayman JW, Fessler JA. Regularization for uniform spatial resolution properties in penalized-likelihood image reconstruction. IEEE Trans Med Imag. 2000 Jun;19(6):601–15. [PubMed]
6. Stayman JW, Fessler JA. Compensation for nonuniform resolution using penalized-likelihood reconstruction in space-variant imaging systems. IEEE Trans Med Imag. 2004 Mar;23(3):269–84. [PubMed]
7. Nuyts J, Fessler JA. A penalized-likelihood image reconstruction method for emission tomography, compared to post-smoothed maximum-likelihood with matched spatial resolution. IEEE Trans Med Imag. 2003 Sep;22(9):1042–52. [PubMed]
8. Qi J, Leahy RM. Resolution and noise properties of MAP reconstruction for fully 3D PET. IEEE Trans Med Imag. 2000 May;19(5):493–506. [PubMed]
9. Shi H, Fessler JA. Quadratic regularization design for 3d cylindrical PET. Proc IEEE Nuc Sci Symp Med Im Conf. 2005;4:2301–5.
10. Fessler JA. Analytical approach to regularization design for isotropic spatial resolution. Proc IEEE Nuc Sci Symp Med Im Conf. 2003;3:2022–6.
11. Shi H, Fessler JA. Quadratic regularization design for fan beam transmission tomography. Proc SPIE 5747, Medical Imaging 2005: Image Proc. 2005:2023–33.
12. Fessler JA. Technical Report 297. Comm and Sign Proc Lab, Dept of EECS, Univ of Michigan; Ann Arbor, MI: Aug, 1995. Resolution properties of regularized image reconstruction methods; pp. 48109–2122.
13. Li Q, Asma E, Qi J, Bading JR, Leahy RM. Accurate estimation of the Fisher information matrix for the PET image reconstruction problem. IEEE Trans Med Imag. 2004 Sep;23(9):1057–64. [PubMed]
14. Clinthorne NH, Pan TS, Chiao PC, Rogers WL, Stamos JA. Preconditioning methods for improved convergence rates in iterative reconstructions. IEEE Trans Med Imag. 1993 Mar;12(1):78–83. [PubMed]
15. Lawson CL, Hanson RJ. Solving least squares problems. Prentice-Hall; 1974.
16. Fessler JA. Matlab tomography toolbox. 2004. Available from http://www.eecs.umich.edu/~fessler.
17. De Man B, Basu S. Distance-driven projection and backprojection. Proc IEEE Nuc Sci Symp Med Im Conf. 2002;3:1477–80.
18. Thibault JB, Sauer K, Bouman C, Hsieh J. A three-dimensional statistical approach to improved image quality for multi-slice helical CT. Med Phys. 2007 Nov;34(11):4526–44. [PubMed]
19. Thibault JB, Bouman CA, Sauer KD, Hsieh J. A recursive filter for noise reduction in statistical iterative tomographic imaging. Proc SPIE 6065, Computational Imaging IV. 2006:60650X.
20. Ahn S, Leahy RM. Spatial resolution properties of nonquadratically regularized image reconstruction for PET. Proc IEEE Intl Symp Biomed Imag. 2006:287–90.