f ( x,Q e=x c o s e+ y s i n e ) d e y) ( JI0
(3)
where Qe(t), called fitered projections, are related to the p r e jections, Pe(t), by Qe<t)=h(t - a) da
(4)
where the filter impulse response h ( t ) is the inverse Fourier transform of Iwl function in the frequency domain over the bandwidth of the system:h ( t )=
l:
Iw I exp (j2nwt) dw
(5)
II.
RECONSTRUCTION
ALGORITHMS
Let a twodimensional functionf (x, y ) represent across section of the human body. (The property of tissue that is r e p resented by this function will be left unspecified at this time.) A line running through the cross section is called a ray (Fig. 1). The integral of f(x,y) along a ray is called a ray integral and a set of ray integrals forms a projection. A ray integral may be definedmathematically as follows. Theequation of line A B in Fig. 1 is given by xcos8+ysin8=tl (1)
W is the frequency beyond which the spectral energy in any projection may be assumed to be zero. These equations suggest t
he following steps for a digital implementation of the algorithm. Step I (Filtering): Let us say that each projection is sampled with a sampling interval of T cm. In order that the sampled projections do not suffer from aliasing distortion[211,'We have used the symbol w to represent the frequency. If tis mea-
sured in centimeten the dimendomof ware cyclea/cm.
Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.
KAK: COMPUTERIZED TOMOGRAPHY
1247
\.
\
mL
:\
\9’
Fig. 1 . This figure illustrates the variables used i the discussion on the filtered-backprojection algon rithm for parallel projection data. The function Pe(r) is a parallel projection of the function f ( x, y ) .
[42], this implies that W= 1/27. Substituting this value of W in (51, we get for the impulse response
tained from (4):Qe(n7)=7
x P e ( m? ) h ( ( n - m)~).m
(8)
Now since the data is measured with a sampling interval of 7, and is correspondingly band limited, for digital processing the impulseresponseneedonly be known with the samebandwidth and hence the same samplinginterval. We get from (6)
[S.
n=O
n is oddwhere n takes both negativeand positive integer values.At the samplingpoints nr thefilteredprojections maybe ob-
For band-limited functions satisfying W= 1/27 the summation in (8) results in an exact evaluationof Q e ( t ) at t= nr. The discrete convolution in (8) may be implemented directly ona generalpurposecomputer. However, it is much faster to implement it in the frequency domain using fast Fourier transform (FFT) algorithms. (By using specially designed hardware direct implementations of (8) can be made as fast or faster than the frequency-domainimplementation.)For the frequency-domain implementation one hasto keep in mind the fact that one can now only perform periodic (or circular) convolutions[ 5 7] . The convolution required in (8) is aperiodic. To minimize the interperiod interference artifacts inherent to periodic convolution we pad both the projection data and the impulse response function with a sufficient number of zeros. For example,let the number rays in each projection Mrays, of be and the numberof data elements inthe impulse response Mirnp
Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.
PROCEEJXNGS OF THE IEEE, VOL. 67, NO. 9, SEPTEMBER 1979
\(a)
D i
01
(b) Fig. 2. (a) This figure illustrates the variablesused in the discussion on the fdtered-backprojection algorithm for fan-beam projections taken with an equiangular set of rays. S is the source of a fanbeam of rays. (b) T i faure corresponds to (a) for the case when the fan-beam projections are hs taken with a set of rays that equispaced on a straight line such as D, D shown here. are,
Now let (MI2be the smallest integer that is a power of 2 and that is greater than Mra
ys+ Mimp. We zero-pad both the projection data and the impulse response so that each is[MI? elements long (assuming that we want to use a base-2 FFT algorithm). Therefore, the frequency-domain implementationmay be expressed as
Qe(n 7 )= 7 X IFFT{FFT (n7 ) with ZP}{Po
x F F T{~ (~ T ) ZP}} with
(9)
where IFFT denotes the inverse of FFT and ZP stands for zero padding. One usually obtainssuperiorreconstructions when some smoothing is also incorporated in (8) or (9). For example, in (9) smoothing may be implemented by multiplying the product the two FFT’s by a Hamming window. of S t e p 2 (Backprojection): The second step deals with reconstructing the image from the filtered projections using a digital approximation to the integral in (3). When the number of projections Mproj is large enough and uniformly distributed
Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.