Image Reconstruction from Projections | 16

ABSTRACT

Mathematically the fundamental problem in CT imaging is that of estimating

an image or object from its projections or integrals measured at dierent

angles see Section for an intro

duction to this topic A projection of an image is also referred to as the

Radon transform of the image at the corresponding angle after the main pro

ponent of the associated mathematical principles In the continuous

space the projections are ray integrals of the image measured at dierent

ray positions and angles in practice only discrete measurements are avail

able The solution to the problem of image reconstruction from projections

may be formulated variously as completing the corresponding Fourier space

backprojecting and summing the given projection data or solving a set of

simultaneous equations Each of these methods has its own advantages and

disadvantages that determine its suitability to a particular imaging applica

tion In this chapter we shall study the three basic approaches to image

reconstruction from projections mentioned above

Note Most of the derivations presented in this chapter closely follow those

of Rosenfeld and Kak with permission For further details refer to Her

man and Kak and Slaney Parts of this chapter are reproduced

with permission from RM Rangayyan and A Kantzas Image reconstruc

tion Wiley Encyclopedia of Electrical and Electronics Engineering Supple

ment Editor J G Webster Wiley New York NY pp

This material is used by permission of John Wiley Sons Inc

Projection Geometry

Let us consider the problem of reconstructing a D image given parallelray

projections of the image measured at dierent angles Referring to Figure

let fx y represent the density distribution within the image Although

discrete images are used in practice the initial presentation here will be in

continuousspace notations for easier comprehension Consider the ray AB

represented by the equation

x cos y sin t

The integral of fx y along the ray path AB is given by

fx y ds

fx y x cos y sin t

dx dy

where is the Dirac delta function and s x sin y cos The mutually

parallel rays within the imaging plane are represented by the coordinates t s

that are rotated by angle with respect to the x y coordinates as indicated

in Figures and with the s axis being parallel to the rays ds is

thus the elemental distance along a ray When this integral is evaluated for

dierent values of the ray oset t

we obtain the D projection p

t The

function p

t is known as the Radon transform of fx y Note Whereas a

single projection p

t of a D image at a given value of is a D function

a set of projections for various values of could be seen as a D function

Observe that t represents the space variable related to ray displacement along

a projection and not time Because the various rays within a projection are

parallel to one another this is known as parallelray geometry

Theoretically we would need an innite number of projections for all to be

able to reconstruct the image Before we consider reconstruction techniques

let us take a look at the projection or Fourier slice theorem

The Fourier Slice Theorem

The projection or Fourier slice theorem relates the three spaces we encounter

in image reconstruction from projections the image Fourier and projection

Radon spaces Considering a D image the theorem states that the D

Fourier transform of a D projection of the D image is equal to the radial

section slice or prole of the D Fourier transform of the D image at the

angle of the projection This is illustrated graphically in Figure and may

be derived as follows

Let F u v represent the D Fourier transform of fx y given by

F u v

fx y expj ux vy dx dy

Let P

w represent the D Fourier transform of the projection p

t that is

t expj wt dt

where w represents the frequency variable corresponding to t Note If x y s

and t are in mm the units for u v and w will be cyclesmm or mm

Let

FIGURE

Illustration of a ray path AB through a sectional plane or image fx y The

t s axis system is rotated by angle with respect to the x y axis system

ds represents the elemental distance along the ray path AB p

is the ray

integral of fx y for the ray path AB p

t is the parallelray projection

Radon transform or integral of fx y at angle See also Figures and

Adapted with permission from A Rosenfeld and AC Kak Digital

Picture Processing nd ed New York NY

Academic Press

FIGURE

Illustration of the Fourier slice theorem F u v is the D Fourier trans

form of fx y F w

w is the D Fourier transform of p

F w

w is the D Fourier transform of p

t Reproduced with

permission from RM Rangayyan and A Kantzas Image reconstruction

Wiley Encyclopedia of Electrical and Electronics Engineering Supplement

Editor J G Webster Wiley New York NY pp

This

material is used by permission of John Wiley Sons Inc

ft s represent the image fx y rotated by angle with the transformation

given by

cos sin

sin cos

Then

ft s ds

t expj wt dt

ft s ds

expj wt dt

Transforming from t s to x y we get

fx y expj wx cos y sin dx dy

F u v for u w cos v w sin

F w

which expresses the projection theorem Observe that t x cos y sin and

dx dy ds dt

It immediately follows that if we have projections available at all angles from

we can take their D Fourier transforms ll the D Fourier space

with the corresponding radial sections or slices and take an inverse D Fourier

transform to obtain the image fx y The diculty lies in the fact that in

practice only a nite number of projections will be available measured at

discrete angular positions or steps Thus some form of interpolation will be

essential in the D Fourier space Extrapolation may also be required

if the given projections do not span the entire angular range This method

of reconstruction from projections known as the Fourier method succinctly

relates the image Fourier and Radon spaces The Fourier method is the most

commonly used method for the reconstruction of MR images

A practical limitation of the Fourier method of reconstruction is that inter

polation errors are larger for higher frequencies due to the increased spacing

between the samples available on a discrete grid Samples of P

w computed

from p

t will be available on a polar grid whereas the D Fourier transform

F u v andor the inversetransformed image will be required on a Carte

sian rectangular grid This limitation could cause poor reconstruction of

highfrequency sharp details

Backprojection

Let us now consider the simplest reconstruction procedure backprojection

BP Assuming the rays to be ideal straight lines rather than strips of nite

width and the image to be made of dimensionless points rather than pixels

or voxels of nite size it can be seen that each point in the image fx y

contributes to only one ray integral per parallelray projection p

t with

t x cos y sin We may obtain an estimate of the density at a point by

simply summing integrating all rays that pass through it at various angles

that is by backprojecting the individual rays In doing so however the

contributions to the various rays of all of the other points along their paths

are also added up causing smearing or blurring yet this method produces a

reasonable estimate of the image Mathematically simple BP can be expressed

fx y

t d where t x cos y sin

This is a sinusoidal path of integration in the t Radon space In practice

only a nite number of projections and a nite number of rays per projection

will be available that is the t space will be discretized hence interpola

tion will be required

Examples of reconstructed images Figure a shows a synthetic

D image phantom which we will consider to represent a crosssectional

plane of a D object The objects in the image were dened on a discrete

grid and hence have step andor jagged edges Figure a is a plot of the

projection of the phantom image computed at

observe that the values

are all positive

a b

FIGURE

a A synthetic D image phantom with eightbit pixels repre

senting a crosssection of a D object b Reconstruction of the phantom

in a obtained using projections from

in steps of

with the

simple BP algorithm Reproduced with permission from RM Rangayyan

and A Kantzas Image reconstructionWiley Encyclopedia of Electrical and

Electronics Engineering Supplement Editor J G Webster Wiley New

York NY pp

This material is used by permission of John

Wiley Sons Inc

Figure b shows the reconstruction of the phantom obtained using

projections from

in steps of

with the simple BP algorithm

While the objects in the image are faintly visible the smearing eect of the

BP algorithm is obvious

Considering a point source as the image to be reconstructed it becomes

evident that BP produces a spokelike pattern with straight lines at all pro

jection angles intersecting at the position of the point source This may be

considered to be the PSF of the reconstruction process which is responsible

for the blurring of details

FIGURE

a Projection of the phantom image in Figure a computed at

b Filtered version of the projection using only the ramp lter inherent to the

FBP algorithm Reproduced with permission from RM Rangayyan and A

Kantzas Image reconstruction Wiley Encyclopedia of Electrical and Elec

tronics Engineering Supplement Editor J G Webster Wiley New York

NY pp

This material is used by permission of John Wiley

Sons Inc

The use of limited projection data in reconstruction results in geometric dis

tortion and streaking artifacts The distortion may

be modeled by the PSF of the reconstruction process if it is linear and shift

invariant this condition is satised by the BP process The PSFs of the simple

BP method are shown as images in Figure a for the case with projec

tions over

and in Figure b for the case with projections from

The reconstructed image is given by the convolution of the original

image with the PSF the images in parts c and d of Figure illustrate the

corresponding reconstructed images of the phantom in Figure a Lim

ited improvement in image quality may be obtained by applying deconvolution

lters to the reconstructed image De

convolution is implicit in the ltered convolution backprojection technique

which is described next

Filtered backprojection

Consider the inverse Fourier transform relationship

fx y

F u v expj ux vy du dv

Changing from the Cartesian coordinates u v to the polar coordinates w

where w

and tan

vu we get

fx y

F w expj wx cos y sin w dw d

F w

exp fj wx cos y sin gw dw d

Here u w cos v w sin and du dv w dw d Because F w

F w we get

fx y

F w jwj expj wt dw

wjwj expj wt dw

with t x cos y sin as before If we dene

wjwj expj wt dw

we get

fx y

t d

x cos y sin d

a b

c d

FIGURE

PSF of the BP procedure using a projections from

in steps

b projections from

in steps of

The images a

and b have been enhanced with Reconstruction of the phantom in

Figure a obtained using c projections as in a with the BP algo

rithm d projections as in b with the BP algorithm Reproduced with

permission from RM Rangayyan and A Kantzas Image reconstruction

Wiley Encyclopedia of Electrical and Electronics Engineering Supplement

Editor J G Webster Wiley New York NY pp

This

material is used by permission of John Wiley Sons Inc

It is now seen that a perfect reconstruction of fx y may be obtained by

backprojecting ltered projections q

t instead of backprojecting the original

projections p

t hence the name ltered backprojection FBP The lter is

represented by the jwj function known as the ramp lter see Figure

Observe that the limits of integration in Equation are for and

for w In practice a smoothing window should be applied to reduce

the amplication of highfrequency noise by the jwj function Furthermore

the integrals change to summations in practice due to the nite number of

projections available as well as the discrete nature of the projections them

selves and of the Fourier transform computations employed Details of the

discrete version of FBP are provided in the next section

An important feature of the FBP technique is that each projection may be

ltered and backprojected while further projection data are being acquired

which was of help in online processing with the rstgeneration CT scanners

see Figure Furthermore the inverse Fourier transform of the lter

jwj with modications to account for the discrete nature of measurements

smoothing window etc see Figure could be used to convolve the projec

tions directly in the t space using fast array processors FBP is the most

widely used procedure for image reconstruction from projections however

the procedure provides good reconstructed images only when a large number

of projections spanning the full angular range of

are available

Discrete ltered backprojection

The ltering procedure with the jwj function in theory must be performed

over w In practice the signal energy above a certain frequency

limit W will be negligible and jwj ltering beyond the limit will only amplify

noise Thus we may consider the projections to be bandlimited to W Then

using the sampling theorem p

t can be represented by its samples at the

sampling rate W as

sin W t

W t

Then

exp

j w

where

w if jwj W

otherwise

FIGURE

a The bandlimited ramp lter inherent to the FBP algorithm b The

ramp lter weighted with a Butterworth lowpass lter The lters are shown

for both positive and negative frequency along the w axis with w at the

center The corresponding lters are shown in the Radon domain t space in

Figure

FIGURE

a The inverse Fourier transform of the bandlimited ramp lter inherent to

the FBP algorithm in the t space b The inverse Fourier transform of the

ramp lter weighted with a Butterworth lowpass lter The functions which

are used for convolution with projections in the CBP method are shown for

both t and t with t at the center The corresponding lters are shown

in the frequency domain w space in Figure

If the projections are of nite order that is they can be represented by a

nite number of samples N then

exp

j w

Let us assume that N is even and let the frequency axis be discretized as

w k

for k

Then

exp

This represents a DFT relationship and may be eval

uated using the FFT algorithm

The ltered projection q

t may then be obtained as

wjwj expj wtdw

exp

j k

If we want to evaluate q

t for only those values of t at which p

t has been

sampled we get

exp

In order to control noise enhancement by the jk

j lter it may be bene

cial to include a lter window such as the Hamming window then

exp

with

cos

Using the convolution theorem we get

where denotes circular periodic convolution and g

is the inverse DFT of

Butterworth or other lowpass

lters may also be used instead of the Hamming window

Observe that the inverse Fourier transform of jwj does not exist due to

the fact that jwj is neither absolutely nor square integrable However if we

consider the inverse Fourier transform of jwj expjwj as we get the

function

for large t p

The reconstructed image may be obtained as

fx y

x cos

y sin

where the L angles

are those at which the projections p

t are available

and the x y coordinates are discretized as appropriate

For practical implementation of discrete FBP let us consider the situation

where the projections have been sampled with an interval of mm with no

aliasing error Each projection p

m is then limited to the frequency band

W W with W

cyclesmm The continuous versions of the ltered

projections are

wHw expj wtdw

where the lter Hw jwjb

w with b

w as dened in Equation

The impulse response of the lter Hw is

sin t

sint

Because we require ht only at integral multiples of the sampling interval

we have

n even

n odd

The ltered projections q

m may be obtained as

nhm n m N

where N is the nite number of samples in the projection p

m Observe

that hn is required for n N N When the l

ter is implemented as a convolution the FBP method is also referred to as

convolution backprojection CBP

The procedure for FBP may be expressed in algorithmic form as follows

Measure the projection p

Compute the ltered projection q

Backproject the ltered projection q

Repeat Steps for all projection angles

l L

The FBP algorithm is suitable for online implementation in a translate

rotate CT scanner because each parallelray projection may be ltered and

backprojected as soon as it is acquired while the scanner is acquiring the next

projection The reconstructed image is ready as soon as the last projection

is acquired ltered and backprojected If the projections are acquired using

fanbeam geometry see Figure one could either rebin the fanbeam

data to compose parallelray projections or use reconstruction algorithms

specically tailored to fanbeam geometry

Examples of reconstructed images Figure b shows a plot of the

ltered version of the projection in Figure a using only the ramp lter

jwj inherent to the FBP algorithm Observe that the ltered projection has

negative values

Figure a shows the reconstruction of the phantom in Figure a

obtained using projections with the FBP algorithm only the ramp lter

that is implicit in the FBP process was used with no other smoothing or low

pass lter function The contrast and visibility of the objects are better than

those in the case of the simple BP result in Figure b however the image

is noisy due to the increasing gain of the ramp lter at higher frequencies

The reconstructed image also exhibits artifacts related to the computation of

the projections on a discrete grid refer to Herman Herman et al

and Kak and Slaney for discussions on this topic The use of additional

lters could reduce the noise and artifacts Figure b shows the result of

reconstruction with the FBP algorithm including a fourthorder Butterworth

lter with the dB cuto at times the maximum frequency present in

the data to lter the projections The Butterworth lter has suppressed the

noise and artifacts at the expense of blurring the edges of the objects in the

image

a b

FIGURE

a Reconstruction of the phantom in Figure a obtained using pro

jections from

in steps of

with the FBP algorithm only the ramp

lter that is implicit in the FBP process was used b Reconstruction of the

phantom with the FBP algorithm as in a but with the additional use of

a Butterworth lowpass lter See also Figure Reproduced with permis

sion from RM Rangayyan and A Kantzas Image reconstruction Wiley

Encyclopedia of Electrical and Electronics Engineering Supplement Editor

J G Webster Wiley New York NY pp

This material is

used by permission of John Wiley Sons Inc

The Radon transform may be interpreted as a transformation of the given

image from the x y space to the t space In practical CT scanning the

projection or rayintegral data are obtained as samples at discrete intervals in t

and Just as we encounter the Nyquist or Shannon sampling theorem in the

representation of a D signal in terms of its samples in time we now encounter

the requirement to sample adequately along both the t and axes A major

distinction lies in the fact that the measurements made in CT scanners are

discrete to begin with and the signal the body or object being imaged

cannot be preltered to prevent aliasing Undersampling in either axis will

lead to aliasing errors and poor reconstructed images

Figure a shows the reconstructed version of the phantom in Fig

ure a obtained using only projections spanning the

range

in sampling steps of

and using the FBP algorithm Although the edges

of the objects in the image are sharper than those in the reconstruction ob

tained using the BP algorithm with the same parameters see Figure c

the image is aected by severe streaking artifacts due to the limited num

ber of projections used Figure b shows the reconstructed image of the

phantom obtained using projections but spanning only the angular range

in steps of

The limited angular coverage provided by the

projections has clearly aected the quality of the image and has introduced

geometric distortion

Algebraic Reconstruction Techniques

The algebraic reconstruction technique ART is related to the

Kaczmarz method of projections for solving simultaneous equations

The Kaczmarz method takes an approach that is completely dierent from

that of the Fourier or FBP methods the available projections treated

as individual ray sums in a discrete representation are seen as a set of

simultaneous equations with the unknown quantities being the discrete pixels

of the image The large size of images encountered in practice precludes the

use of the usual methods for solving simultaneous equations Furthermore

in many practical applications the number of available equations may be far

less than the number of pixels in the image to be reconstructed the set of

simultaneous equations is then underdetermined The Kaczmarz method of

projections is an elegant iterative method that may be implemented easily

Note The Kaczmarz method uses the term projection in the vectorial

or geometric sense and individual ray sums are processed one at a time

Observe that a set of ray sums or integrals is also known as a projection The

distinction should be clear from the context

a b

FIGURE

a Reconstruction of the phantom in Figure a obtained using projec

tions from

in steps of

with the FBP algorithm b Reconstruc

tion of the phantom obtained using projections from

in steps

with the FBP algorithm The ramp lter that is implicit in the FBP

process was combined with a Butterworth lowpass lter in both cases See

also Figure Reproduced with permission from RM Rangayyan and A

Kantzas Image reconstruction Wiley Encyclopedia of Electrical and Elec

tronics Engineering Supplement Editor J G Webster Wiley New York

NY pp

This material is used by permission of John Wiley

Sons Inc

Let the image to be reconstructed be divided into N cells f

denoting the

value in the n

cell The image density or intensity is assumed to be constant

within each cell Let M be the number of ray sums available expressed as

m M

where w

is the contribution factor of the n

image element to the m

ray

sum equal to the fractional area of the n

cell crossed by the m

ray path

as illustrated in Figure Note An image and its Radon transform are

each represented using only one index in this formulation of ART Observe

that for a given ray m most of w

will be zero because only a few elements

of the image contribute to the corresponding ray sum Equation may

also be expressed as

A grid representation with N cells gives the image N degrees of freedom

Thus an image represented by f f

may be considered to

be a single point in an N dimensional hyperspace Then each of the above

raysum equations will represent a hyperplane in this hyperspace If a unique

solution exists it is given by the intersection of all the hyperplanes at a single

point To arrive at the solution the Kaczmarz method takes the approach of

successively and iteratively projecting an initial guess and its successors from

one hyperplane to the next

Let us for simplicity consider a D version of the situation with N M

as illustrated in Figures and Let f

represent vectorially the

initial guess to the solution and let w

represent vectorially

the series of weights coecients in the rst ray equation The rst ray sum

may then be written as

f p

The hyperplane represented by this equation is orthogonal to w

Consider

two images or points f

and f

belonging to the hyperplane We have w

and w

Hence w

Therefore w

is orthogonal to

the hyperplane

With reference to Figure Equation indicates that for the vector

OC corresponding to any point C on the hyperplane its projection on to the

vector w

is of a constant length The unit vector OU along w

is given by

FIGURE

ART treats the image as a matrix of discrete pixels of nite size xy

Each ray has a nite width The fraction of the area of the n

pixel

crossed by the m

ray is represented by the weighting factor w

area of ABCDxy for the n

pixel f

in the gure Adapted with

permission from A Rosenfeld and AC Kak Digital Picture Processing nd

ed New York NY

Academic Press

The perpendicular distance of the hyperplane from the origin is

kOAk OU OC

Now

and

kGHk kOFk kOAk f

OU kOAk

Because the directions of GH and OU are the same GH kGHkOU

Thus

Therefore

In general the m

estimate is obtained from the m

estimate as

That is the m

estimate on hand is projected on to the hyperplane of

them

ray sum and the deviation from the true ray sum p

is obtained This

deviation is normalized and applied as a correction to all the pixels according

to the weighting factors in w

When this process is applied to all theM ray

sum hyperplanes given one cycle or iteration is completed Note Because

the image is updated by altering the pixels along each individual ray sum the

index of the updated estimate or of the iteration is equal to the index of the

latest ray sum used However as the entire process is iterated the index of

the estimate is reset at the beginning of each iteration

Depending upon the initial guess and the organization of the hyperplanes

a number of iterations may have to be completed in order to obtain the solu

tion if it exists The following important characteristics of ART are worth

observing

ART proceeds ray by ray and is iterative

FIGURE

Illustration of the Kaczmarz method of solving a pair of simultaneous equa

tions in two unknowns The solution is f

The weight vectors

for the two ray sums straight lines are w

and w

The equations of the straight lines are w

f f

and

f f

The initial estimate is f

The

rst updated estimate is f

the second updated estimate is

Because two ray sums are given two corrections constitute

one cycle or iteration of ART The path of the second cycle of ART is also

illustrated in the gure Reproduced with permission from RM Rangayyan

and A Kantzas Image reconstructionWiley Encyclopedia of Electrical and

Electronics Engineering Supplement Editor J G Webster Wiley New

York NY pp

This material is used by permission of John

Wiley Sons Inc

FIGURE

Illustration of the algebraic reconstruction technique f

is an improved

estimate computed by projecting the initial guess f

on to the hyperplane

the straight line AG in the illustration corresponding to the rst ray sum

given by the equation w

f p

Adapted with permission from A Rosen

feld and AC Kak Digital Picture Processing nd ed New York NY

Academic Press

If the hyperplanes of all the given ray sums are mutually orthogonal

we may start with any initial guess and reach the solution in only one

cycle

On the other hand if the hyperplanes subtend small angles with one

another a large number of iterations will be required The number of

iterations may be reduced by using optimized rayaccess schemes

If the number of ray sums is greater than the number of pixels that is

M N but the measurements are noisy no unique solution exists

the procedure will oscillate in the neighborhood of the intersections of

the hyperplanes

If M N the system is underdetermined and an indenite or innite

number of partial solutions exist It has been shown that unconstrained

ART converges to the minimumvariance estimate

The major advantage of ART is that any a priori information available

about the image may be introduced easily into the iterative procedure

for example upper andor lower limits on pixel values and the spatial

boundaries of the image This may help in obtaining a useful solution

even if the system is underdetermined

Approximations to the Kaczmarz method

We could rewrite the reconstruction step in Equation at the n

pixel

level as

where q

This equation indicates that

when we project the m

estimate on to the m

hyperplane the cor

rection factor for the n

cell is

Here p

is the given true ray sum for the m

ray and q

is the computed

ray sum for the same ray for the estimated image on hand p

is the

error in the estimate which is normalized and applied as a correction to all

the pixels with appropriate weighting Because the correction factor is added

to the current image this version of ART is known as additive ART

In one of the approximations to Equation the weights w

are simply

replaced by zeros or ones depending upon whether the center of the n

image

cell is within the m

ray of nite width or not Then the coef

cients need not be computed and stored we may instead determine the

pixels to be corrected for the ray considered during the reconstruction pro

cedure Furthermore it follows that

the number of pixels

crossed by the m

ray The correction applicable to all of the pixels along

the m

ray is p

Then

Because the corrections could be negative negative pixel values may be en

countered Because negative values are not meaningful in most imaging ap

plications the constrained and thereby nonlinear version of ART is dened

max

The corrections could also be multiplicative

This version of ART is known as multiplicative ART In this case no pos

itivity constraint is required Furthermore the convex hull of the image is

almost guaranteed subject to approximation related to the number of ray

sums available and their angular coverage because a pixel once set to zero

will remain so during subsequent iterations It has been shown that the mul

tiplicative version of ART converges to the maximumentropy estimate of the

image

A generic ART procedure may be expressed in the following algorithmic

form

Prepare an initial estimate f

of the image All of the pixels in the

initial image could be zero for additive ART however for multiplicative

ART pixels within at least the convex hull of the object in the image

must have values other than zero

Compute the ray sum q

for the rst ray path m for the estimate

of the image on hand

Obtain the dierence between the true ray sum p

and the computed

ray sum q

and apply the correction to all the pixels belonging to the

ray according to one of the ART equations for example Equation

or Apply constraints if any based upon the a priori

information available

Perform Steps and for all rays available m M

Steps constitute one cycle or iteration over all available ray sums

Repeat Steps as many times as required If desired compute a

measure of convergence such as

Stop if the error or dierence is less than a prespecied limit else go

back to Step

For improved convergence a simultaneous correction procedure Simulta

neous Iterative Reconstruction Technique SIRT has been proposed

where the corrections to all the pixels from all the rays are rst computed

and the averaged corrections are applied at the same time to all the pixels

that is only one correction is applied per pixel per iteration Guan and

Gordon proposed dierent rayaccess schemes to improve convergence

including the consecutive use of rays in mutually orthogonal directions In

the absence of complete projection data spanning the full angular range of

ART typically yields better results than the FBP or the Fourier

methods

Examples of reconstructed images Figure b shows the recon

struction of the phantom shown in part a of the same gure obtained us

ing parallelray projections with three iterations of constrained additive

ART as in Equation Ray sums for use with ART were computed

from the phantom image data using angledependent ray width given by

maxj sinj j cosj The ART reconstruction is better than that

given by the BP or the FBP algorithm The advantages of constrained ART

due to the use of the positivity constraint that is the a priori knowledge

imposed and of the ability to iterate are seen in the improved quality of

the result

Figure shows reconstructed images of the phantom obtained using

only projections spanning the angular ranges of a

in steps

and b

in steps of

respectively The limited number

of projections used and the limited angular coverage of the projections in

the second case have aected the quality of the reconstructed images and

introduced geometric distortion However when compared with the results

of BP and FBP with similar parameters see Figures and ART has

provided better results

a b

FIGURE

a A synthetic D image phantom with eightbit pixels repre

senting a crosssection of a D object the same image as in Figure a

b Reconstruction of the phantom obtained using projections from

in steps of

with three iterations of constrained additive ART See

also Figure Reproduced with permission from RM Rangayyan and A

Kantzas Image reconstruction Wiley Encyclopedia of Electrical and Elec

tronics Engineering Supplement Editor J G Webster Wiley New York

NY pp

This material is used by permission of John Wiley

Sons Inc

a b

FIGURE

Reconstruction of the phantom in Figure a obtained using a pro

jections from

in steps of

with three iterations of constrained

additive ART b projections from

in steps of

with three

iterations of constrained additive ART See also Figures and Re

produced with permission from RM Rangayyan and A Kantzas Image

reconstruction Wiley Encyclopedia of Electrical and Electronics Engineer

ing Supplement Editor J G Webster Wiley New York NY pp

This material is used by permission of John Wiley Sons Inc

Imaging with Diracting Sources

In some applications of CT imaging such as imaging pregnant women Xray

imaging may not be advisable Imaging with nonionizing forms of radiation

such as acoustic ultrasonic and electromagnetic optical or thermal

imaging is then a valuable alternative Xray imaging is also not suitable

when the object to be imaged has poor contrast in density or atomic number

distribution An important point to observe in acoustic or electromagnetic

imaging is that these forms of energy do not propagate along straightline ray

paths through a body due to refraction and diraction When the dimensions

of the inhomogeneities in the object being imaged are comparable to or smaller

than the wavelength of the radiation used geometric propagation concepts

cannot be applied it becomes necessary to consider wave propagation and

diractionbased methods

When the body being imaged may be treated as a weakly scattering object

in the D sectional plane and invariant in the axial direction the Fourier

diraction theorem is applicable This theorem states that the D Fourier

transform of a projection including the eects of diraction gives values of

the D Fourier transform of the image along a semicircular arc Interpolation

methods may be developed in the Fourier space taking this property into

account for reconstruction of images from projections obtained with diracting

sources Backpropagation and algebraic techniques have also been proposed

for the case of imaging with diracting sources

Display of CT Images

Xray CT is capable of producing images with high density resolution on the

order of one part in For display purposes the attenuation coecients

are normalized with respect to that of water and expressed as

HU K

where is the measured attenuation coecient and

is the attenuation

coecient of water The K parameter used to be set at in early models

of the CT scanner It is now common to use K to obtain the CT

number in Hounseld units HU named after the inventor of the rst

commercial medical CT scanner This scale results in values of about

for bone for water about for air to for soft tissue

and about for lung tissue Table shows the mean and SD of the

CT values in HU for several types of tissue in the abdomen

TABLE

Mean and SD of CT Values in

Hounseld Units HU for a Few

Types of Abdominal Tissue

Tissue Mean HU SD

Air

Fat

Bile

Kidney

Pancreas

Blood aorta

Muscle

Necrosis

Spleen

Liver

Viable tumor

Marrow

Calcication

Bone

Based upon Mategrano et al

Estimated from CT exams with con

trast see Section based upon pixels in each category The

contrast medium is expected to increase the CT values of vascularized tissues

by HU The CT number for air should be the estimated

value is slightly dierent due to noise in the images Reproduced with per

mission from FJ Ayres MK Zuo RM Rangayyan GS Boag V Odone

Filho and M Valente Estimation of the tissue composition of the tumor

mass in neuroblastoma using segmented CT images Medical and Biological

Engineering and Computing

IFMBE

The dynamic range of CT values is much wider than those of common dis

play devices and that of the human visual system at a given level of adaptation

Furthermore clinical diagnosis requires detailed visualization and analysis of

small density dierences For these reasons the presentation of the entire

range of values available in a CT image in a single display is neither practi

cally feasible nor desirable In practice small windows of the CT number

scale are selected and linearly expanded to occupy the capacity of the display

device The window width and level center values may be chosen inter

actively to display dierent density ranges with improved perceptibility of

details within the chosen density window Values above or below the window

limits are displayed as totally white or black respectively This technique

known as windowing or density slicing may be expressed as

gx y

if fx y m

fx ym if m fx y M

N if fx y M

where fx y is the original image in CT numbers gx y is the windowed

image to be displayed mM is the range of CT values in the window to be

displayed and N is the range of the display values The window width is

M m and the window level or center is M m the display range is

typically with bit display systems

Example Figure shows a set of two CT images of a patient with

head injury with each image displayed using two sets of window level and

width The eects of the density window chosen on the features of the image

displayed are clearly seen in the gure either the fractured bone or the brain

matter are seen in detail in the windowed images but not both in the same

image See Figures and for more examples of density windowing

in the display of CT images See also Figures and as well

as Section for more examples of Xray CT images Examples of images

reconstructed from projection data from other modalities of medical imaging

such as MRI SPECT and PET are provided in Sections and

A dramatic visualization of details may be achieved with pseudocolor tech

niques Arbitrary or structured color scales could be assigned to CT values

by LUTs or grayscaletocolor transformations Some of the popular color

transforms are the rainbow VIBGYOR violet indigo blue green yel

low orange red and the heated metal color black red yellow white

sequences Diculties may arise however in associating density values with

dierent colors if the transformation is arbitrary and not monotonic in in

tensity or total brightness Furthermore small changes in CT values could

cause abrupt changes in the corresponding colors displayed especially with

a mapping such as VIBGYOR An LUT linking the displayed colors to CT

numbers or other pixel attributes may assist in improved visual analysis of

image features in engineering and scientic applications

a b

c d

FIGURE

Two CT images of a patient with head injury with each image displayed

with two sets of window level and width Images a and b are of the same

section images c and d are of another section The window levels used to

obtain images a d are and HU respectively the window

widths used are and HU respectively The windows in b

and d display the skull the fracture and the bone segments but the brain

matter is not visible the windows a and c display the brain matter in

detail but the fracture area is saturated Images courtesy of W Gordon

Health Sciences Centre Winnipeg MB Canada

Agricultural and Forestry Applications

Cruvinel et al and Vaz et al developed a portable Xray and

gammaray minitomograph for application in soil science and used the scan

ner to measure water content and bulk density of soil samples Soilrelated

studies address the identication of features such as fractures wormholes and

roots and assist in studies of ow of various contaminants in soil

Forestry applications of CT have appeared in the literature in the form

of scanning of live trees to measure growth rings and detect decay using a

portable Xray CT scanner and monitoring tree trunks or logs in the

timber and lumber industry Figures and show a portable CT scan

ner in operation and images of a utility pole

FIGURE

A portable CT scanner to image live trees Reproduced with permission from

AM Onoe JW Tsao H Yamada H Nakamura J Kogure H Kawamura

and M Yoshimatsu Computed tomography for measuring annual rings of a

live tree Proceedings of the IEEE

IEEE

FIGURE

CT images of a Douglas r utility pole and photographs of the correspond

ing physical sections The sections in a demonstrate normal annual growth

rings The sections in b indicate severe decay in the heartwood region

Reproduced with permission from AM Onoe JW Tsao H Yamada H

Nakamura J Kogure H Kawamura and M Yoshimatsu Computed to

mography for measuring annual rings of a live tree Proceedings of the IEEE

IEEE

Microtomography

The resolution of common CT devices used in medical and other applications

varies from the common gure of mm to about m in cross

section and mm between slices Special systems have been built to image

small samples of the order of cm

in volume with resolution of the order of

m in crosssection

See Section for illustrations of the LSF and MTF of a CT system Such

an imaging procedure is called microtomography microCT or CT being a

hybrid of tomography and microscopy Most CT studies are performed with

nely focused and nearly monochromatic Xray beams produced by a particle

accelerator such as a synchrotron Stock provides a review of the basic

principles techniques and applications of CT Whereas most CT studies

have been limited to small excised samples Sasov discusses the design

of a CT to image whole small animals with resolution of the order of m

Umetani et al present the application of synchrotronbased CT in

a microangiography mode to the study of microcirculation within the heart

of small animals Shimizu et al studied the eects of alveolar hem

orrhage and alveolitis pneumonia on the microstructure of the lung using

synchrotronbased CT Johnson et al developed a microfocal Xray

imaging system to image the arterial structure in the rat lung and studied

the decrease in distensibility due to pulmonary hypertension

Shaler et al applied CT to study the ber orientation and void

structure in wood paper and wood composites Illman and Dowd

studied the xylem tissue structure of wood samples and analyzed the loss

of structrural integrity due to fungal degradation

Example of application to the analysis of bone structure Injury

to the anterior cruciate ligament is expected to lead to a decrease in bone

mineral density Posttraumatic osteoarthritis is known to cause joint space

narrowing and fullthickness cartilage erosion It has been established that

the cancellous architecture of bone is related to its mechanical properties and

strength

Conventional studies on bone structure have been performed through histo

logical analysis of thin slices of bone samples Spatial resolution of the order

of m can be realized with optical microscopy or microradiography How

ever in addition to being destructive this procedure could cause artifacts due

to the slicing operation Furthermore limitations exist in the D information

derived from a few slices Boyd et al applied CT techniques to analyze

the morphometric and anisotropic changes in periarticular cancellous bone

due to ligament injury and osteoarthritis

In the work of Boyd et al anterior cruciate ligaments of dogs were tran

sected by arthrotomy After a specic recovery period at euthanasia the

femora and tibia were extracted Cylindrical bone cores of diameter mm

and length mm were obtained from the weightbearing regions of

the medial femoral condyles and the medial tibial plateaus Contralateral

bone core samples were also extracted and processed to serve as internal con

trols The cores were scanned at a resolution of m and D images with

diameter of voxels and length of up to voxels were reconstructed

Morphometric parameters such as bone volume ratio relative surface density

and trabecular thickness were estimated from the images It was postulated

that the anisotropy of the bone fabric or trabecular orientation is related to

mechanical anisotropy and loading conditions

Figure shows two sample bonecore sectional images each of a case

of transected anterior cruciate ligament of a dog after weeks of recovery

and the corresponding contralateral control sample Figure shows D

renditions of the same samples The bone core related to ligament transection

demonstrates increased trabecular spacing and hence lower bone density than

the contralateral sample Boyd et al observed that signicant periarticular

bone changes occur as early as three weeks after ligament transection the

changes were more pronounced weeks after transection

a b

FIGURE

a Two sample sectional images of a bone core sample in a case of anterior

cruciate ligament transection after weeks of recovery b Two sample

sectional images of the corresponding contralateral bone core sample The

diameter of each sample is mm See also Figure Images courtesy of

SK Boyd University of Calgary

Example of application to the study of microcirculation in the

heart Umetani et al developed a CT system using monochromatized

synchrotron radiation for use as a microangiography tool to study circulatory

disorders and earlystage malignant tumors Two types of detection systems

were used an indirect system including a uorescent screen optical coupling

and a CCD camera and a direct system with a beryllium faceplate a photo

conductive layer and an electronbeam scanner

FIGURE

D renditions of CT reconstructions of a bone core sample in a case of

anterior cruciate ligament transection after weeks of recovery left and

the corresponding contralateral bone core sample right The diameter of

each sample is mm See also Figure Images courtesy of SK Boyd

University of Calgary

In one of the experiments of Umetani et al the left side of the heart of a

rat was xed in formalin after barium sulphate was injected into the coronary

artery Figure a shows a projection image a microradiograph of the

specimen obtained at a resolution of m The image clearly shows the left

anteriordescending coronary artery Figure b shows a D visualization

of a part of the specimen of diameter mm and height mm reconstructed

at a resolution of m The D structure of small blood vessels with diameter

of the order of m was visualized by this method Visualization of

tumorinduced small vessels that feed lesions was found to be useful in the

diagnosis of malignant tumors

Application Analysis of the Tumor in Neuroblasto

Neuroblastoma

Neuroblastoma is a malignant tumor of neuralcrest origin that may arise

anywhere along the sympathetic ganglia or within the adrenal medulla

There are three types of ganglion cell lesions that form a spectrum of

neoplastic disease Neuroblastoma is the most immature and malignant form

of the three usually presenting before the age of ve years Ganglioneuro

blastoma is a more mature form that retains some malignant characteristics

with peak incidence between ve and years of age Ganglioneuroma is well

dierentiated and benign typically presenting after years of age

Neuroblastoma is the most common extracranial solid malignant tumor

in children it is the third most common malignancy of childhood It

accounts for ! of all childhood cancers and ! of all deaths

related to cancer in the pediatric age group The median age at diagnosis is

two years and ! of the diagnosed cases are in children under the age of

ve years

In the US about children and adolescents younger than years of

age are diagnosed with neuroblastoma every year Neuroblastoma is the

most common cancer of infancy with an incidence rate that is almost double

that of leukemia the next most common malignancy occurring during the

rst year of life The rate of incidence of neuroblastoma among infants

in the US has increased from per million in the period to

per million in the period Gurney et al estimated an

annual increase rate of ! for extracranial neuroblastoma Although some

countries instituted screening programs to detect neuroblastoma in infants

studies have indicated that the possible benet of screening on mortality is

small and has not yet been demonstrated in reliable data

FIGURE

a Projection image of a rat heart specimen obtained using synchrotron ra

diation b D CT image of the rat heart specimen Reproduced with per

mission from K Umetani N Yagi Y Suzuki Y Ogasawara F Kajiya T

Matsumoto H Tachibana M Goto T Yamashita S Imai and Y Kajihara

Observation and analysis of microcirculation using highspatialresolution

image detectors and synchrotron radiation Proceedings SPIE Medical

Imaging Physics of Medical Imaging pp

SPIE

Sixtyve percent of the tumors related to neuroblastoma are located in the

abdomen approximately twothirds of these arise in the adrenal gland Fif

teen percent of neuroblastoma are thoracic usually located in the sympathetic

ganglia of the posterior mediastinum Ten to twelve percent of neuroblastoma

are disseminated without a known site of origin

Staging and prognosis The most recent staging system for neuroblas

toma is the International Neuroblastoma Staging System which takes into

account radiologic ndings surgical resectability lymphnode involvement

and bonemarrow involvement The staging ranges from Stage for a

localized tumor with no lymphnode involvement to Stage with the disease

spread to distant lymph nodes bone liver and other organs

The main determinant factors for prognosis are the patient"s age and the

stage of the disease The survival rate of patients diagnosed with

neuroblastoma under the age of one year is ! whereas it is only ! for

patients over the age of three years Whereas the survival rate of patients

diagnosed with Stage neuroblastoma is in the range ! it is only

! for patients diagnosed with Stage of the disease

The site of the primary tumor is also said to be of relevance in the overall

prognosis Tumors arising in the abdomen and pelvis have the worst prognosis

with adrenal tumors having the highest mortality Thoracic neuroblastoma

has a better overall survival rate ! than abdominal tumors !