# An analysis of the fan beam CT reconstruction kernel.

AbstractAn analysis of the integral kernels for two well-known fan beam CT reconstruction algorithms has been performed. It was found that the kernel for the standard O([[OMEGA].sup.4]) fan beam tomography reconstruction has an essential bandregion larger than prior numerical analyses had revealed. Here, denotes the essential bandwidth in the reconstructed image. It is also shown that the approximate kernel used in the 0([[OMEGA].sup.3]) convolution-back-projection reconstruction algorithm has an essential bandregion which is significantly different from that of the fan beam data. These calculations imply that properly sampled implementations of the reconstruction algorithms require the fan beam projection data to be available at higher sampling densities than is needed for the recovery of the fan beam data. The fan beam data must either be resampled to the finer mesh before applying the reconstruction algorithms or else acquired on the finer mesh.

Key words and phrases: Fan beam tomography, sampling.

2000 AMS Mathematics Subject Classification--94A20, 65R10, 94A08, 44A12, 92C55.

1 Introduction

In tomographic reconstruction an object is reconstructed from measurements of projections of the object. The maximum resolution attainable in the reconstruction of the object is limited by the discrete nature of the data and by the resolution inherent to the measurement system. Once the measurement system's characteristics are known, it is natural to ask how finely the data must be sampled to achieve the maximum resolution. The resolution in a data set or image can be quantified by the size of the bandregion, or region of support of the Fourier transform. Roughly speaking, a larger bandregion corresponds to a higher resolution and, correspondingly, a finer sampling is required.

In a practical tomographic setting, because the object to be recovered and its projections are compactly supported, the bandregion is unbounded. Instead, one works with the essential bandregion, which, intuitively, is the region where most of the support of the Fourier transform is located. That is, outside the essential bandregion, the Fourier transform is vanishingly small and for all practical purpose can be neglected. A precise treatment of these ideas appears in [7].

The essential bandregion for the fan beam transform of an essentially bandlimited function was determined in [8, 11]. The Petersen-Middleton[12] lattice sampling theory is applied to determine the sampling density necessary to accurately recover the continuous fan beam data from samples. To recover the original function, the fan beam inversion formula, which is an integral transform of the fan beam data, must be applied to the discretely measured data. To determine the sampling density required for an accurate computation from discrete data of this integral transform, the essential bandwidth of the transform kernel must be taken into account. In the definitive paper addressing sampling issues for the fan beam transform[8], it is asserted that numerical experiments indicate that the essential bandregion of the kernel is no bigger than that of the fan beam data. As a consequence, the size of the fan beam data bandregion determines the sampling density necessary for the fan beam reconstruction algorithm.

In this article we will use analytic and numerical techniques to compute the essential bandregion for the kernels used in two well-known fan beam CT reconstruction algorithms[9]. The first reconstruction algorithm is the one appearing in [8]. That algorithm, which will be referred to as the exact algorithm, is of computational complexity O([[OMEGA].sup.4]), where [OMEGA] > 0 denotes the essential bandwidth of the function to be recovered. It will be shown that the essential bandregion of the kernel is larger than that of the fan beam data. This implies that the discrete computation of the continuous integral transform must be performed at a finer sample density than is required for the recovery of the continuous fan beam projection data from its samples. Thus, to properly implement the algorithm, the fan beam data must be available on the finer sample set. This can be accomplished by either resampling the fan beam data to the finer grid required or, alternatively, that interpolation step can be avoided by acquiring the fan beam data on the finer grid.

On the other hand, from our numerical studies it will be apparent that, although non-negligible in the sense of [7], the Fourier transform of the kernel is indeed relatively small on the portion of the kernel essential bandregion not inside the essential bandregion of the data. This explains the reasonable results which are obtained when the algorithm is applied at the sampling density sufficient for the recovery of the fan beam data.

In the second algorithm, referred to as the approximate algorithm, the kernel is modified to rewrite the inversion as a convolution backprojection. This approximation reduces the computational complexity to an attractive O([[OMEGA].sup.3]). However, it will be shown that the essential bandregion of the approximate kernel differs significantly from that of the data. Moreover, unlike the exact kernel, the Fourier transform of the approximate kernel is not small on the portion of the kernel essential bandregion not inside the data bandregion. Thus, the discrete computation of the continuous integral transform necessarily requires either a resampling of the the fan beam projection data or a finer data acquisition.

[FIGURE 1 OMITTED]

The organization of the paper follows. First, the fan beam kernels are presented. Sampling on lattices is reviewed briefly. The essential bandregion for the each of the two kernels is computed. Implications for the implementation of the exact and approximate reconstruction algorithms are discussed. In the last section numerical experiments are presented to illustrate the previous results.

2 The Fan Beam Reconstruction Kernels

ACT scanner provides two-dimensional images of a cross section of an object by measuring integrals of the object's x-ray attenuation through the cross section. For modern CT scanners the x-ray source moves on a circle (the source circle) of radius r > 0 exterior to a circle (scan circle) of radius [rho] > 0 containing the object being imaged. The source angle [beta] measures the position of the x-ray source with respect to x axis. The fan angle a measures the detector position with respect to the ray from the source through the origin. By convention, is positive if, viewed from the source, the detector position is to the left of the central ray. See Figure 1.

It is natural to use the fan beam parameterization for lines in [R.sup.2]. The measurements can be modeled by the fan beam transform of the x-ray attenuation.

The fan beam transform of a function f on [R.sup.2] is

Df([beta], [alpha]) = [f.sup.[infinity].sub.0] f(b + tc)dt, (1)

where b = [(r cos [beta], r sin [beta]).sup.T] and c = [(cos([pi] + [alpha] + [beta]), sin([pi] + [alpha] + [beta])).sup.T]. Df can be regarded as a 2[pi] x [pi]-periodic function.

In a real measurement system only a finite set of data is available. The goal in fan beam tomography is to reconstruct f within [absolute value of x] [less than or equal to] p with as much resolution as possible from discrete samples of Df.

In order to achieve a desired resolution in the reconstructed image, the fan beam transform must be sampled at an appropriate density. A function f [member of] [L.sub.2]([R.sup.n]) is [OMEGA]-bandlimited if its Fourier transform,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (2)

satisfies [??]([xi]) = 0 for [absolute value of [xi]] > [OMEGA]. Such functions can represent details no smaller than [pi]/[OMEGA]. In practice, as the reconstruction region has compact support, the functions to be recovered cannot be bandlimited. As a result, one works with essentially bandlimited functions. A function f on [R.sup.n] is essentially [OMEGA]-bandlimited if [??]([xi]) is small for [absolute value of [xi]] > [OMEGA]. More precisely, it is assumed that [[epsilon].sub.0](f, [OMEGA]) = [[integral].sub.[absolute value of [xi]] > [OMEGA] [absolute value of [??]([xi])d[xi]] is negligible. Refer to [7] for a thorough treatment of these ideas.

Let f [member of] S([R.sup.2]), the Schwartz space of rapidly decaying functions. An essentially [OMEGA]-bandlimited approximate reconstruction of f from fan beam data is achieved by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

where [PHI] is the unit vector associated with the angle [phi] = [beta] + [alpha] - [pi]/2, [v.sub.[OMEGA]] is an [OMEGA]-bandlimited approximation to the ramp filter, and [V.sub.[OMEGA]] is the corresponding [OMEGA]-bandlimited approximate 5 function. That is, for a window [??]([sigma]) which is supported on [absolute value of [sigma]] [less than or equal to] 1, even, and which is approximately one on its support,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (4)

Popular choices for [v.sub.[OMEGA]] include the RamLak kernel, the Shepp-Logan kernel, and the cosine kernel which, for [sigma] [member of] [-1, 1], are given by [??]([sigma]) = 1, [??]([sigma]) = sinc([sigma][pi]/2), and [??]([sigma]) = cos([sigma][pi]/2), respectively.

Let [gamma]([beta]; x) denote the fan angle from the source b to the reconstruction point x = [(s cos [psi], s sin [psi]).sup.T]. Let

[b.sub.[perpendicular to] = (r cos([beta] + [pi]/2), [r sin([beta] + [pi]/2)).sup.T]. (5)

It follows that

cos [gamma] ([beta]; x) = (b - x)b/[absolute value of b - x] [absolute value of b], with [gamma]([beta]; x) </> 0 when [xb.sub.[perpendicular] >/< 0. (6)

The scaling properties of [v.sub.[OMEGA]] together with (6) allow (3) to be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (7)

The straight forward computation of (7) by the trapezoidal rule leads to an O([[OMEGA].sup.4]) algorithm. This implementation is equivalent to the reconstruction given in equation (5.26) of [9]. This algorithm will be referred to as the exact reconstruction algorithm. Accordingly, the exact reconstruction kernel [h.sub.E] is defined by

[h.sub.E]([beta], [alpha]; x) = [[absolute value of b - x].sup.-2][v.sub.[absolute value of b - x][OMEGA]](sin ([gamma]([beta]; x) - [alpha])) cos [alpha]. (8)

In order to improve computational efficiency, Natterer and Wiibbeling[9] suggest an approximation that allows the reconstruction formula to be written as a convolution-backprojection. When [rho] [much less than] r, say [rho]/r [less than or equal to] 1/3, we have [absolute value of x] [much less than] r. Hence, [OMEGA][absolute value of b - x] [approximately equal to] [OMEGA]r, so the bandwidth [OMEGA][absolute value of b - x] of the ramp filter in the kernel can be replaced by [OMEGA]r, which is constant in x. Now the reconstruction formula becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (9)

This integral can be viewed as a convolution-backprojection. Using interpolation to compute values of [v.sub.r[OMEGA]](sin (*)) between lattice points results in an algorithm of complexity O([OMEGA].sup.3]). This algorithm is referred to as the approximate reconstruction algorithm. The corresponding approximate reconstruction kernel [h.sub.A] is given by

[h.sub.A]([beta], [alpha]; x) = [[absolute value of b - r].sup.-2] [v.sub.r[OMEGA]](sin ([gamma]([beta]; x) - [alpha])) cos [alpha]. (10)

2.1 Sampling on Lattices

In order to accurately compute either (7) or (9) by a trapezoidal rule, the integrand must be suitably sampled. To discuss the sampling requirements, it is necessary to review sampling theory on lattices. Given n linearly independent vectors [w.sub.i] [member or] [R.sup.n], denote by W the invertible matrix with columns [w.sub.i]. The lattice [L.sub.W] is defined by

[L.sub.W] = {[W.sub.z] | z [member of] [Z.sup.n]} = W[Z.sup.n]. (11)

Although the lattice [L.sub.W] does not uniquely determine the generating matrix W, [absolute value of det W] is independent of the choice of the generating matrix for the lattice.

The dual lattice [L.sup.[perpendicular to].sub.W] is defined to be [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [W.sup.[perpendicular to]] = 2[pi][([W.sup.-1]).sup.T]. Any non-singular n x n matrix M with integer entries defines a sublattice [L.sub.P] [subset not equal to] [L.sub.W], where P = WM. If [L.sub.P] [subset not equal to] [L.sub.W], then [L.sup.[perpendicular to].sub.P] [contains not equal to] [L.sup.[perpendicular to].sub.W], and [absolute value of det P/det W] = [absolute value of det [W.sup.[perpendicular to]]/det [P.sup.[perpendicular to]]] = [absolute value of det M].

A Q-periodic function f can be sampled on a lattice [L.sub.W] if and only if [L.sub.Q] [subset not equal to] [L.sub.W]. The sample points can be treated as elements of [L.sub.W]/[L.sub.Q]. The Fourier and inverse Fourier transforms in the periodic setting are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (12)

f(x) = [summation over [eta][member of][L.sup.[perpendicular to].sub.Q]) [??]([eta])[e.sup.ix[eta]]. (13)

Definition 1 Let [L.sub.Q] [subset not equal to] [L.sub.W], and f [member of] [L.sub.2]([R.sup.n]/[L.sub.Q]). f is [W.sup.[perpendicular to]]-bandlimited if [??] = 0 outside some compact set K [subset not equal to] [L.sup.[perpendicular to].sub.Q], and (K + [xi]) [intersection] (K + [xi]') = 0 for [xi], [xi]' [member of] [L.sup.[perpendicular to].sub.W], and [xi] [not equal to] [xi]'.

Parseval's formula in the periodic setting follows.

Theorem 2 Let [f.sub.1], [f.sub.2] [member of] [L.sup.2]([R.sup.n]/[L.sub.Q]) be [W.sup.[perpendicular to]-bandlimited and continuous, both with the same bandregion K. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (14)

It is Theorem 2 which specifies the requirements for continuous integrals of bandlimited functions to be computed from sampled data. Note that the Fourier transforms of both [f.sub.1] and [f.sub.2] must vanish outside of K. If the bandregions of [f.sub.1] and [f.sub.2] are [K.sub.1] and [K.sub.2], respectively, then K = [K.sub.1] [union] [K.sub.2] must be used when applying Theorem 2.

In the event that [f.sub.1], [f.sub.2] [member of] [L.sup.2]([R.sup.n]/[L.sub.Q]) are continuous but do not meet the bandlimited hypotheses of Theorem 2, an explicit expression for the trapezoidal rule approximation can be given[8]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

The error in computing the integral using the trapezoidal rule is given by the aliasing, or overlapping of the bandregion of [f.sub.1] with shifts of the bandregion of [f.sub.2]. Estimates on the size of this error are available[7]. For essentially bandlimited functions, this error is indeed small and can be ignored. Hence, for all practical purposes, the hypotheses of Theorem 2 can be broadened to include essentially [W.sup.[perpendicular to]]-bandlimited functions.

3 The Exact Kernel

In this section the essential bandregion, [K.sup.E.sub.[OMEGA]], of the exact fan divergent beam reconstruction kernel will be computed analytically. It will be shown that [K.sup.E.sub.[OMEGA]] is a proper superset of [K.sub.[OMEGA]], the essential bandregion of the fan beam data obtained by projecting an essentially [OMEGA]-bandlimited function with compact support. It will also be shown numerically that, while not completely negligible in the sense of [7], the Fourier transform of the exact divergent beam reconstruction kernel is indeed relatively small on [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]].

Let [B.sub.[rho]] = {x [member of] [R.sup.2] | [absolute value of x] [less than or equal to] [rho]}. For large [OMEGA] it was shown in [8, 11] that the essential bandregion [K.sub.[OMEGA]] of Df for an essentially [OMEGA]-bandlimited f [member of] [C.sup.[infinity].sub.c] ([B.sub.[rho]]) is slightly larger than

[K.sub.[OMEGA] = {(k, 2a) [member of] Z x 2Z | [absolute value of 2a - k] [less than or equal to] [OMEGA]r, [absolute value of k]r [less than or equal to] [absolute value of k - 2a][rho]}. (16)

[K.sub.[OMEGA]] is the black "bowtie" illustrated in Figure 2.

The exact kernel [h.sub.E] is a periodic function on 2[pi]Z x [pi]Z, so its Fourier transform is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

and is defined for (k, 2a) [member of] Z x 2Z. The cos [alpha] multiplication widens the essential bandregion by a negligible amount. Define [[h.bar].sub.E] by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (18)

[FIGURE 2 OMITTED]

Before the Fourier transform of [[h.bar].sub.E] can be computed, the following lemma is needed.

Lemma 3 Let v = s/r and [gamma] be as in (6). Recall that (s, [psi]), resp. (r, [beta]), are the polar coordinates of x, resp. [beta]. Then,

[e.sup.-i2a[gamma]([beta];x)] = [(1 - v[e.sup.i([beta]-[psi])] / 1 - v[e.sup.-i([beta]-[psi])]).sup.a]. (19)

Proof Prom the definition of [gamma] it follows that

sin [gamma]([beta]; x) = v sin ([beta] - [psi]) / [(1 + [v.sup.2] - 2v cos([beta] - [psi])).sup.1/2]. (20)

Substitute (20) into

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (21)

After simplification, we obtain

[e.sup.-i[gamma]([beta];x)] = [(1 - v[e.sup.i([beta]-[psi])] / 1 - v[e.sup.-i([beta]-[psi])]).sup.1/2]. (22)

Raising (22) to the 2a power completes the proof.

The Fourier transform of [[h.bar].sub.E] can now be explicitly computed.

Theorem 4

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (23)

Proof

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

Recall, x = [(s cos [psi], s sin [psi]).sup.T], and the source point b = (r cos [beta], r sin [beta]). Hence,

[absolute value of b - x] = [([s.sup.2] + [r.sup.2] - 2sr cos([beta] - [psi])).sup.1/2]. (25)

Lemma 3 and Graf's generalization of Neumann's addition theorem for Bessel functions [14], equation 11.3(1), together give

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (26)

Substituting (26) into the last line of (24) and evaluating the [beta] integral gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (27)

Theorem 5 Let [absolute value of x] [less than or equal to] [rho] < r. Then, the essential bandregion of the exact kernel [h.sub.E]([beta], [alpha]; x) is slightly larger than the set

[K.sup.E.sub.[OMEGA]] = {(k, 2a) [member of] Z x 2Z / [absolute value of 2a - k] [less than or equal to] [OMEGA]r, [absolute value of k] [less than or equal to] [rho][OMEGA]}. (28)

Proof As mentioned before, the cos [alpha] term has an negligible effect. Hence, the behavior of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] determines the essential bandregion of [h.sub.E].

Debeye's estimates[7] imply that for [J.sub.v](z) decays exponentially fast once v increases slightly beyond z. Hence, for all m in the range of integration in (27), the magnitude of the first Bessel function, [absolute value of [J.sub.k](sm)], is small once [absolute value of k] is slightly larger than s[OMEGA], and the second Bessel function, [absolute value of [J.sub.k-2a](rm)], is small once [absolute value of k - 2a] > r[OMEGA]. Taken together, this implies that [absolute value of [[[??].bar].sub.E] (k; 2a; x)] decays rapidly for (k, 2a) in a set slightly larger than {[absolute value of 2a - k] [less than or equal to] r[OMEGA], [absolute value of k] [less than or equal to] s[OMEGA]}. For s < [rho] each of these sets is a proper subset of the set [K.sup.E.sub.[OMEGA]] arising in the s = [rho] case. Thus, given any [absolute value of x] [less than or equal to] [rho], [[[??].bar].sub.E] (k, 2a, x) is vanishingly small on a set slightly larger than [K.sup.E.sub.[OMEGA]].

Strictly speaking, we have only shown that the essential bandregion of [h.sub.E] is contained in a set slightly larger than [K.sup.E.sub.[OMEGA]]. The possibility that [[[??].bar].sub.E] may be small within [K.sup.E.sub.[OMEGA]] has not yet been eliminated. At first glance, numerical experiments (see Figure 3a) appear to indicate that [K.sup.E.sub.[OMEGA]] is indeed too large, and that [K.sub.[OMEGA]] should be used as the essential bandregion instead of [K.sup.E.sub.[OMEGA]]. [K.sub.[OMEGA]] differs from [K.sup.E.sub.[OMEGA]] in that [K.sup.E.sub.[OMEGA]] contains the finite region

[K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]] = {(k, 2a) [member of] Z x 2Z / [absolute value of k] [less than or equal to] [rho] [OMEGA], [absolute value of k - 2a] [greater than or equal to] [OMEGA]r}. (29)

[K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]] is the gray region illustrated in Figure 2.

Since the sampling requirements for the reconstruction algorithm proposed in [8] rely on the presumption that both [??] and [[??].bar].sub.E] are both small on [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]], the behavior of [[??].bar].sub.E] on this region must be examined. That [??] is small for [absolute value of k - 2a] slightly larger than [OMEGA]r follows from

[[integral].sup.[OMEGA].sub.-[OMEGA]] [J.sub.k] ([rho]m)[J.sub.k-2a] (rm) dm (30)

being small on that region. This is a consequence of the discussion immediately following (3.7) of [8].

A closer look at the numerical examples in Figures 3b and 3c reveals that although the maximum values of the modulus of the discrete Fourier transform of [[h.bar].sub.E] are smaller on [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]] than on [K.sub.[OMEGA]] by about 3%, the modulus of the discrete Fourier transform of [[h.bar].sub.E] does not exhibit rapid decay as (k, 2a) exits [K.sub.[OMEGA]] and proceeds through [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]]. In fact, as [absolute value of k] increases in [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]], the relative maxima increase until [absolute value of k] [approximately equal to] [rho][OMEGA]. Figure 4 shows a log plot of slice through the modulus of the discrete Fourier transform of [[h.bar].sub.E] with [absolute value of x] = 0.95, [OMEGA] = 200, and 2a fixed at r[OMEGA] = 600.

This numerical example illustrates that in a strict sense, [K.sub.[OMEGA]] is not large enough to be taken as the essential bandregion for [[h.bar].sub.E]. [K.sup.E.sub.[OMEGA]] must be used as the essential bandregion. On the other hand, the numerical experiments in Section 6 confirm that for practical purposes, it is not unreasonable to proceed as if [K.sub.[OMEGA]] is the essential bandregion.

[FIGURE 3 OMITTED]

[FIGURE 4 OMITTED]

4 The Approximate Kernel

In this section the essential bandregion [K.sup.A.sub.[OMEGA]] for the approximate kernel [h.sub.A] is determined. The computation of the essential bandregion [K.sup.A.sub.[OMEGA]] is more complicated than that for [K.sup.E.sub.[OMEGA]]. As with [h.sub.E], the cos [alpha] multiplication widens the essential bandregion by a negligible amount. So for the rest of this section, we will work with

[[h.bar].sub.A]([beta], [alpha]; x) = [[absolute value of b - x].sup.-2] [v.sub.r[OMEGA]](sin([gamma]([beta]; x) [alpha])) = h([beta]; x)[q.sub.r[OMEGA]]([beta], [alpha]; x), (31)

where h([beta]; x) = [[absolute value of b - x].sup.-2], and [q.sub.r[OMEGA]]([beta], [alpha]; x) = [v.sub.r[OMEGA]](sin([gamma]([beta]; x) - [alpha])).

The essential support for [[??].bar].sub.A] is found from the convolution in k of [??](k; x) with [??](k, 2a; x).

4.1 The Distance Weighting h([beta]; x)

In polar coordinates, x = [(s cos [psi], s sin [psi]).sup.T] and b = (r cos [beta], r sin [beta]), where r is the fixed source radius, and [beta] is the fan source coordinate. Again, [upsilon] = s/r, where 0 [less than or equal to] s [less than or equal to] [rho] and [theta] = [rho]/r < 1. In practice, for modern scanners 0 is approximately 1/3, and often smaller. In polar coordinates,

h([beta]; x) = 1/[r.sup.2] + [s.sup.2] - 2rs cos([beta] - [psi]). (32)

h is 2[pi]-periodic, so its spectrum is defined on Z and is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (33)

This integral can be evaluated after the substitution z = [e.sup.i[beta]] by the use of residue theory. Alternatively, consider

T(s, [psi]) = ([r.sup.2] - [s.sup.2] [??](k; x) = [r.sup.2] - [s.sup.2]/2[pi] [[integral].sup.2[pi].sub.0] [e.sup.ik[beta]/[r.sup.2] + [s.sup.2] - 2rscos([beta] - [psi])d[beta]. (34)

In [3], pg 451, it is shown that T(s, [psi]) is the unique solution to the differential equation [[nabla].sup.2]T = 0, with boundary condition T(r, [psi]) = [e.sup.-ik[psi]]. Hence,

T(s, [psi]) = [(s/r).sup.k] [e.sup.-ik[psi]] = [[upsilon].sup.k][e.sup.-ik[psi]] (35)

and

[??](k; x) = [([r.sup.2] - [s.sup.2]).sup.-1][(s/r).sup.k] [e.sup.-ik[psi]] = [r.sup.2] [(1 - [v.sup.2]).sup.-1] [v.sup.k][e.sup.-ik[psi]. (36)

Because [??](k; x) decays geometrically in k, uniformly in {[absolute value of x] [less than or equal to] [rho]}, the essential support of [??](k; x) is small. Thus, the weighting of [q.sub.r[OMEGA]]([beta], [alpha]; x) by h([beta]; x) becomes a convolution in the Fourier domain of [??](k, 2a; x) with a function of small essential support in the k direction. The effect is a slight broadening in the k direction of the essential bandregion of [q.sub.r[OMEGA]]([beta], [alpha]; x).

4.2 The Fan Beam Convolution Kernel [v.sub.r[OMEGA]](sin [alpha])

To obtain the bandregion of [q.sub.r[OMEGA]]([beta], [alpha]; x), the Fourier integrals in both the a and [beta] directions must be evaluated. The a integral is handled first.

Recall that [[upsilon].sub.[OMEGA]] is the convolution kernel appearing in parallel beam tomography. It is bandlimited by construction. In the fan beam tomography convolution kernel, [v.sub.r[OMEGA]] appears as with the argument modified by sine. Specifically, the fan beam convolution kernel is [v.sub.r[OMEGA]](sin [alpha]).

For [alpha] [member of] [-[pi]/2, [pi]/2], define [v'.sub.r[OMEGA]] by

[v'.sub.r[OMEGA]]([alpha]) = [v.sub.r[OMEGA](sin [alpha]). (37)

[v'.sub.r[OMEGA]] is extended to R as a periodic function with period [pi].

Lemma 6 [v'.sub.r[OMEGA]] is essentially bandIimited with essential bandregion [-r[OMEGA], r[OMEGA]].

Proof Since [v'.sub.r[OMEGA]] is [pi]-periodic, its spectrum is discrete and defined on 2Z (in contrast with the continuous spectrum for [v.sub.r[OMEGA]]).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (38)

Denote the inner integral by I(a, k).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (39)

The first integral is computed using the well-known expansion [2], (8.514(5)),

cos (k sin [alpha]) = [[infinity].summation over (j=0) [[epsilon].sub.j][J.sub.2j](k) cos (2j[alpha]), (40)

where [[epsilon].sub.0] = 1 and for j [not equal to] 0, [[epsilon].sub.j] = 2.

The second integral of equation (39) is odd with respect to k. As [??] is even with respect to k, and the integral in k of (38) is over a symmetric interval,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (41)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (42)

By [7], III.2.5, [??](2a) decays exponentially for 2a > r[OMEGA] and, hence, [v'.sub.r[OMEGA]]([alpha]) = [v.sub.r[OMEGA]](sin [alpha]) is essentially bandlimited to [-r[OMEGA], r[OMEGA]].

In practice, to implement a parallel beam convolution backprojection, [v.sub.[OMEGA]] is discretely sampled. The spectrum of the sampled [v.sub.[OMEGA]] consists of the spectrum of [v.sub.[OMEGA]] repeated periodically. At the coarsest sampling consistent with the bandwidth [OMEGA], the period is 2[OMEGA], so the right-hand side of the spectrum continues into the periodic repeat of the left-hand side. To reduce the Gibbs phenomenon associated with the non-smoothness of [v.sub.[OMEGA]] at [OMEGA], one can use the Shepp-Logan kernel, which matches both [v.sub.[OMEGA]] and [v'.sub.[OMEGA]] at [+ or -] [OMEGA]. See the upper left corner of Figure 5. The continuity of the derivative of the Shepp-Logan kernel at the cut off (wrapping) frequency reduces the the Gibbs phenomenon artifacts in parallel beam reconstructions. The right figure shows the magnitude of the discrete Fourier transform of the corresponding fan beam kernel, [v'.sub.512]([alpha]), sampled at the same density, again with the periodic repeat of the negative portion of the spectrum plotted shown to the right of the positive portion of the spectrum. Because of the significant oscillation in the Shepp-Logan and Ram-Lak kernels at the cut off (wrapping) frequency, the spectrum is effectively discontinuous there. This will introduce artifacts into the reconstructed images. For the cosine kernel, the oscillation throughout the spectrum is smaller and there is still continuity near the cut off frequency. There will be reduced artifacts in the fan beam reconstructed images.

For fan beam tomography the sine in the argument of the kernel alters the behavior. As seen in the upper right corner of Figure 5, especially for frequencies near [+ or -] r[OMEGA], the behavior of [v'.sub.r[OMEGA]] becomes somewhat erratic. This behavior is explained by equation (41). The continuity across the boundary between repetitions of the spectrum is lost. There appears to be no advantage to using the Shepp-Logan reconstruction kernel over the Ram-Lak kernel (for which the window [eta] is the characteristic function of the interval [-1, 1]), as both exhibit similar behavior. See the middle row of Figure 5. One way to correct for this erratic behavior is to employ a window, such as the cosine window for which [??]([sigma]) drops smoothly to zero [sigma] = [+ or -] [OMEGA]. This is illustrated in the bottom row of Figure 5. Unfortunately, a kernel with this property reduces the high spatial frequency content of the reconstructed image, somewhat lowering reconstruction resolution (in a FWHM sense). However, eliminating the erratic behavior will reduce the artifacts in the reconstructed image.

[FIGURE 5 OMITTED]

4.3 The Behavior of [??](k, 2a; x)

Since [q.sub.r[OMEGA]] is realized as a [beta] dependent shift of the fan beam convolution kernel, the [alpha] integral in the definition of [??] yields in a [beta] dependent phase factor. This phase factor incorporates the entire [beta] dependence. Evaluating the [beta] direction Fourier transform of this phase shift finishes the computation of [??]. Specifically,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (43)

where p([beta], 2a; x), the phase factor, is denoted by

p([beta], 2a; z) = [e.sup.-i[gamma]([beta]; x)2a]. (44)

The essential support for [[??]'.sub.r[OMEGA]](2a) has already been computed. So, it remains to compute the essential support for [??](k, 2a; x).

From Lemma 3 with x = (s cos [psi], s sin [psi]) and [upsilon] = s/r,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (45)

Substituting z = [e.sup.-i[beta]] gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (46)

where C is the unit circle in C, oriented counterclockwise. From a similar substitution with z = [e.sup.i[beta]], we obtain

[??](k, 2a; x) = [??](-k, -2a; x). (47)

To study the behavior of [??](k, 2a; x), the lattice Z x 2Z is divided into four regions. For a [greater than or equal to] 0, the regions are

1. a = 0,

2. 0 < a < k,

3. 0 < a and k [less than or equal to] 0,

4. 0 < a and 0 [less than or equal to] k [less than or equal to] a,

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

Regions 2, 3, and 4 are extended to a < 0 by reflection through the origin. See Figure 6. By (47) it suffices to look at only one of the two components of each region.

For region 1 it follows immediately from from (46) that

[??](k, 0; z) = [e.sup.-ik[psi]] [[delta].sub.k0]. (48)

For region 2, when a > 0 and k - a > 0, the integrand of (46) has no poles in the interior of the unit disk and, therefore,

[??](k, 2a; x) = 0. (49)

Region 3 is the first interesting case is. Here, a > 0 and k [less than or equal to] 0. Equation (47) permits us to work with the equivalent a < 0 and k [greater than or equal to] 0.

So now assume a < 0 and set n = -a. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (50)

As k [greater than or equal to] 0 and n > 0, the pole at z = v is the only pole. To evaluate (50) it is necessary only to compute the residue at z = v. Since z = v + (z - v),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (51)

Also,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (52)

Hence, the integrand of (50) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (53)

For the residue at z = v only the terms where j - l = -1 are relevant. Denote by R(k, n) the residue at v. Select these terms and apply equations (4.1.3) and (4.3.2) from [13]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (54)

Equation (4.22.1) of [13] was used for the last step in equation (54) and [sub.2][F.sub.1] denotes Gauss' hypergeometric function, (4.21.3) in [13].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (55)

Therefore, on region 3,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (56)

For region 4 we have 0 [less than or equal to] k [less than or equal to] a. Inside of C, the integrand of (46) has a pole only at z = 0; so to compute the residue, [(z - v/1 - vz).sup.a] is expanded in powers of z.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (57)

Reindex, letting n = l + j, and interchange the order of summation. Then,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (58)

Since a [greater than or equal to] a - k [greater than or equal to] 0, the residue of the integrand in (46) is the coefficient of [z.sup.a-k] in (58). Hence,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (59)

where equation (4.22.1) of [13] was used for the last step.

As a side comment, note that both (56) and (59) simultaneously hold for k = 0. Accordingly, the identity

[P.sup.(0,-1).sub.a] (1- 2[v.sup.2]) = (1 - [v.sup.2]) [P.sup.(0,1).sub.a-1] (1- 2[v.sup.2]) (60)

must hold. Indeed, this identity follows directly from equations (4.1.3) and (4.22.2) in [13].

The following theorem[4] will be used to explicitly compute the decay of [??](k, 2a; x).

Theorem 7 Let [mu] > -1, 0 < [theta] < 1, [alpha] [greater than or equal to] -1, [beta] [greater than or equal to] -1, and [theta] < [mu]/[mu]+2. Then, for n [member of] N, as n [right arrow] [infinity],

1. when -1 < [mu] < 2[theta]/1-[theta], [[theta].sup.[mu]n][P.sup.([alpha] + [mu]n, [beta]).sub.n] (1 - 2[[theta].sup.2] = O ([n.sup.-1/2]) (61)

2. when [mu] > 2v/1 - v, [v.sup.[mu]n][P.sup.([alpha] + [mu]n, [beta]).sub.n] (1 - [2v.sup.2] = C [n.sup.-1/2] [t.sup.n] (1 + O ([n.sup.-1])), (62)

where for fixed v, C = C([mu]) is uniformly bounded in [mu] > 2v/1 - v, and t = t([mu], v) < 1, with t [right arrow] 1 as [mu] [right arrow] 2v/1 - v.

Lemma 8 For 0 [less than or equal to] < v < 1, define [K.sub.v] by

[K.sup.v] = {(k, 2a) [member of Z x 2Z | [absolute value of k] [greater than or equal to [absolute value of k - 2a]v}. (63)

1. [??](k, 2a; x) = O ([k.sup.1/2]) as (k, 2a) recedes along the rays k = 2av emanating from the origin within [K.sub.v].

2. [??](k, 2a; x) decays exponentially in v as (k, 2a) recedes along the rays k = 2av emanating from the origin outside [K.sup.[upsilon]].

Proof On region 2 and on region 1 away from the origin [??](k, 2a; x) is identically 0. Moreover, both regions (excluding the origin) are outside of [K.sub.v].

On the a > 0 portion of region 4, for any lattice point (k, 2a), k = [mu]a, with 0 [less than or equal to] [mu] [less than or equal to] 1. Consider [??], along the ray determined by [mu]:

[absolute value of [??](k, 2a; x)] = [absolute value of [??]([mu]a, 2a; x)] = [absolute value of [v.sup.[mu]a] [P.sup.([mu]a, -1).sub.a(1 - [mu]) (1 - [2v.sup.2])]. (64)

Let v = a(1 - [mu]). Note that v [member of] N. Then,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (65)

which by Theorem 7 decays exponentially in v when [mu]/1 - [mu] > 2v/1 - v, or equivalently [mu] > 2v/1 + v. For [mu]/1 - [mu] < 2v/1 - v (equivalently [mu] < 2v/1 + v], [absolute value of [??](k, 2a; x)] decays like O([v.sup.-1/2]). That is, for (k, 2a) in region 4 with a > 0, [absolute value of [??](k, 2a;x)] decays exponentially when receding from the origin with k > v/1 + v (2a), and decays slowly otherwise. We remark that as [mu] increases away from the critical value 2[theta]/1 + [theta] the rate of decay increases, and that restricting [mu] to [mu] [greater than or equal to] [[mu].sub.0] > 2[theta]/1 + [theta] allows for a uniform estimate on the decay rate to be obtained.

Observe that [mu] < 2v/1 + v within the first quadrant part of region 4 fills the the interior of the first quadrant portion of [K.sub.v]. Correspondingly, [mu] > 2v/1 + v within the first quadrant part of region 4 together with the first quadrant part of regions 1 and 2 fill the exterior of [K.sub.v] within the first quadrant. The boundary ray, k = v/1 + v (2a), coincides with one of the boundary segments of [K.sub.v]. The behavior in the third quadrant follows by reflection through the origin.

Restrict attention now to the fourth quadrant part of region 3, where a > 0 and k [less than or equal to] 0. For any lattice point (k, 2a) in this region, k = [mu](-a), where 0 [less than or equal to] [mu]. Consider [??] along this ray. With n = -a,

[absolute value of [??](k, 2a; x)] = [absolute value of [??]([mu]n, 2a; x)] = [absolute value of [v.sup.[mu]n](1 - [v.sup.2]) P.sup.([mu]n, 1).sub.n-1] (1 - 2[v.sup.2])], (66)

which by Theorem 7 exhibits exponential decay in n for [mu] > 2v/1 - v, and O([n.sup.-1/2]) decay for [mu] < 2v/1 - v. That is, for (k, 2a) in portion of region 3, with a < 0, [absolute value of [??](k, 2a; x)] decays exponentially when receding from the origin with k > - v/1 - v (2a), and decays slowly otherwise. Again we note that a uniform estimate on the exponential rate of decay can be obtained by setting a lower bound for [mu] strictly larger than 2v/1 - v. The region with exponential decay fills the exterior of [K.sub.v] within the fourth quadrant, and the slow decay region fills the interior. Again, the boundary ray, k = - v/1 - v(2a), coincides with the boundary segment of [K.sup.v].sub.[OMEGA]]. Reflecting onto the second quadrant completes the picture.

Recall that [??](k, 2a; x) = [[??]'.sub.r[OMEGA]](2a)[??](k, 2a). It was just shown that the factor [??](k, 2a; x) decays outside of of {[absolute value of k] [less than or equal to] [absolute value of k - 2a]v}.

Theorem 9 For all [absolute value of x] [less than or equal to] [rho], [[??].sub.A](k, 2a) decays exponentially as the distance from the origin of (k, 2a) increases outside a set slightly larger than [K.sup.A.sub.[OMEGA]],

[K.sup.A.sub.[OMEGA]] = {(k, 2a) [member of] Z x 2Z | [absolute value of 2a] [less than or equal to] [OMEGA]r, | k [absolute value of k - 2a][rho]}, (67)

Proof For all x in the scan circle [absolute value of x] [less than or equal to] [rho] we have v [less than or equal to] [theta]. Hence, for any reconstruction point x, [??](k, 2a; x) will decay on the complement of the union over all [K.sub.v] with v [less than or equal to] [theta]. Since [K.sub.[theta]] = [[universal].sub.0[less than or equal to]v[less than or equal to][theta]][K.sub.v], [??](k, 2a; x) will decay on the complement of [K.sub.[theta]]. For all x in the scan circle [??](k, 2a; x) decays slowly along some rays inside [K.sub.[theta]].

Recall that [??](k, 2a; x) = [[??]'.sub.r[OMEGA]](2a)[??](k, 2a; x). In Section 4.2 it was shown that [[??]'.sub.r[OMEGA]](2a) has essential support on [absolute value of 2a] [less than or equal to] r[OMEGA]. Together, these imply that [??](k, 2a; x) decays exponentially on a set slightly larger than the complement of [K.sub.A.sub.[OMEGA] = [K.sub.[theta]] [intersection] {(k, 2a)| [absolute value of 2a] [less than or equal to] r[OMEGA]}.

Finally, it was observed that the multiplications in (31) of [q.sub.r[OMEGA]] by cos [alpha] and by the distance weighting h([beta]; x) have the effect of broadening the essential support of [[[??].bar].sub.A] by a negligible amount. Thus, the essential bandregion of [[[??].bar].sub.A] is given by [K.sup.A.sub.[OMEGA]].

Numerical studies (see Figure 7) confirm the behavior predicted in Theorem 9. In particular, the rapid decay of [[[??].bar].sub.A] exterior to the entire "bowtie" is apparent. Compare with the behavior of [[[??].bar].sub.E] on [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]] in Figure 3.

5 Implications

In the previous two sections, the essential bandregions [K.sup.E.sub.[OMEGA]] and [K.sup.A.sub.[OMEGA]] of the kernels were shown to be different from the essential bandregion [K.sub.[OMEGA]] of the fan beam projection data. In this section, the implications for implementation of the exact and approximate reconstruction algorithms are explored.

Both the exact and approximate algorithms require the computation by a trapezoidal rule of continuous integrals from discrete data. For the exact algorithm (7) is computed and for the approximate algorithm the equation is (9). In either case Theorem 2 is used. Because the essential bandregions of the kernel and the fan beam data are not the same, the union of the two essential bandregions must used for the set K in Theorem 2. Should the union be larger than than [K.sub.[OMEGA]] (as is indeed the case here), the sample density needed to apply the trapezoidal rule will be finer than is required to accurately recover the projection data itself. Since the kernels themselves are known analytically, obtaining them at the higher density needed for the trapezoidal rule is easy. Handling the projection data is more difficult. Either the projection data must be resampled via interpolation, or else the projection data can itself be acquired at the higher density. It is important to note that we are not asserting that the sampling requirements from [8] are insufficient, but rather that the implementation of the reconstruction algorithms require the data to be available at a finer density.

In the reference implementation of the exact algorithm[8], it is asserted (with slightly different notation) that the essential support for [h.sub.E] is [K.sub.[OMEGA]], the essential support of the fan beam data. A consequence of this assertion is that the sample density sufficient for recovering the fan beam data [D.sub.f] should be sufficient for use in the implementation of the exact reconstruction algorithm.

In particular, [V.sub.[OMEGA]]* f is computed from discrete samples of [D.sub.f] on a lattice [L.sub.w]. The sample lattice [L.sub.w] is determined by requirements on the dual lattice [L.sup.[perpendicular to].sub.W]. Theorem 2 with K = [K.sub.[OMEGA]] is applied to (7). The translates of Ka by non-zero elements of the dual lattice [L.sup.[perpendicular to].sub.W] must not overlap [K.sub.[OMEGA]]. This can be accomplished by selecting the matrix [W.sup.[perpendicular to]] to generate the dual lattice as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (68)

where

[[xi].sub.1] > 2[rho]r[OMEGA]/r + [rho], [[xi].sub.2] > 2r[OMEGA]. (69)

This translates to the sampling requirements

[DELTA][beta] < r + [rho]/[rho] [pi]/r[OMEGA], [DELTA][alpha] < [pi]/r[OMEGA] (70)

which will be referred to as the "standard sampling requirements."

However, in Section 3 it was shown that [absolute value of [[??].sub.E]] does not effectively vanish on [K.sub.E.sub.[OMEGA] - [K.sub.[OMEGA]. Ignoring this leads to aliasing artifacts in the reconstruction. That is, applying the trapezoidal rule at the sample density suggested in [8] will not give an alias-free recovery of [V.sub.[OMEGA]* f. To properly recover [V.sub.[OMEGA]* f, the dual lattice must be chosen so that translations of [K.sub.E.sub.[OMEGA] = [K.sub.[OMEGA] [union] [K.sub.[OMEGA] by non-zero elements do not overlap. To properly avoid overlaps, the bound on [[xi].sub.1] of (69) must be replaced by

[[xi].sub.1] > 2[rho][OMEGA]. (71)

This in turn implies, at least theoretically, that to apply the trapezoidal rule, the sampling density in the [beta] direction must be increased by a factor of 2[rho]/2[rho]r/r + [rho] = 1 + [theta] to avoid aliasing artifacts arising from the non-vanishing of [[[??].bar].sub.E] on the region [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]]. No resampling is required in the a direction.

On the other hand, from a practical perspective, although not effectively vanishing in the sense of [7], [absolute value of [[[??].bar].sub.E]] is still relatively small on [K.sup.E.sub.[OMEGA]] - [K.sub.[OMEGA]] as compared with its values on Ka, especially when [theta] [much less than] 1. Thus, it is expected that the induced artifacts would be weak compared to the ideal reconstruction.

For the approximate algorithm the computational sample density is even finer than that for the exact algorithm. In this case [(a is not a subset of the kernel essential bandregion [K.sup.A.sub.[OMEGA]]. The set [K.sub.[OMEGA]] [universal] [K.sup.A.sub.[OMEGA]] is larger than either component in both the k and 2a directions. See Figure 8. As a result, the fan beam data must be resampled in both the [beta] and [alpha] directions. The shifts of the dual lattice to avoid overlaps of [K.sub.[OMEGA]] [universal] [K.sup.A.sub.[OMEGA]] are given by the following bounds:

[[xi].sub.1] > 2[r.sup.2][rho][OMEGA]/[r.sup.2] - [[rho].sup.2], [[xi].sub.2] > (2r + [rho])[OMEGA]. (72)

This corresponds to

[DELTA][beta] < [r.sup.2] - [[rho].sup.2]/r[rho] [pi]/r[OMEGA], [DELTA][alpha] < 2[pi]/(2r + [rho])[OMEGA] (73)

which will be referred to as the "extra fine requirements."

It is interesting to note that for each of [K.sub.[OMEGA]] and [K.sup.A.sub.[OMEGA]] individually, a vertical shift of [[xi].sub.2] > 2r[OMEGA] is sufficient to avoid overlapping, but for [K.sub.[OMEGA]] and [K.sup.A.sub.[OMEGA]], that shift is not sufficient. Compared with the sampling requirements for the fan beam data, the working density requirements are a factor of 2[r.sup.2][rho]/[r.sup.2] - [[rho].sup.2] - [[rho].sup.2] /2[rho]r/r+p = [(1 - [theta]).sup.-1] finer in the [beta] direction and a factor of (2r + [rho])/(2r) = 1 + [theta]/2 in the [alpha] direction.

6 Numerical Experiments

In this section we report on the results of numerical experiments performed to illustrate the results of the previous sections.

Unless otherwise indicated, for the numerical experiments shown here the scan radius is [rho] = 1 and the source radius is r = 3. The RamLak kernel was used as it exhibits the sharpest roll off in Fourier space and better highlights artifacts in the reconstructed images. The exact and approximate kernels [[h.bar].sub.E]([beta], [alpha]; x), resp. [[h.bar].sub.A]([beta], [alpha]; x), were computed by sampling the analytic expressions, equations (18) and (31) with [OMEGA] = 200 and [absolute value of x] = 0.95. The standard, resp., extrafine sampling distances are the coarsest rates consistent with equations (70), resp., (73) evaluated at ft = 200. Figure 9 shows the modulus of the discrete Fourier transforms of the exact and approximate kernels sampled at the standard rate. Sampling the kernels induces the spectrum to repeat periodically. As predicted theoretically, when the exact kernel is sampled at the standard rate, the periodicity does not introduce any apparent overlap, indicating that the exact kernel has been adequately sampled.

[FIGURE 8 OMITTED]

[FIGURE 9 OMITTED]

On the other hand, from the rightmost image of Figure 9, it is apparent that the sample rate appropriate for the exact kernel (and Df) is too coarse for the approximate kernel. There will be aliasing in the high spatial frequency components in the reconstructions. The aliasing can be alleviated by increasing the sampling density. Although increasing the sampling density only in the [beta] direction will eliminate the overlap seen in Figure 9, it is still not sufficient to accurately compute equation (9), as the sampling must be dense enough to ensure that the periodic repeats of [K.sub.[OMEGA]] [union] [K.sup.A.sub.[OMEGA]] do not overlap.

To observe the the aliasing this undersampling introduces, consider an ideal bandlimiting phantom approximating a delta function, centered at (0.50, 0) in a scan circle of radius 1, scaled to have a maximum value of 1 at its center:

[f.sub.[OMEGA](x) = 2 [J.sub.1]([OMEGA][absolute value of x - (0.5, 0)]) / ([OMEGA][absolute value of x - (0.5, 0)]. (74)

This phantom, centered at the origin, was used by Kruse[5] in his studies of resolution in CT. The Fourier transform of the phantom is simply a phase shifted characteristic function of the disk of radius [OMEGA],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. (75)

Although this function is not spatially limited, for x outside of [rho] = 1 and for [OMEGA] large, [f.sub.[OMEGA]] (x) can effectively be ignored. This phantom is indeed bandlimited, its projections can be computed analytically, and it can be used to demonstrate the Fourier properties of the fan beam transform.

Using the projection slice theorem, [Df.sub.200] was computed analytically, and sampled at the standard rate for [OMEGA] = 200. The modulus of the discrete Fourier transform of the sampled [Df.sub.200] is shown in Figure 10 to the right of the modulus of the sampled exact kernel. For both the essential bandregion is contained within the periodically extended [K.sub.[OMEGA]]. Hence, for a reconstruction with the exact kernel, Theorem 2 can be applied with essentially no error. Figure 11 shows the phantom [f.sub.200] on the upper left, and the reconstruction using the exact kernel with [OMEGA] = 200 in the upper right. The bottom row is a plot of the difference between the two. The reconstructions were computed on a 257 x 257 grid, twice the Nyquist rate in each linear direction. The peak value of the phantom corresponds to attenuation of water, a level of zero CT units and the zero of the phantom, corresponding to air, is assigned a level of-1000 CT units. The phantom and reconstructions are shown with a window of 20 CT units, about -1000 CT unit background. In terms of the original phantom, this corresponds to a window of [-0.01, 0.01]. The maximum difference between the two is about 1.2 CT units.

[FIGURE 10 OMITTED]

For the approximate kernel, the standard sampling rate is not sufficient. The error term, the last term in equation (15), will be non-zero. In particular, the portion of the data spectrum in the upper left and lower right corners (corresponding to the regions of overlap in kernel spectrum) will appear in the error term. See the top row of Figure 12. The introduced aliasing error can be removed by use of the finer sampling set. This is illustrated in the bottom row of Figure 12, where the modulus of the discrete Fourier transform of the sampled approximate kernel with x = (0.95, 0) and the sampled data set are displayed next to each other. The top row was computed at the standard sampling rate, while the bottom row was computed on the extra fine grid. The error arising from bandregion overlap will be larger at reconstruction points further from the origin. This happens because the kernel bandregion is narrower for reconstruction points closer to the origin. Indeed, for points close to the origin, the standard sampling grid is fine enough to avoid aliasing error. See Figure 13

When reconstructing with the approximate kernel, artifacts may arise from several sources. If the working reconstruction sampling density is not fine enough, then aliasing errors will result. A related artifact source is the interpolation step in an O([[OMEGA].sup.3]) approximate algorithm, especially when there is undersampling in the a direction. Even when the sampling density is fine enough, and an O([[OMEGA].sup.4]) algorithm is used to eliminate any interpolation error, there remains an underlying reconstruction error stemming from the replacement of exact kernel with the approximate kernel. This is illustrated in the table of approximate reconstructions shown in Figure 13. Each of the reconstructions shown in the top row was performed at the standard sampling density; those in the the bottom row were performed with samples on the extra fine grid. The upper left image shows the result of applying the fast O([[OMEGA].sup.3]) algorithm at the standard sampling density. The artifact at the far left arises from undersampling. The two arcs on the right, emanating from the object center, appear to arise from the undersampled interpolation. They are missing from the upper left image, obtained by applying the slow O([[OMEGA].sup.4]) algorithm, which uses the approximate kernel, but does not use interpolation, at the standard sampling density. The width of the left edge artifact is reduced, though the circular artifact in the center appears slightly stronger. In both images a faint mirror image of the phantom can be seen. The bottom row shows reconstructions performed at the finer density. The bottom left image shows faint residual undersampling artifacts. These are probably due to a slight undersampling in the interpolation which reflects the non-sharp cutoff at the edge of the kernel bandregion. Slightly increasing the sampling density beyond the extra fine density proposed here eliminates that residual artifact. The lower right image was obtained by applying the slow O([[OMEGA].sup.4]) algorithm with the approximate kernel at the extra fine density. All interpolation artifacts, and essentially all undersampling artifacts, have been eliminated. The remaining artifact appears to be the consequence of replacing the exact kernel with the approximate kernel. This image is essentially the best reconstruction possible with the approximate kernel. Oversampling beyond the fine grid proposed here does not improve the reconstruction.

[FIGURE 11 OMITTED]

[FIGURE 12 OMITTED]

[FIGURE 13 OMITTED]

Despite the presence of the artifacts seen in the Figure 14, in most practical situations the fast O([[OMEGA].sup.3]) algorithm applied at the standard sampling density results in reconstructions in which all main features are faithfully recovered, with most of the errors confined to high spatial frequency components. To understand the reason for this phenomenon, consider Figure 15 in which the modulus of the discrete Fourier transform of [[h.bar].sub.E] - [[h.bar].sub.A] is shown for [absolute value of x] = 0.95 at both the standard and extra fine sampling densities. There is remarkable agreement on [K.sub.[OMEGA]] [intersection] [K.sup.A.sub.[OMEGA]]. This indicates that, except for the high spatial frequency content in the region outside of [K.sub.[OMEGA]] [intersection] [K.sup.A.sub.[OMEGA]], [[h.bar].sub.A] is indeed a good substitute for [[h.bar].sub.E], and confirms that the errors arising in replacing [[h.bar].sub.E] by [[h.bar].sub.A] induces high frequency errors. It should be noted, however, that the agreement in the low spatial frequency domain is not perfect and a residual, low-level disagreement between the two kernels remains.

[FIGURE 14 OMITTED]

[FIGURE 15 OMITTED]

References

[1] M. Abramowitz, I. Stegun, Eds., Handbook of mathematical functions with formulas, graphs, and mathematical tables, United States Department of Commerce, U.S. Government Printing Office, 1972.

[2] I.S. Gradshteyn, I.M. Ryzhik, Tables of Integrals, Series, and Products, 4th. ed., Academic Press, New York, 1980.

[3] F. Hildebrand, Advanced Calculus for applications, 2nd. ed., Prentice Hall, Englwood Cliffs, New Jersey, 1976.

[4] S. Izen, Refined estimates of Jacobi polynomials, J. Approx. Theory, Academic Press, in press.

[5] H. Kruse, Resolution of reconstruction methods in computerized tomography, SIAM J. Sci. Stat. Comput. 10(1989), pp, 447-474.

[6] A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, SIAM, Philadelphia, PA, 2001.

[7] F. Natterer, The Mathematics of Computerized Tomography, Wiley, New York, 1986.

[8] F. Natterer, Sampling in Fan Beam Tomography, SIAM J. Appl. Math, 53 (1993), pp. 358-380.

[9] F. Natterer and F. Wubbeling, Mathematical Methods in Image Reconstruction, SIAM, Philadelphia, PA, 2001.

[10] F. Olver, Asymptotics and Special Functions, A.K. Peters, Natick, MA, 1997.

[11] V. Palamodov, Localization of harmonic decomposition of the Radon transform, Inverse Problems, 11 (1995), pp. 1025-1030.

[12] D. P. Petersen and D. Middleton, Sampling and Reconstruction of wave number-limited functions in n-dimensional Euclidean spaces, Inform. Contr., 5 (1962), pp. 279-323.

[13] G. Szego, Orthogonal Polynomials, AMS, Providence, RI, 1975.

[14] C. Watson, A Treatise on the Theory of Bessel Functions, 2nd ed., Cambridge University Press, Cambridge, 1952.

Steven H. Izen Department of Mathematics, Case Western Reserve University, Cleveland, Ohio shi@cwru.edu

Printer friendly Cite/link Email Feedback | |

Author: | Izen, Steven H. |
---|---|

Publication: | Sampling Theory in Signal and Image Processing |

Geographic Code: | 1USA |

Date: | Sep 1, 2008 |

Words: | 9727 |

Previous Article: | A memorial tribute to Isao Someya, 1915-2007. |

Next Article: | Error of sinc approximation of analytic functions on an interval. |

Topics: |