KAK: COMPUTERIZED TOMOGRAPHY
1249
over 180°, equation(3) may be approximated by?(x, y )= MprojA
?T
Mproj
Qei(x cos Bii=1
+y
sin B i ) .
(10)
The function f (x, y ) is the reconstructed approximation to the original function f ( x, y ) . I the event one has only a n limitednumber of projections there maybe better approximations[21]. Clearly, thecontribution of the filtered projection at 8, to the reconstruction at (x, y ) is Qei(x cos 8i+ y sin B i ) . The calculation of these contributions for all pixels from one Qei is teed backprojection. The sum of all the backprojections is f (x, y ) . Inbackprojecting a Qei(t) to a point (x, y ), we need to know it for t= x cosdi+ y sin B i . This value of t may not correspond to one of the discrete values for which Qei(t) is known. Clearly, interpolation is necessary. Often, linear interpolation is adequate. In order to eliminate the computations required for interpolation, preinterpolation of thefunctions Q e ( t ) isalso used. In this technique prior to backprojection the function Qe(t) can be preinterpolated onto 10 to 100 times the number of points at which it was originally From this dense set of points one simply retains the nearest neighbor to obtain the value of Qei at x cos Bi+ y sin B i . With preinterpolationand with appropriate programming, backprojection for parallel data can be implemented with virtually no multiplications. We would like to draw the reader’s attention to the article by Brooks and Weiss[ 161whohave shown that greater accuracy in a reconstructed imageis obtained by using linear interpolationthan if the exact filtered projection value was available at every x cos B i+ y sin Bi. Also, Rowland 11041 has shown that linear interpolation has the optimal noise sensitivity among the Lagrange interpolating functions.
A reconstruction algorithm forthe case when projection data is measured with an equiangular set of rays hasbeen rigorously derived by mathematicians Herman and Naparstek[ 661.
A less mathematically rigorous derivation of the same algorithm was recently given by Scudder[ 1061. (The reader is cautioned that there are a couple of misprints that can be misleading in Scudder’s derivation. The reader is also referred to[ 811 .) To present the formulas on which the digitalimplementation is based, we first decide to represent, forthe sake of convenience, the image in polar coordinates. Therefore an image will now be denoted by f ( r,$). (The reconstruction will still be done on a rectangular array.) The projection data will now be denoted by Rp(7) where the angle 7 gives the angular location of a ray in the projection taken at angle 0 (Fig. 2(a)). To facilitate ourpresentation wewill also need to define two new parameters L and 7‘.In a given projection at angle 0, L is the distance from the source S to the pixel at(I,$):
L(0, r, 4)= dD2+ r2+ 2Dr sin (0 -$).
(1 1)
The variable 7’will give us the angular location of the ray that passes through a given pixel ( r,$) in a given projection at 0:
The image f ( r,$) and the fan-beam projections Rp(7) can be shown to be related by
‘(”;’ - 7)2 sin (7
h ( r’ - T ) d y d@ (13)
where -7m and rm are angles for the extreme rays in each B. Reconstruction Algorithm for Fan-Beam Projections projection, and where the function h ( 7 ) is the same as in (6) Generated by Equiangular Rays (with argument t replaced by 7). The relationship in (1 3) suggests a weightedfiltered-backAlmost all fast CT sanners today do reconstructions from the implementation. This wbe l l i fan-beam projections. In this and the next two subsections we projection algorithm for will only discuss algorithms that directly reconstruct an image demonstrated by the following steps. Step 1 (Modify Each Projection): Let us now assume that from fan-beam projections. Note that it is also possible to reeach fan-beam projection Rp(7) is sampled with an angular arrange the fan projection data into parallel projections and then use the algorithm described before[47],[49],[ 971, sampling interval of a radians. We will again assume that the[ 1211. T i rearrangement can be done“on the fly” provided sampled data R g ( n a ) is free of aliasing errors. The first step hs the angular interval between fan projections and the sampling consists of obtaining for each projection R p ( n a ) a corresponding R$ n a ) as follows interval in each projectionsatisfy certain conditions. There are two types of fan-beam projections: those that take R b ( n a )= R p ( n a ) D cos n a . (14) the projection data with an equiangular set of rays, and those that utilize a set of rays that result in equispaced detector lo- Note that n= 0 covesponds to the center ray in each projecl call l i cations on a straight line. In this subsection we will present tion. We w R g ( n a )modified projections. Step 2 (Filtering): I this step we filter each modified pron a digital implementation for the former case. Note that when
the rays are equiangular, the detectors for measuring projec- jection by convolving it with an impulse response g ( 7 ) defined tion data areequispaced along the arcDl D2 shown in Fig. 2(a). by The radius of this arc is 2 0 where D is the source to center distance. where h ( 7 ) is given by (6) with t replaced by 7. In the continuous domain t i operation of filtering is described by hs