High-Resolution Computed Tomography

This paper investigates high-resolution reconstructions from eeciently sampled data in parallel-beam tomography, in particular local tomography. A class of sampling schemes is deened and characterized, and it is shown that the standard scheme and the interlaced scheme of Cormack are most promising. An error analysis for the ltered backprojection algorithm for both global and local tomography is presented. The analysis provides insights on how to realize the theoretically superior resolution of the interlaced scheme in practice. A numerical experiment with real data demonstrates the feasibility of high-resolution local tomography using the interlaced scheme.


Introduction
Computed tomography (CT) entails the reconstruction of a function f from a nite number of line integrals of f.Recent applications include microtomography, where very small structures are investigated and the highest possible resolution is demanded.This paper addresses the question of how to sample the data e ciently and reconstruct from e ciently sampled data in order to achieve maximum resolution.It continues a line of inquiry begun by Cormack 2] and further developed by Lindgren and Rattey 24], Natterer 19,20], Kruse 14], Klaverkamp 13], and one of the authors 5,4,8,9].A class of potentially reasonable sampling schemes will be considered, and it will be shown that two schemes, the standard scheme and the interlaced scheme of Cormack, are most e cient.The interlaced scheme in turn allows in principle for almost twice the resolution as the standard scheme.It will be investigated how much of this advantage can be realized in practice, using the ltered backprojection algorithm.This is done by means of a detailed error analysis which yields speci c guidance on the choice of the parameters in the algorithm and explains the advantages and drawbacks of using the interlaced lattice.
Ordinary tomography is global inasmuch as reconstruction at a point requires integrals over lines far from that point.Local tomography as originally developed in 27,29,30] Dept. of Mathematics, Oregon State University, Corvallis, OR 97331.faridani@math.orst.eduy Dept. of Physiology and Biophysics, Mayo Clinic, 200 First Street S.W., Rochester, MN 55905.This work was supported NIH grant R01 RR 11800-4 and by NSF grant DMS-9803352.uses only integrals over lines close to the point of interest, but reconstructs a related function Lf = f + ?1 f instead of the function f itself.However, images of Lf look similar to images of f and, e.g., preserve boundaries and the direction of density jumps across boundaries 6].In this paper we are primarily interested in local reconstructions from e ciently sampled data.For recent developments in local tomography see, e.g., 7,8,12,15,23,25,26].
The paper is organized as follows.In the next section some notation is introduced and standard facts about the Radon transform are reviewed.Section 3 contains a characterization of a class of sampling schemes which we call 'admissible sampling lattices'.It is shown that among all admissible lattices the commonly used standard lattice and the so-called interlaced lattice rst introduced by Cormack are the only ones which fully exploit a certain symmetry of the Radon transform and allow to limit the angular range of directions of exposure to 180 degrees.The following section brie y describes the implementation of the ltered backprojection algorithm for data sampled on an arbitrary admissible sampling lattice.Section 5 contains the main results of this paper.We review the sampling conditions for the standard and interlaced lattices which are derived by using Shannon sampling theory and give the conditions so that an essentially bandlimited function can be determined reliably by the sampled data.In order to obtain insight on how to obtain good recontructions from properly sampled data we derive an error estimate for the ltered backprojection algorithm.This estimate extends previous work 5] inasmuch as it applies to local tomography as well as to global tomography.The estimates show that if the parameters of the algorithm are chosen properly it should be possible to realize the higher theoretical resolution of the interlaced lattice in practice.The last section contains a numerical experiment with real data using local tomography.It demonstrates that the interlaced scheme can be successfully used in practice for high resolution tomography.

Standard notation and de nitions
Let Z Z; IR; C denote the integers, real and complex numbers, respectively.The characteristic function M of a set M is de ned by M (x) = 1 for x 2 M and M (x) = 0 otherwise.
The two-dimensional Radon transform maps a density function f into its line integrals.Throughout this paper we will assume that f 2 C 1 0 ( ), i.e., f is in nitely di erentiable and vanishes outside the unit disk of IR 2 .This assumption simpli es the mathematical proofs, and although the density functions occurring in practice are not necessarily smooth, we will see that our theoretical results describe the phenomena observed in practice well.
Let = (cos '; sin ') be the unit vector in IR 2 with polar angle ', and ?= (?sin '; cos ').For f 2 C i.e., Rf('; s) is the integral of f over the line in direction ?with signed distance s from the origin.Sometimes Rf is considered as a function of s for xed '.In this case we write R ' f(s) for Rf('; s).
The Fourier transform of a function g 2 C 1 0 (IR n ) is de ned by ĝ( ) = (2 ) ?n=2 Z IR n g(x)e ?i<x; > dx and is extended to larger classes of functions or distributions by continuity or duality.Here < x; > denotes the usual inner product in IR n .
In particular, the Fourier transform of R ' f is given by (R ' f) ^( ) = (2 ) ?1=2 Z IR R ' f(s)e ?is ds: The following well-known relation between the Fourier transforms of R ' f and f is very useful: (2) Equation ( 2) is called the projection-slice theorem.It can be used to prove the following approximate inversion formula: Theorem 2.1 Let e 2 L 2 (IR 2 ) be a radial function such that j j 1=2 ê( ) 2 L 2 (IR 2 ), and the even function of one variable given by ê( ) = (2 ) ?1 (j j).Let the associated convolution kernel k be given by k( ) = 1 2 (2 ) ?3=2 j j ( ).Then e f(x) = Z 2 0 Z IR k(< x; > ?s)Rf('; s) ds d': ( Proof: The relation (3) can be veri ed by writing e f as e f(x) = R ê( ) f( )e i<x; > d , expressing the integral in polar coordinates, and using the relation (2); see, e.g., 8].The formula is also valid under more general conditions 16,18,27].
If e is an approximate -function, (3) gives an approximate reconstruction formula for f.In this case the convolution kernel cannot have compact support.This means that in order to reconstruct e f(x) one needs integrals over lines far from the point x, i.e., the reconstruction is not local.A local reconstruction formula is obtained by letting e = e 0 , with e 0 an approximate -function and Calder on's operator given by ( g) ^( ) = j jĝ( ): Then, because of ( e 0 ) f = e 0 ( f) the reconstruction is an approximation for f.One can show 6] that images of f show the same boundaries and the same direction of jumps across boundaries as images of f, but f is cupped in regions where f is constant.This cupping can be somewhat ameliorated by reconstructing a linear combination of f and ?1 f 7].Using , we may write the relation between the point-spread function e and the convolution kernel k as k = (4 ) ?1 R ' e. Concrete examples for k and e 0 can be found in 6].
Discretizing equation (3) yields the ltered backprojection algorithm, the most popular reconstruction method.In order to describe it, we need to discuss how the Radon transform is sampled.

Sampling the Radon Transform
From (1) we see that Rf is a function with domain 0; 2 ) IR.The subsequent analysis of sampling and resolution will make use of Fourier analysis.This requires both the domain of Rf as well as the sampling sets to have a group structure.Equipped with addition modulo 2 the interval 0; 2 ) becomes a group, called the circle group, which we denote by T. Then the domain of Rf may be identi ed with the group T IR.The addition on T IR can be viewed as the usual addition in IR 2 but modulo 2 in the rst component.
The task of tomography is to reconstruct f from nitely many measurements of Rf.
In the parallel-beam sampling geometry, a set of angles f' j ; j = 0; : : :; P ?1g is selected and for each angle ' j a number line integrals Rf(' j ; s jl ) are measured.We require the set f' j ; s jl g of all points where Rf is measured to be a subgroup of T IR, and for practical reasons there should be more than one measured line for each occurring angle ' j .A sampling set satisfying these two requirements is called an admissible sampling lattice.The admissible lattices may be parameterized as follows: Lemma 3.1 Let L be an admissible sampling lattice.Then there is d > 0 and integers N; P, such that 0 N < P and L = L(d;N;P) = f(' j ; s jl ) : ' j = 2 j=P; s jl = d(l + jN=P); j = 0; : : : ; P ?1; l 2 Z Zg : (4) We defer a formal de nition of admissible sampling lattices and the proof of Lemma 3.1 to the end of this section.>From (4) we see that there are P directions of exposure corresponding to the equidistant angles ' j = 2 j=P.For each direction integrals over an equidistant set of lines with spacing d are measured.This collection of equidistant parallel lines is shifted by an amount djN=P which varies with the angle ' j .
The most important lattices are the standard lattice L S = f(' j ; s l ) : ' j = 2 j=P; s l = d l; j = 0; : : : ; P ?1; l 2 Z Zg which is obtained by letting N = 0, and the interlaced lattice L I = f(' j ; s jl ) : ' j = 2 j=P; s jl = d(l + j=2); j = 0; : : :; P ?1; l 2 Z Zg : where P is even and N = P=2.We see that for the standard lattice the detector positions s l do not change with the angle of view.For the interlaced lattice the detector array is shifted by one-half of a detector spacing when going from one angle of view to the next.
In practice one chooses P = 2p for both lattices, and for the interlaced lattice one lets p be even.Then, because of the symmetry relation Rf('; s) = Rf(' + ; ?s); (5) only the angles ' j 2 0; ) need to be measured.It turns out that among all admissible lattices the standard and interlaced lattices are the only ones which fully exploit this symmetry.Indeed, one would need that for 0 ' < , (' + ; s) 2 L implies that also ('; ?s) 2 L. This gives the requirement that P be even and that for 0 j < P=2 and l 2 Z Z there is l 0 2 Z Z with s jl 0 = ?sj+P=2;l , i.e., dl 0 + djN=P = ?dl?d(j + P=2)N=P: This is the case if and only if 2jN P + N 2 2 Z Z; j = 0; : : : ; P=2 ?1: Letting j = 0 gives that N must be even.But then the condition for j = 1 reads 2N=P 2 Z Z which is true only if N = 0 or N = P=2.N = 0 gives the standard lattice, and N = P=2 the interlaced lattice.
We conclude this section with a formal de nition of admissible sampling lattices and the proof of Lemma 3.1: De nition 3.2 An admissible sampling lattice L is a subset of T IR with the following properties: i) L is a subgroup of T IR. ii) L is discrete, i.e., there is > 0 such that if ('; s); (' 0 ; s 0 ) 2 L then either ('; s) = (' 0 ; s 0 ) or max(j' ?' 0 j; js ?s 0 j) .
iii) There is s 0 > 0 such that (0; s 0 ) 2 L. Since L is a subgroup and (0; d) 2 L, we have (0; ld) 2 L for all l 2 Z Z.If (0; s) 2 L with s 6 2 dZZ then there is l 0 2 Z Z such that l 0 d < s < (l 0 + 1)d, and it follows that (0; s ?l 0 d) 2 L. Since 0 < s ?l 0 d < d, this contradicts the minimality of d.
We show rst that M is a nite set.It follows from Claim 2i) that for each ' 2 M there is 0 s d such that ('; s) 2 L. If M has in nitely many elements, then there are in nitely many ('; s) 2 L with 0 ' 2 and 0 s d.Hence there must be an accumulation point which is impossible since L is discrete.
To prove the lemma observe that if M = f0g the Lemma holds with P = 1 and N = 0 because of Claim 1.For M 6 = f0g Claim 3 yields that ('; s) 2 L implies that ' = 2 j=P for some 0 j < P, P > 1.Let s 1 be the smallest non-negative s such that (2 =P; s) 2 L. By Claim 2i) we have 0 s 1 < d.Adding P copies of (2 =P; s) gives (0; Ps 1 ) 2 L which by Claim 1 implies Ps 1 =d 2 Z Z.It follows that s 1 has the form s 1 = dN=P for some integer N with 0 N P ?1. Adding j copies of (2 =P; s 1 ) yields that (2 j=P; jNd=P) 2 L, and then Claim 2 implies that (2 j=P; s) 2 L if and only if s = jNd=P + ld for some l 2 Z Z. 2 4 The ltered backprojection algorithm The ltered backprojection algorithm is the most popular reconstruction method.It is a computer implementation of the approximate inversion formula (3).We derive the algorithm assuming the data are values Rf(' j ; s jl ) of the Radon transform of f, sampled on an admissible sampling lattice L(d;N;P) as described in (4).Discretizing (3) with the trapezoidal rule gives e f(x) ' 2 P P?1 X j=0 Q j (< x; j >); Q j (t) = d X l k(t ?s jl )Rf(' j ; s jl ); with ' j ; s jl as in (4), and j = (cos ' j ; sin ' j ).The reconstruction is usually computed for values of x on a rectangular grid x m 1 m 2 = (m 1 =M 1 ; m 2 =M 2 ).Since computing the discrete convolution Q j (< x; j >) for each occurring value of < x; j > would take too long, one rst computes Q j (iH), jij 1=H, and then obtains an approximation I H Q j (< x; j >) for Q j (< x; j >) by linear interpolation with stepsize H.We assume that H = d=(N 0 m); with 0 < m; N 0 2 Z Z; and N 0 N=P 2 Z Z: (6) This gives H = d=m for the standard lattice (N 0 = 1) and H = d=(2m) for the interlaced lattice (N 0 = 2).Then the e ect of interpolating the convolution is the same as replacing the kernel k with the piecewise linear function I H k which interpolates k at the points Hl, l 2 Z Z; see, e.g., 5, p.84].Hence the algorithm computes the function f R (x) = 2 P P?1 X j=0 I H Q j (< x; j >) = 2 d P P?1 X j=0 X l2ZZ I H k (< x; j > ?s jl ) Rf (' j ; s jl ) As will be seen in the next section, the interpolation step can introduce signi cant errors in certain cases.It was recently shown 21] that the interpolation can be avoided by chosing the points x where the reconstruction is computed on a polar grid rather than on a rectangular grid, and interchanging the order of the two summations.
For references on the ltered backprojection algorithm see 17].Numerous other reconstruction algorithms have been developed; see 19, Ch.V] for a survey.The interlaced sampling lattice, which allows for larger values of d, was rst suggested in 2] based on a geometrical argument.It was rediscovered in 24] by using Shannon sampling theory.More general parallel-beam sampling schemes were found in 4] and applied in 3].E cient sampling schemes for the fan-beam geometry were derived in 20].

Sampling conditions and error estimates
In this section we state the sampling conditions which determine the theoretical resolution of the standard and interlaced lattices, and derive error estimates for the parallel-beam ltered backprojection algorithm.In order to do this we need some Fourier analysis for T IR .The Fourier transform on T IR is de ned by ĝ(k; ) = ( 2  With these preparations we are now ready to analyze the resolution of sampling lattices on T IR.Clearly, the quality of the reconstruction will depend to some degree on the amount of measured data, i.e., the density of points in the sampling lattice.In the subsequent analysis this requirement is expressed as a sparsity condition on the reciprocal lattice, namely that for some compact set M the translates M + , 2 L ? have to be disjoint.This set M plays the role of an essential support of the Fourier transform of the data, in the sense that j d Rf( )j is su ciently small for 6 2 M. A suitable set of this kind was given by Natterer 19] based on results by Lindgren and Rattey 24]: Theorem 5.The parameter b plays the role of a cut-o frequency.If j f( )j is su ciently small for j j > b then the right-hand side of (9) will be small.In this case Theorem 5.3 imposes the condition that the sampling lattice L be such that the translated sets M 0 (#; b)+ , 2 L ? are disjoint.Below we will express this condition in terms of the lattice parameters d; N and P.But rst we will state the error estimate for the ltered backprojection algorithm.Theorem 5.4 Let f 2 C 1 0 ( ), g = Rf, f R be as in ( 7)with e and k as in Theorem 2.1 and in addition such that k 2 L 1 (IR) and P l2ZZ jk(Hl)j 2 < 1.Let M Z Z IR be compact and let g be sampled on an admissible sampling lattice L with parameters d; N; P such that the translates M + , 2 L ? are disjoint.Then, for b > 0, Again the parameter b can be viewed as a cut-o frequency: G H has bandwidth b and so G H e f is a bandlimited approximation to f.Before proving the theorem we discuss its application to the standard and interlaced lattices.This requires to specify a set M, to nd so-called sampling conditions for the lattice parameters d; N; P which ensure that the translated sets M + , 2 L ? are disjoint, and to nd estimates for R Z Z IRnM jF( ; x)jd and for R (Z Z IR)nM jĝ( )j d .We will rst consider the standard lattice.
The last error E 4 can be kept small by chosing the convolution kernel k such that R j j>b j k( )jd is small.In global tomography one can use bandlimited k so that E 4 vanishes entirely.In local tomography this is not possible since k needs to have small support.This concludes our discussion of the errors for the standard lattice.In summary, we expect a good reconstruction as long as the sampling conditions (10) are met.The interpolation stepsize H may be as large as d.
For the interlaced lattice we let M = M 0 (#; b) as in (8).For this lattice P = 2p and N = P=2 = p.We always let p be even, so that because of the symmetry relation (5) only the angles ' j 2 0; ) need to be measured.The reciprocal lattice is L ? = f(p(2k 1 ?k 2 ); 2 k 2 =d); k 1 ; k 2 2 Z Zg.It turns out 19, 5] that the sets M 0 (#; b)+ , 2 L ? are disjoint if either the conditions (10) Comparison of ( 10) and (15) shows that the interlaced lattice allows to double the detector spacing d with only a small increase of p.It is therefore up to almost twice as e cient as the standard lattice.An estimate for R (Z Z IR)nM 0 jĝ( )j d which controls the error term E 2 has already been given in (11).Observe that for j j #jlj the relation ( 12) again implies the estimate (13).Therefore, with b 0 = (1 ?#)b=#, Z ZZ IRnM 0 (#;b) jF( ; x)jd (max j j b j k( )j) X jlj>b 0 Z #jlj ?#jlj jJ l ( jxj)jd (max j j b j k( )j) 4# X l>b 0 le ?l (max j j b j k( )j) 4#e ?b 0 1 + b 0 1 ?e ?+ e ?(1 ?e ? ) 2 !b 0 = (1 ?#)b=#; = (1 ?# 2 jxj 2 ) 3=2 =3: (16) As in the case of the standard lattice, the error E 1 decreases exponentially with b, although at a slower rate.Also, for #jxj close to 1 we expect this error to be signi cantly larger than for the standard lattice, due to the term involving (1 ?e ? ) 2 .Of greater concern, however, is the error E 3 , which unlike as in the case of the standard lattice is now critical.Consider the choice of parameters d = 2 =b, H = =b.Now the sum over l in the estimate for E 3 may have 3 signi cant terms for j j < b: X l2ZZ f(( + 2 l=d) ) = X l2ZZ j f(( + bl) )j ' j f(( ?b) )j + j f( )j + j f(( + b) )j: As discussed before the contribution of the term f( ) is largely cancelled by the factor (1 ?sinc 2 (H =2)).However, this is not the case for the other two terms.E.g., let be close to b.Then, assuming again that f is large near the origin, j f(( ?b) )j will be large and is not attenuated by the factor (1?sinc 2 (H =2)) which will be close to 1. Therefore we expect considerable reconstruction errors for this choice of parameters.That this is indeed the case has been demonstrated for global tomography in 14, 5].Hence when using the interlaced lattice one should choose H b, so that (1 ?sinc 2 (H =2)) is small for j j < b.Typical choices in practice are H = =(16b) or smaller.
As in the case of the standard lattice the error E 4 is kept small be chosing the convolution kernel k appropriately.
In summary, we expect good reconstructions for the interlaced lattice for d close to the optimal value 2 =b if the sampling conditions (15) are satis ed and if H is su ciently small.The error E 1 will be larger than in the standard case but not problematic unless # is close to 1 and the point x where the reconstruction is computed is very close to the boundary of the unit disk .
Proof of Theorem 5.
Using this estimate gives f Combining ( 20) and ( 21) gives the estimate of the theorem for the error E 4 (x).
In order to estimate jf (2)   R (x)j we apply the Poisson summation formula of Theorem 5.1 to the function G of (19).This gives G(' j ; ) = d X where the Poisson summation formula of Theorem 5.2 and the disjointness of the sets M + , 2 L ? was used.Combining the estimates for e 1 and e 2 and observing that sup 2ZZ IR jc q x ( )j p 2 k k(1) k 1 gives the estimates for E 1 and E 2 and nishes the proof. 2

A Numerical Experiment
Detailed numerical experiments for global tomography and simulated data have been reported in 5].Here we demonstrate the possibility of obtaining higher resolution with the interlaced lattice using local tomography and real-world data from the Multidisciplinary Micro CT 3D Imaging Facility at the Mayo Clinic.The data sets consisted of 1440 angles of view distributed over the full range from 0 to 2 .The standard data set had 734 line integrals (rays) per angle of view, while the interlaced data set contained 366 rays per angle.The corresponding reconstructions are shown in the upper left (standard lattice) and upper right (interlaced lattice) of Figure 1.Clearly, the images have approximately the same resolution and are of comparable quality.The image shown in the lower left of Figure 1 shows that using only 366 rays per view with the standard lattice does result in a substantial loss of resolution.As predicted by our error estimates, the interpolation stepsize H had to be chosen very small for the interlaced case.While this increases the computation time, it is hardly a serious concern.More restrictive is the fact that violating the sampling condition on P leads to large errors in interlaced sampling while it has often only very moderate e ects when using standard sampling.The interlaced scheme seems therefore most useful when the detector spacing is the main factor limiting resolution.
+ 2 l=d) j )e i2 ljN=P Splitting the integral over Z Z IR into integrals over M and M c = (ZZ IR)nM gives E(x) = e 1 (x) + e 2 (

Figure 1 :
Figure 1: Top left: Local reconstruction of f with the standard lattice, with 1440 angles of view and 734 rays per view.Top right: Reconstruction with the interlaced lattice, 1440 angles of view and 366 rays per view.Bottom left: Reconstruction using the standard lattice with 1440 views and 366 rays per view.
IR we note that an admissible sampling lattice L has a corresponding \reciprocal lattice" L ? in the Fourier domain.L ? is the set of all 2 Z Z IR such that < y; >2 2 Z Z for all y 2 L. From (4) it follows that L 10, a proof see, e.g.,10, Theorem (8.36)].In order to state the analogous formula for T ?(d; N; P) = f(Pk 1 ?Nk 2 ; 2 k 2 =d) ; k 1 ; k 2 2 Z Zg : Theorem 5.
4: We begin with computing the Fourier transform of the interpolated kernel I H k. Writing I H k as