ProcIEEE_Kak_computerized_tomography_with_xray_emission_ultrasoundsource

PROCEEDINGS OF THE IEEE,

VOL. 67, NO. 9, SEPTEMBER 1979

1245

Computerized Tomography with X-Ray, Emission, and UltrasoundSourcesAVINASH c. KAK, MEMBER, IEEE

A b m t - T h i s paper reviews the major developments that have taken images with this resolution the more recent scanners use from phce during the hst three y e m in imaging with computed tomography 500 to 1000 projections, and from 500 to 1000 rays in each (cr)wiug x-my, emi on, md u d l sourws. space limitations projection. The reconstruction time has also gone down conhave resulted m some selection of topics by the authar. There are four major sections dealing with algorithms, X-ray CT, emission CT, and siderably over the years. By exploiting parallelism inherent in ultrasound CT. S e most of the cmreatiy USXI are of the reconstruction process and array processing the per slice i n algxithms wtered*kprojdon tp, we have coacentnted ( 1 tfrese m the reconstruction times have beenbroughtdown to less than ye 1 section on algorithm (with emphasison their implementation a p cs. 30 s. s et) I n X - m y C T m i m p o r C P n t q~ n i s e d d u r i n g t h e k s t f e w y e~ b r s A great deal morecould be said abouttheperformance concernedthepaune6er~byaCTsanner,gnenthetrctthat

~ X - r r y s~ i n C T s a n w r s r r e~ d u o m a t i c m d t h e f r c t t h . t t i characteristics andother capabilities of scanners on the market. si sue attenuation coefficients are energy dependent. A m e to this However, the a m of this paperis not totabulate such informar yn questionarereviewedmthe~onX-my~wherewehaverlso tion (which after all is readily available from the manufacturdiscussed the a cIp&d by the pdychr~maticityd the X-my ers) but to review the recent advances in the basic science and photons.Methods f a the removal of these artifacts have .Is0 been r ve e . InemissionCTtheb&stdevelopmentofthehtthree engineering of computerized tomography. eiwd so is m alth~& sl#~e m - Although the results far have been particularly spectacular~ yerr~ the W t m t in p i t r o n -hy, have dictated an essentially mlmductay treatment rad not rll with X-ray transmission,thetechnique of computerizedtorspeCg of the h g photoa and positron tomognphy have been sur- mography has also been successfully applied in nuclear medik veyed. we have reviewed recent developments in ultrasound cine imaging. The result of X-ray CT is a two dimensional map CT. W have pointed out that becam of the sensitivity of this tech- of the linear attenuation coefficient values (at some effective e nique to refmcticm, it ie currently limited to soft tissue structures, with ultmsonic detection of tumors m the fern& breast a significant photon energy). WithCT in nuclear medicine, on the other application. hand, the image yields information about a location and concentration of aphotonorpositronemittingradioisotope, I. INTRODUCTION which is administered to the patient either by inhalation or OUNSFIELD'S invention of the computed tomography

by injection. For this reason such scanners in nuclear medi(CT)scanner, revealed tothe public in 1972, was a cine arecalled emission CT (ECT) scanners. major breakthrough in biomedical imaging[69],[70]. The clinical usefulness ofX-ray and emissionCT systems To be sure, even before this invention there were researchers has already been established. X-ray derives CT its clinical who were concerned with the problem of constructing cross- significance from the fact that it shows the anatomical morsectional images of an object from its projections. Theearliest phology of any desired cross section of the human body withmathematical work on this subject is by Radon[ 101I. More out superposition, with great sensitivity (around 1 percent)2 recently, some of the fmt investigators to examine this p r o b and with good resolution. The clinical usefulness of emission lemtheoretically and/or experimentally (and often indepenCT lies in the fact that after a gamma-ray emitting isotope is dently) include (in roughly a chronological order): Bracewell administered to a patient, its distributionin the body gives[lll,Oldendorf[961,Cormack[401,[411,KuhlandEdwards physicians information regarding the functional state of[841, DeRosier and Klug[451, Tretiak e t u1.[1161, Rowley various organs. Andby examining the temporal dependence[ 1051, Berry and Gibbs[ 91, Ramachandran and Lakshmina- of thisdistribution,dynamicstudies may be performed on rayanan[102], Bender et al.[8], Hermanand Rowland[641, various organs. and Bates and Peters[ 71. Just on the horizon there is yet another tomographic imagDuring the last few years great improvements have been ing mode currently undergoing clinical t i l . This is the techras madein X-ray tomography, with the result that the patient nique of ultrasonic computed tomography (UCT). The basic scan time has been reduced to less than 10 s from over 4 min. On some of the latest scanners thehigh contrast spatialresolujs recon'Thisdoesnotimplythatlinearattenuationcoefficient with accuracy. Because of hardening beam and tion is beginning to approach10line pairs per centimeter structed 1-percent (discussed here in Section 111-B) the of andthelowcontrast(1-percentcontrast) spatialresolution other effects attenuation coefficient (at the effective numerical valuesthe the computed energyof is now around 3 to 5 line pairs per centimeter, the resolution CT scanner) may be distorted by as much as 5 percent or more[83]. criterion in both cases being visual detectability.' To present Insome cases thiserror might be much greater. In adramaticex-

strain*

F, e

w

Manuscript received December 18, 1978;revised March 15,1979. The author is with the School of Electrical Engineering, Purdue University, West Lafayette, IN 47907. 'The low contrast resolution dependa upon the signal-to-noise ratio in the reconstructed image and, therefore, is a function of radiation dosage.

I . ample DiChiro et Q(461 comparedthetomogramsofaheadfor a slice taken

near the apex (where the beam hardening effects are the greatest) with a slice takenatalower level. Theyshowedthatthe attenuation value computed for cerebrospinal fluid was 30 (CT units) in the former case and 11 in the latter. Apparently, the severe beam hardening caused this large change in the computed CT number. The often quoted 1-percent sensitivity onlyrefers to the ability of a CT scanner to distinguish between tissues(over a few square millimeters area) whose attenuation coefficients may differ by only 1 percent.

0018-9219/79/0900-1245$00.75 1979 IEEE 0

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1246

PROCEEDINGS OF THE IEEE. VOL. 67, NO. 9, SEPTEMBER 1979

goal here is the same as with X-rays, namely, to obtain quantitative crosssectional images depicting morphological detail of the humanbody. However, incontrastwith X-rays,CT imagingwith ultrasound is madedifficultby thefactthat rays of sound'.energy do not necessarily travel along straight lines intissue,andmayundergoseriousrefraction at interfaces between soft and hard tissues. On account of this,a p plications of ultrasound CT appear t o be limited for the foreseeable future to tho& parts of human body that are free of hard tissues such as bone.Thusmost of thecurrenteffort in t i area is focused on the clinically important problem of hs tumor detection in the female breast. Common to all threedifferent CTimaging modalities described above are the computer algorithms for reconstructing tomograms projection from data. Al reconstruction algol rithms are based on the assumption that the projections are line integrals of the quantity of interest. The extent to which ti assumption is violated determines the magnitude of artihs facts visible in the tomogram. Examples of such artifacts include those caused by beam hardening in the case of X-ray C T, those caused byinadequate or inaccurateattenuation compensation for emission CT, and those caused by refractive effects in ultrasound CT. Study of artifacts and ways to minimize them has assumed importance during the last threeyears, and a number articles hasbeen published on this subject. of Our aim in this paper isto do a review of the major develop mentsthat have taken place in CT imaging during the last three years. We have focused on this period for two reasons: i) the relevant science and engineering have grown enormously during this period; and ii) review papers by other authorshave morethanadequately covered thetime period before1976 (excepting ultrasound CT which has not been reviewed so far). The review article by Brooks and DiChiro[ 141 surveys X-ray CT for work published till 1976. Emission CT hasbeenreviewedbyTer-Pogossian[ 1131, Phelps[ 991, and Budinger[271. Because of space limitations there had to be some selection in the topicspresentedhere. In emission CT the biggest development of the last three years is the great

interest that has developed inpositrontomography. While have we reviewed most of the recent advances in single photon emission CT, our treatment of positron emission CT is essentially introductory. Positron emission CT presents aseries of unique technical problemsthat are not addressed here. In what follows, Section I1 deals with recent developments in reconstruction algorithms. Section 111 is on X-ray CT and examines the degree to which the projection data generated on CT scanners conform to the assumptions made in the development of algorithms. Section IV is on emission CT and includes positron both tomography tomography and using single gamma-ray counting. Section V is on ultrasound CT.

where tl is the perpendicular distance of the linefrom the origin. The integral of the function f(x, ) along this line may y be expressed as

1=

ray AB

I

f(x,y)ds

where 6(*) the Dirac delta function. The function P e ( t ) as is a function of t (for a given value of 0) defines theparallel projection of f (~y, ) for angle 8. The two dimensional function P e ( t ) is also called the Radon transform of f(x, y ) . One may also generate projections by integrating a function along a set of lines emanating from a point source as shown in Fig. 2. Such projections are called fan-beam projections. A reconstructionalgorithm tells us how to reconstructa function f(x, y ) from its projections. Over the past few years many such algorithms have been developed[ 141. The algorithms that are currently being used on most if not all CT scanners are of the filtered-backprojection type. We will now present and discuss the steps involved in the implementation of these algorithms. We will assume that the projections are uniformly sampled, which is the case most often encountered in practice. For the case of nonuniform sampling the reader is referredto[681.A . Filtered-Backprojection Algorithm for Parallel Projection Data T i algorithm is based on following relationships[ 121,1791, hs 11021,[lo71

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.

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

‘Sometimes,because of its speed, the FFT interpolation method is used (see, for example, reference[ 5 7, p. 1901). Thismethodconsistsof takingthe FFT of aprojection,paddingthetransformwith a large number of zeros and taking the IFFT. It has recently been shown[ 82 1 that this method can lead to large errors especially near the beginning and the end of a data sequence.

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1250

PROCEEDINGS OF THE IEEE,

VOL. 67, NO. 9, SEPTEMBER 1979

where Qs(7)is the filtered projection correspondingto Rb(y). we need to introduce two more parameters U and t'. ParamIts digital implementation is given by (exactforthe band- eter U is for each pixel ( x,y ) (or equivalently ( r, 9)) the ratio of (Fig. 2(b))tothesourceto origin distance.Note that limited case and when a satisfies the Nyquist condition): E is theprojection of thesource to pixel distanceonthe Qg(na)= a~; ( m g((n - m) a). a) (1 7) central ray.m

The discrete impulse response g ( n a ) is given by the samples of (1 5: )n=O

-D+rsin(o-$) D

(21)

t' The other parameter is the value of t that correspondsto the iay passing through the pixel (r,$). In Fig. 2(b) t' is equal to the distance OF The convolution in (1 7) possesses a frequency-domain implementation very similar to (9). The comment we made there With all the parameters defied as above, theimage f ( r, 6) and regarding the desirability of incorporatingsomesmoothing its projectionsR g ( t ) can be shown to be related by also applies here. Step 3 (Fan Backprojection): This last step consists ofimplementingthe remaining integration (13). in In terms of filtered projectionsti step may be expressed as hs (23) where the function h ( t ) is the same as in ( 6 ) . As was the case with (13),equation(23) possesses a weighted filtered-backWhen the number of projections Mproj is large and uniformly projection implementation. Except different for weighting distributed over 360°, for digital implementation (19) may be factors and fiter function, 'the steps are almost identical to expressed as those in the preceding case and will not be repeated here. .

111. X-RAYTRANSMISSION COMPUTED TOMOGRAPHYwhere x= r cos$ and y= r sin$. The function?(x, y ) is the reconstructedapproximation tothe image f ( r,$). To find the contribution of the filtered projection Qsi to the reconstruction at a pixel ( x,y ) we must first find L and y' using (1 1) and (1 (This is in sharp contrast to the case of parallel 2). projections wh

ere in order to compute such a contribution we only had to find x cos Bi+ y sin B i . ) Computation of such contributions at all pixels from one Qai is called a fan backprojection. The sum o,f all such fanbackprojectionsmultiplied by 2n/Mproj gives f ( x,y ) .C. Reconsttuction Algorithm forFan-Beam Projections Taken with Equispaced Detectors on a Straight Line

Consider a hypothetical parallel beam of monoenergetic X-rays propagatingthrough tissue as shownin Fig. 3. Let Ni, be thetotalnumber of photonsthatentertheobject (withinthemeasurementtime interval) through the beam from side A; andlet Nd be thetotalnumber of photons exiting (within the same time interval) through the beamon side B. In the plane shown in the figure the tissue everywhere may be characterized by its linear attenuation coefficient. Let p ( x, y ) be the attenuation coefficient at location ( x, y ) in the plane shown. When the width T of the beamis sufficiently small, the numbers N d and Ni, are related by[ 621,[ 1121:=Nin Nd

exP[ - J P ( x, Y ) d~]

(24)

or, equivalently,Nd Jayp(x,y)ds=1n- Nin

A reconstructionalgorithmfor t i casewas f i t derived hs by Lakshminarayanan[ 851 (see also[ 651 ). It can also be more simply derived[ 81 I by using a method similar to that of Scudder[ 1061 for the equiangular case. Fig. 2(b) shows the geometry for data collection. Although theprojections are measuredon a linesuch as DlD2, for theoreticalpurposes it is moreefficient ty assume the existence of an imaginary detector line DlDz passing through the origin. We now associate a ray integral along SB with point F on D; Dh as opposed to point B on Dl Dz. Each projection is now denoted by R&) where the distance t along the imaginary detector line D; Dh identifies the location of a ray in the fan projection at angle fl. Again, let D be the source to center distance; and f ( r,$) the image we seek to reconstruct in polar coordinates. For ease in presentation

(25)

where ds is an element of length and where the integration is carried out along line A B shown in the figure. The left-hand side precisely constitutes a ray integral for a projection. Now an ideal detector placed on side B will produce a signal directly proportional to N d . Therefore, measurementslike In (Ni,/Nd) taken for different raysat different angles may be used to generate projection data for the function p ( x, y ) . We would like to remind the reader that this is strictly true only under the assumptionthatthe X-ray beam consists of monoenergetic photons. This assumption is necessary because the linear attenuation coefficient is in general a function of photon energy.

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

KAK: COMPUTERIZEDTOMOGRAPHY

1251

P8(k7)=[

p(s, y)dr= an N i nray path AB

NdX-ray source

B\

\

Detector Array Rotating instrument Plate

Fig. 4 . The figure depicts a fan-beam rota

tional type scanner.

p\\

Ninwidth o f one detectorT

A

Fig. 3. A parallel beam o f X-rays propagating through tissue.

Later in this section we will address thefactthat since, in practice, this assumption is never completely satisfied, artifacts are created in X-ray CT images. In this section, we wdiscuss three topics concerningX-ray l i CT: two scan configurations that have been introduced during the last three years, polychromaticity artifacts and the properties of noise in the reconstructed images. We have not discussed X-rayCTusingvideo recorded fluoroscopic data[ 31 -[ 51,[ 791. Ideally, a CTimage should represent a cross-sectional slab of constant thickness. We have not discussed here the variations in this thickness caused by the collimation geometry[191.A . Different Scan Configurations in X-Ray CAT There are two scan configurations that lead t o rapid data collection. These are i) fan-beam rotational type; and ii)fixed detector ring with rotating source type. 1) Fan-Beam Rotational Scanners: Here a fan-beam of Xas shown in rays is used to illuminate a multidetector array Fig. 4. Both the source and the detector array are mounted on ayoke which rotatescontinuouslyaround thepatient. Data collection time for such scanners range from 1 to 20 s. In t i time up to a 1000 projections maybe taken. If the hs projections are taken 'on the fly' there is a rotational smearing present in the data; however it is usually so small that its effects are not noticeable in the final image. Most such scanners use fan-beams with fan angles ranging from 3' to 45'. 0 The detector bank usuallyhas 500 to 700 detectors. Images are reconstructed on 256 X 256 or 5 12 X 5 12 matrices. Most such scannersuse xenon gas ionization detectors. Three such detectors are shown in Fig. 5. Each detector consists of a central collecting electrode with a high-voltage strip on each side. X-ray photons that enter a detector chamber cause ion-

L_ L

Ionization chrnber forxenonone detector

x-ray

photons

A~Thigh wltage surface alumlnm entrance windm a t grovld potential hlgh voltage

'(-vel

electrode collecting

electrode

Fig. 5. Xenon gas detectors for fan-beam rotational type Scam'wrs.

izations with high probability (which depends upon the length I of the detector and the pressure of the gas). The resulting current through the electrodes is a measureof the incident X-ray intensity. In the EM1 6000 scanner the collector plates are made of copper and the highvoltage strips of tantalum. In the same scanner the length I (shown in Fig. 5) is 8 cm, the voltage applied between the electrodes 170 V and the pressure of the gas 10 atm. The overall efficiency of this particular detector is around 60percent. The primary advantages of xenon gas detectors are that they can be packed ciosely and they are inexpensive. The entrance width T in Fig. 5 may be as small as 1 mm. This is in direct contrast t o the scintillator-photomultiplier detector (to be described later) which

suffers from the disadvantage that the smallest commercially available photomultiplier tube has a diameter of 12 mm. Yaffe et al.[ 1251 have discussed in detail the energy absorption efficiency,linearity of response and the sensitivity t o scattered and off-focus radiation for xenon gas detectors. Williams[ 1231 has discussed their u e in commercial CT syss

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1252

PROCEEDINGS OF THE IEEE, VOL. 67, NO. 9, SEPTEMBER 1979

(c) Fig. 6. Three examples of reconstructions obtained on a fan-beam rotational scanner. Number of projections: 504. Number of rays in each projection: 504. (Courtesy of EMI.) (a) Noncontrast CTscan at thelevel of thethalamus. Theanterior horns of the lateral ventricles, superior cerebellar cistern and slightly calcified pineal gland are seen. There is differentiation -of gray and white matter in the region of the head of thecaudate nucleus and the internal capsule. (b) High skull CTscandemonstrating corticalatrophy. There is differentiation of gray and white matter in the individual gyri. The falx cerebri i also seen. (c) CT scan of the abdomen revealing a right renal cyst. Vascular s structures demonstrated include the aorta, celiac axis, hepatic and splenic arteries, inferior vena cava and portal vein.

tems. In Fig. 6 we have shown three examples of reconstructions made on a fan-beam rotational scanner. 2 ) Fixed Detectors Rotating Source Scanners: In these scanners a large number of detectors are mounted on a fixed ring as shown in Fig. 7. Inside this ring is an X-ray tube that continuallyrotatesaroundthepatient. During this rotation the output of the detector integrators facing the tube is sampledevery few milliseconds. A l such samples foranyone l detector constitute what is known as a detector fan. Of course, as far as a reconstruction algorithm is concerned there is no difference between a detector fan and the more conventional source fan. Note that the detectors do not have to be packed closely (more on this toward the end of this section). This observation together with the fact that the detectors are spread all around on a ring dictates the use of scintillation detectors as opposed to ionization gas chambers. Most scintillation detectors currently in use are made of sodium iodide, bismuth germanate and cesium iodide crystals. (See[43] for a comparison of sodiumiodideandbismuth germanate.) Thecrystal of which a scintillation detector is madeserves two purposes. First it traps most of the X-ray photons which strike it, with a degree of efficiency which depends upon the photon energy and

the size of the crystal. The X-ray photons then undergo photoelectric absorption Compton (or scatter subsequent with photoelectric absorption) resulting in the production of secondary electrons. The second function of the crystal is that of a phosphor-a solid which can transform the kine

tic energy of the secondary electrons into flashes of light. The geometrical design and the encapsulation of the crystal are such that most of these flashes of light leave the crystal through a side where they can be detected by a photomultiplier tube. Such scanners were fust developed by A.S.& E. Corporation(whicharenow sold by Pfizer Corporation). An example of a scanner of this type is also the EM1 7000 system which uses 1088 cesium iodide detectors and in each detector fan 1356 samples are taken. This system differs from the one depicted in Fig. 7 in one respect, namely the X-ray source rotates around the patient outside the detector ring. This makes it necessary to nutate the detector ring so that measurements figure may be made 1631. Fig. 8 likethoseshowninthe shows two reconstructions made such a scanner. on 3) Discussion o f Scan ConfigurationDifferences: An i m portant difference exists between the two scan configurations discussed here.The data in the fan-beam rotatingdetector scanners is essentially limited in the number of rays each pro-

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

KAK: COMPUTERIZED TOMOGRAPHYa f i x e dr i n g ofdetectors

1253

anx-raysource r o t a t i n ga r o u n d t h ep a t i e n t

Fig. 7. A fixed detector-ring rotating source type scanner.

Fig. 8 . Two examples of reconstructions obtained on a fixed detector-ring rotating source type scanner. Note the superior resolution and contrast. (Courtesy of EMI.) (a) Noncontrast CT scan of the mid-brain demonstrating the third ventricle, frontal horns of the lateral ventricles, and quadrageminal cistern. Note the folia of the cerebellar vermis. (b) CT scan at the lumbarlevel. The psoas muscles are well demonstrated al as are the aorta and vena cava. Leaves of the mesentery are dearly seen. The w l of the colon is outlined by air and feces. Note the lumbar theca and anterior paraspinal veins. No contrast was used on the scan.

jection can have, although there is no limit on the number of projections. One can have only as many rays in each projection as the number of detectors in the detector array. On the

other hand, the data collected in the fixed detector scannersis limited in the number of projections that may be generated, while there is no limit on the number of rays in each projec-

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1254

PROCEEDINGS OF THE IEEE,1

VOL. 67, NO. 9, SEPTEMBER 1979I

tion.’ (It is now known that for good quality reconstructions the number of projections should be comparable to the number of rays in each projection.) In a fan-beam rotating detector scanner if one detector is defective, the same ray in every projection gets recorded incorrectly. Such correlatederrorsin all theprojectionsform ring artifacts[ 1081. On the othe

r hand when one detector fails in the fixed detector ring type scanners, it implies a loss or partial recording of one projection. When a large number of projections are measured, a loss of one projection usually does not noticeably degrade thequality of reconstruction[ 1091. The very nature of construction of the gas ionization detectors in the fan-beam rotational scanners lends them a certain degree of collimation whichis a protection againstreceiving scatter radiation. On theotherhand,thedetectorsinthe fixed ring scanners cannot be collimated since they must be capable of receiving photons from a large number of directions as the X-ray tube is rotating around the patient. This makes Fig. them more vulnerable to scatter radiation 861.[ B. Difficulties with the Interpretation of CTNumbers The reader will recall from the introduction to this section that the ray integral in (25) is given by ln(Ni,/Nd) only under the assumption that photons are monoenergetic. In practice this assumption is violated. Fig. 9 shows an example of an experimentally measured X-ray tube spectrum taken from Epp and Weiss[ 5 11 for an anode voltage of 105 kVp. For polychromatic photons (24)has to be replaced by6 (26) where&,(E) represents the incident photon number density (also called energy spectral density of the incident photons).&,(E) dE is the total number of incident photons in the energy range E and E+ dE. This equation incorporates the fact that the linear attenuation coefficient p at a point ( x, y ) is also a function of energy. The reader may note that if we were to measure the energy spectrum of exiting photons (on side B in Fig. 3) it would be given bySexit(E)=Sin(E)exP

Energy I n K e V

9. An example of an experimentally measuredX-ray trum. (From Epp and Weiss[ 5 1].)

tubespec-

[-I]P(x,y,E)ds

(27)

In the energyrangesused for diagnostic examinations the linear attenuation coefficient for many tissues decreases with energy. For a propagating polychromatic X-ray beam this causes the low energy photons to be preferentially absorbed, so that the remaining beam becomes proportionately richer’Althoughone may takea verylarge number of samplesineach projection,“useful information” would now be limited by the width of thefocalspot on the X-ray tube andby thesize of thedetector aperture. In discussing polychromatic X-ray photons one has to bear i mind n typesof detectors 1881. The thatthere are basicallythreedifferent output of a detector may be proportional to the total number of on it; or it may be proportional to total photon photons incident energy; or it may respond to energydeposition per unit Most countingtypedetectors are of thefirst type, mostscintillation type detectors are of the second type, while most ionization detectors are of the third type. In determining the output of a detector one must also take into account the dependence of detector sensitivity on photon energy. In this paper, we w assume, for the sake of simplicity, that

l i over the energy range of interest the detector sensitivity is constant.

in high energy photons. In other words, the mean energy associated with the exit spectrum Sexit(E), is higher than that associated with the incident spectrum Si,(E). This phenomenon is called beam hardening. Given the fact that X-ray sources in CT scanning are polychromaticand that the attenuation coefficient is energydependent, the following question arises:What parameter does an X-ray CT scanner reconstruct?To answer t i question hs McCullough[ 871,[ 881 has introduced the notion effective of energy o f a CT scanner. It is d e f i e d as that monochromatic energy a t which a given material will exhibit the same attenuation coefficient as is measured by the scanner. McCullough e t al.[ 871 showed empirically that for the original EM1 head scanner the effective energy is 72 keV when the X-ray tube is operated at 120 kVp. (See[ 921 for a practical procedure for determining the effective energy of a CTscanner.) The concept of effective energy is valid only under the condition that the exit spectra are the same for all the rays used in the measurement of projection data. (When the exit specfra are not the same, the result is the appearance o f beam hardening artifacts discussed in thenextsubsection.) It follows from the work byMcCullough[881 that it is a good assumption thatthe measured atknuation coefficient, h e a s u r e d, at a point in a cross section is related to the actual attenuationcoefficient p(E) at that pointby

T i expression applies only when the output of the detectors hsis proportional to the total number of photons incident on them. McCullough has given similar expressions when detectors measure total photon energy and when they respond to total energy depositionlunit mass. Effective energy of a scanner depends not only on the X-ray tube spectrum but a s lo on the nature of photon detection. Although it is customary to say that a CT scanner calculates the linear attenuation coefficient of tissue (at some effective

mass.

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

1255

energy), the numbersactually putoutbythecomputerattached to the scanner are integers that usually range in values from - 1000 to 1000. These integers have been given the name Hounsfield units and are denoted by H . The relationship between the linear attenuation coefficient and the corresponding Hounsfield unit isH=

Cl - Clwater x loooPwater

(29)

where p water is the attenuation coefficient of water, the value of both p and ClWter are taken at the effective energy of the scanner. The value H= 0 corresponds to water; and thevalue H= - 1000 corresponds to p= 0, which is assumed t o be the attenuation coefficient of air. Clearly, if ascanner wasperfectly calibrated it would give a valueof zero for water and - 1000 for air. Under actual operatin

g conditions this is rarely the case. However, if the assumption of linearity between the measured Hounsfield units and the actual view of the attenuation coefficient (at the effective energy of the scanner) is valid, one mayuse the following relationship to convertthe measured numberH, into theideal numberH H=Hm

- Hm, water(30)

x

1000

H m, water- H w, air

where Hm,water and Hm,& are, respectively, the measured Hounsfield units for water and air, respectively. (This relationship may easily be derived by assuming that/I= aHm+ b, cal. culating a and b in terms ofHm, water, Hm, air, and Pwater, and then using (29).) Brooks[ 181 has used (28) to show that theHounsfield unit H a t a point in a image may be expressed as CT

A:

Polychromatic case

II

I

1

i~

'\\

/

/

8: Honochranatlc case where H, and H p are the Compton and photoelectric coefficients of the material being measured, expressed in Hounsfield units. The parameter Q, called spectral factor, depends only upon the X-ray spectrum used and may be obtained by performinga scan on acalibratingmaterial. A noteworthy feature of H, and H p is that they are both energy independent. Fig. 10. (a)Reconstructionfrompolychromaticprojection data ofa waterphantominsidea skull. Notethewhitening effect near the Equation (31) leads to the important result that if two difA quantitativeillustration of the effect in(a). Curve A ferent C images are reconstructed using two different incident skull. (b)thereconstructedvalues T represents on the middle horizontal line spectra (resulting in two different values of Q ), from the re(a): Curve E represents reconstructed the values on the same lme w t h monochromatic X-rays. sulting two measured Hounsfield units for a given point in the cross section, one may obtain H, and H p leading to a degree ofchemical identification of the material at that point. Inwas done stead of performing two different scans, one may also per- Phelps e t al.[981.Reconstructionfromthisdata using the filtered backprojection algorithm (Section 11-A) form only one scan with split detectors for this purpose[20]. with 10 1 projections and 101 parallel rays in each projection. C. Polychromaticity Artifacts in X-Ray CT Notethe"whitening"effect near the skull in Fig. lO(a). This is more quantitatively illustrated in Fig. 10(b) where the Beam hardening artifacts, whose cause was given above, are most noticeable in the CT images of the head, andinvolve two elevationof the reconstructed values near the skullboneis different types of distortions. Many investigators[ 151,[461, quite evident. When CT imaging was in itsinfancy this whiten[ 541,[go] have shown that beam hardening causes an eleva- ing effect used t o be mistaken for gray matter of the cerebral wehave also shown in Fig. 10(b) tionin CT numbersfor tissues close to the skull bone. To cortex.Forcomparison, illustrate t i artifact we have presented in Fig. 10 computer the reconstructed values obtained when the projectiond

ata hs simulation reconstructions of a water phantom inside a skull. was generated for monochromatic X-rays. The other artifact caused by polychromaticity is theappearThe projection data was generated on the computerusing Epp and Weiss[ 5 1] 105 keVp X-ray tube spectrum(Fig. 9). The en- ance of streaks and flares in the vicinity of thick bones and ergy dependence of the attenuation coefficients of the skull between bones[SO],[781,[831. (Note that streaks can also bone was taken from theICRU report[ 731 and of water from be caused by aliasing 121 I,[421.) This artifact is illustrated in

i

e

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1256

PROCEEDINGS OF THE IEEE, VOL. 67, NO. 9, SEPTEMBER 1979Ideal case

(no beam hardening)

case

Thickness o f a homDgencous absorber

Fig. 12. The solid curve shows that the experimental measurement o f a ray integral depends nonlinearly on the thickness o f a homogeneous absorber.

cessing of the projection data and 2) postprocessing of the reconstructed image. Preprocessing techniques are based on the following rationale. If the assumption of the photons being monoenergetic was indeed valid, a ray integral would then be given by (25). For a homogeneous absorber of attenuation coefficient p this implies

(b) Fig. 1 1 . (a) Reconstruction frompolychromatic projection data o f a phantom that consists o f a skull with.four circular bones inside. The rest of the“tissue” inside the skull is water. The dark and wide streaks are caused by polychromaticity the of X-rays. (b)Reconstruction of the same phantom as in (a) from projections generated with monochromatic X-rays. The variations in the gray levels outside t o bone areas within the skull are less than 0.1 percent of the mean value. The image was displayed with a narrow window t o bring out these variations. Note the absence of dark streaks present in (a).

Fig. 11. The phantom used is a skull with water and four circular bones inside. Polychromaticprojectiondata was generated as before using the105 keVp X-ray spectrum. The reconstruction using this data is shown in Fig. 1l(a) with the same number of rays andprojections as before. Note the wide dark streaks between the bones inside the skull. Compare this image with the reconstruction shown in Fig. 1l(b) for the case when the X-rays are monochromatic. In X-ray CT of the head similar dark and wide streaks appearin those cross sections that include the petrous bones, andaresometimes called the interpetrouslucency artifact. Various schemes have been suggested for making these artifacts less apparent. These fall into two categories: 1) prepro-

where I is the thickness of the absorber. This equation says that under ideal conditionsthe experimental measurement ln(Nin/Na) should belinearly proportional to the absorber thickness. This is depicted in Fig. 12. However, under actual conditions a result like

the solid curve in the figure is obtained. Most preprocessing corrections simply specify the selection of an“appropriate” absorber and then experimentally obtain the solid curve in Fig. 12. Thus should a ray integral be measured to be at A, it is simply increased to A’ for tomographic reconstruction. This procedure has the advantage of very rapid implementation and works well forsoft tissue cross sections because differences in the composition of various soft tissues are minimal (they are all approximately water like from the standpoint of X-ray attenuation). For preprocessing corrections the reader is also referred to references[ 151,[ 891,[ 901, and for a technique that combines preprocessing with image deconvolution see[ 341. Preprocessing techniques usually fail when bone is present in a cross section. In such cases it is necessary to postprocess the CT image. Methods for doing so may be subdivided into twotypes:i) iterative[781,[83], and i dual energy techi ) niques of Macovski and his associates[ 1],[ SO]. In the iterative scheme one first does a reconstruction (usually incorporating the linearization correction mentioned above) from .the projectiondata. This reconstruction is then thresholded to get an image that shows only the bone areas. This thresholded image is then“forwardprojected” to determine the contributionmade by bone to each ray integral in each projection. Based on this contribution a correction is applied to each ray integral. Joseph and Spital[78] and Kijewski and Bjangard[ 831 have obtained very impressive results with this technique.

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

KAK: COMPUTERIZEDTOMOGRAPHY

1257

The dual energy technique proposed by Alvarez and Macovski[ 11 is theoreticallythe most elegant approach to eliminating the beam hardening artifacts. Their approach is based on modeling the energy dependence of the linear attenuation coefficient byP ( X, Y, E )= a1 ( x, y ) g(E)+ a2(x, Y ) fm(E).

and I2(-41,A2)=[s2(E)eXP[-(Alg(E)+AZfKN(E))] dE(40)

which gives us two (integral) equations for the twounknowns Thepart u1 (x, y)g(E) describes the contribution made by A 1 and A 2 . The two source spectraS1( E ) and S2 ( E ) may for photoelectricabsorption to the attenuation at point (x, y ); example be obtained by simply changing the tube voltage on a l ( x, y ) incorporatesthe material parameters at ( x, y ) and the X-ray source. This, however, does require that two scans g(E) expresses the (material independent) energy dependence be made for each tomogram. In principle, one can obtain equivalent results from a single scan with split detectors[ 201. of this contribution. The function g(E)is given by Alvarez and Macovski[2] have shown that statistical fluctua1 g(E)= 3. ( 3 4 ) tions in a1 ( x, y ) and a 2 ( x,y ) caused by the measurement erE rors in I and I2 are small com

pared to the differences of these 1 quantities for bodytissues. (See also Brooks[ 181. He has concluded thatg(E)=E-2’8.) ) The second part of ( 3 3 ) given by a 2 ( x, y )~ K N ( Egives the Compton contribution scatter to the attenuation. Again a 2 ( x, y ) depends upon the material properties whereas fKN(E), the Klein-Nishina function, describes the(material independent) energy dependence of this contribution. The function fKN(E) is given byc

(33)

D. Statistical Propertiesof Noise in CT ImagesLet again us consider the hypothetical parallel beam of monoenergetic X-ray photons propagating through tissue as shown in Fig. 3 . If the width 7 of the beam is small enough, the integral of p ( x, y ) along the dotted lineisgivenby (2) and (25): Pe(t)E

1+-ln(1 2a

+ 2a) -

(1+ 3 a ) (35) (1+ 2aI2

ray path A B

J

P ( x, y ) dS= ln Ni, -

ln Ne(k7)(41)

with a= E/ 5 10 975. The energy E is in kiloelectronvolts. The importance of ( 3 3 ) lies in the fact that all the energy dependence has been incorporated in known and material independent functions g ( E ) and fKN(E). Substituting this equation in (26) we get

where Ne(k7) denotes the value of N d for the ray at location shown in the figure. Randomness in the measurement of PO ( t ) is introduced by statistical fluctuations inNe(k7). Note that in practice only Ne(k7) is measured directly. The value of Ni, for all rays is inferred by monitoring the X-ray source with a reference detector and from the knowledge of the spatial distribution of emitted X-rays. It is usually safe to N d= s o ( E ) exP[-(AIg(E)+ A~~ K N ( E ) ) I d E( 3 6 ) assume that the reference X-ray flux islarge enough so that Ni, may be considered to be known with negligible error. In where the rest of the discussion here wewill assume that for each rayintegral measurement Ni, is a known deterministic conA1= al(x,y)ds ( 3 7 ) stant, while on the other hand the directly measured quantity ray path Ne(k7) is a random variable. The randomness of Ne(k7) is statistically described by the Poisson probability function and[601, 11081:

(e, k7)as

s

J-

A2

=

ray path

J-

a 2 ( x, Y ) ds.

(38)

A 1 and A2 are, clearly, ray integrals for the functions ( x, y ) a1 and a2 ( x, y ) . Now if we could somehow determineA 1 and A2 for each ray, fromthis information then the functions ( x, y ) a1 and a2 (x, y ) could be separately reconstructed. And, once we know a l ( x, y ) and a 2 ( x, y ), using ( 3 3 ) an attenuation coefficient tomogram could be presented at any energy free from beam hardening artifacts. A few words about the determination of A 1 and A2: Note that it is the intensity Nd that is measured by the detector. Now suppose instead of making one measurement we make two measurements for each ray path for two differentsource spectra. Let us call these measurements I I and Z2, then

where p{ *} denotes the _probability, Ne(k7)! denotes the factorial of Ne(k7), and Ne(k7) the expected value of the measuremen

t

N (k7)= E{Ne(k7)) e

(43)

where E denotesstatisticalexpectation.Notethatthe ance of each measurement is given by variance (Ne (kr)}= Y (k7). o

vari(44)

Because of the randomness in Ne(k7) thetrue value of Pe(k7) will differ from its measured valuewhich will be denoted by P r ( k 7 ) . To bring out this distinction we reexpress

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1258

PROCEEDINGS OF THE IEEE, VOL. 67, NO. 9, SEPTEMBER 1979

and(45)

We will now define a relative-uncertainty image as follows:’ variance{?(x, y ) ) relative-uncertainty at (x, y )= Ni,[?(x, r ) l‘

By interpreting exp[-Pe(kT)] as the probability that (along a ray such as the one shown in Fig. 3) a photon entering the object from side A will emerge (without scattering or absorption) at side B, one can show thatN e ( k T )= N i, exp[ - P e ( k~ )] .-

(55)

(46)

We will now assume that all fluctuations (departures from the mean) in Ne(k7) that have a significant probability of occurrence are much less than the mean. With this assumption and using (41) and (42) it is easily shown that E and

{PF(k~)) PO ( k~ )=

(47)

In computer simulation studies with this definition the relative-uncertainty image becomes independent of the number of incident photons used for measurements, and is completely determined by the choice of the phantom. Fig. 13(c) shows the relative-uncertainty image for the Shepp and Logan phantom (Fig. 13(a)) for Mpr0j= 120 and T= 2/101 and for h ( t ) given by (7). Fig. 13(d) shows graphically the middle horizontal line through Fig. 13(c). Relative-uncertainty at (x, y ) gives us a relative measure of how much confidence an observer might place in the reconstructed value at the point (x, y ) vis-$-vis those elsewhere. We will now show that the results that have been obtained by Shepp and Logan[ 1071 and Chesler[35] are special cases of our (54) or, equivalently (51). Suppose we want to determine the variance of noise at the origin. From (51) we

From the statistical properties of the measured projections P r ( k T ), we will now derive those of the reconstructed image. By combining (8) and (IO), the relationship between the reconstruction at a point (x, y ) and the measured projections is given bynTMproj

can writevariance{P(~,= o))

z Mproj

(56)

?(x, ),

P;(kT) h(x cos BiMmoji=1

+,

sin Bi - kT).(49)

k

where we haveused the fact that h ( t ) is an even function. Chesler et al.[ 3 5] have argued that since h(k7) drops rapidly with k (see (7)), it is safe to make the following approximation for objects that are approximately homogeneous:

Using (47), (48), and (49), we get

which, when T is small enough, may also be written as+ y sin Bi

- kr) ( 5 0 )

and(58)

. h Z ( x cos B i+ y sin Bi

-k

~ ) (51)

where we have used the assumption that fluctuations in P& (~ T )areuncorrelatedfordifferent rays. Equa

tion ( 5 0 ) shows that the expected value of the reconstructed image is equal to that made from the ideal projection data. Before we interpret (51) we will rewrite it as follows. In terms of the ideal projections Pe(kT), we define new projections as

v ( k~= exp[PO(k7)I e )and a new filter function h, ( t ) ash”(t)= h2(t).

(52)

(53)

Note again that i O i ( O ) are the mean number of exiting photons measured forthecenter ray in each projection. Using (58) Chesler et al.[ 3 5] have arrived at the very interesting result that(forthe same uncertaintyinmeasurement)the total number of photons per resolution element required for X-ray CT (using the filtered-backprojection algorithm) is the same as inthe measurement of attenuation of an isolated (excised) piece of tissue with dimensions equal to those of the resolutionelement. Now consider the case where the cross section forwhich the image is being reconstructed is circularly symmetric. The Nei(0)’s for all i’s will be equal, let’s call t i value No. That hs is,

e

variance{f(x,y))=

A

( C o J i*

-

n

ve(kT)ik

h,(x cos Bi

+ y sin d j - kr).

(54)

result only applies when compensators (such as wedges) are not used. They reduce the dynamic range of the detector output signal. In noise andysed their effect can be approximately modeled by using different N~’ for afferent rays. S

‘ hs Ti

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

1259

and Logan headphantom[ 1071. (b) Reconstruction of thephantomfrom 120 projections and 101 raysineach Fig. 13. (a)TheShepp parallel projection. Display matrix: 64 X 64. (c) The relative-uncertainty image for the reconstruction in (b). (d) A graphical depiction of of the relative-uncertainty values through the horizontal middle line (c).

tional to the area under the square of the filter function used for reconstruction. This does not imply that t h i s area could be made arbitrarily small since any major departure from the variance{f(o,O)}= (60) Iwl function, will introduce spatial distortion in the image even thoughit may beless noisy.Equations(60)and(61) By Parseval's theorem this result may be expressed in the fre- were F i t derived by Shepp and Logan[ 1071. None of the domain quency as construed be above equations should the that to imply varis ance approaches zero as r i made arbitrarily small. Note from Fig. 3 that 7 is also the width of themeasurement beam. In variance{ f(0,O))= -"" IH(w)12 n2 dw (61) any practical system as r is reduced, No will decrease also. MprojNo 4/ 2 7 The preceding discussion has resulted in expressions for the where r is the sampling interval for the projection data. This variance of noise in reconstructions made with a filtered-backresult saysthatthe variance of noise at the origin ispropor-projectionalgorithmfor parallel projectiondata. As menh

The expression (5

8) for the variance may now be written as

'-

I

Authorized licensed use limited to: Illinois Institute of Technology. Downloaded on January 30, 2010 at 11:33 from IEEE Xplore. Restrictions apply.

1260

PROCEEDINGS OF THE IEEE,

VOL. 67, NO. 9, SEPTEMBER 1979

tioned before in the introduction filtered-backprojection algorithms have become very popular becauseof their accuracy. Still thequestion arises that given a set of projections can there by an algorithm that might reconstruct an image with a smaller error. The answer t o t h i s question has been supplied by some very theoretically elegantwork by Tretiak[ 1171. Tretiak hasderived an algorithmindependentlowerbound for the mean squared error in a reconstructed image and has argued that for the case of reconstructions from parallel projection data this lower bound is very close to the error estimates obtained by Brooks and DiChiro[ 17] for the filteredbackprojection algorithm, which leads t o the conclusion that very little improvement can be obtained over the performance of such analgorithm. The reader is also referred t o[ 61 and[ 11 1] for discussions on noise properties of the X-ray CT images.

本文来源：https://www.bwwdw.com/article/e8fm.html