Deconvolution, Deblurring, and Restoration | 17

ABSTRACT

Deconvolution Deblurring and Restoration

Image enhancement techniques are typically designed to yield better look

ing images satisfying some subjective criteria In comparison with the given

image the processed image may not be closer to the true image in any sense

On the other hand image restoration

is dened as image quality improvement under objective evaluation

criteria to nd the best possible estimate of the original unknown image from

the given degraded image The commonly used criteria are LMS MMSE and

distance measures of several types Additional constraints based upon prior

and independent knowledge about the original image may also be imposed to

limit the scope of the solution Image restoration may then be posed as a

constrained optimization problem

Image restoration requires precise information about the degrading phe

nomenon and analysis of the system that produced the degraded image

Typical items of information required are estimates or models of the impulse

response of the degrading lter the PSF or equivalently the MTF the PSD

or ACF of the original image and the PSD of the noise If the degrading

system is shiftvariant then a model of the variation of its impulse response

across the eld of imaging would be required The success of a procedure

for image restoration depends upon the accuracy of the model of degradation

used and the accuracy of the functions used to represent the image degrading

phenomena In this chapter we shall explore several techniques for image

restoration under varying conditions of degradation available information

and optimization

Linear Spaceinvariant Restoration Filters

Assuming the degrading phenomenon to be linear and shiftinvariant the

simplest model of image degradation is

gx y hx y fx y x y

Gu v Hu v F u v u v

where f is the original image h is the impulse response of the degrading

LSI system g is the observed degraded image and is additive random

noise that is statistically independent of the imagegenerating process The

functions represented by uppercase letters represent the Fourier transforms

of the imagedomain functions represented by the corresponding lowercase

letters A block diagram of the image degradation system as above is given

in Figure

FIGURE

Image degradation model involving an LSI system and additive noise

The image restoration problem is dened as follows Given g and some

knowledge of h f and nd the best possible estimate of f When the

degrading phenomenon can be represented by an LSI system it is possible to

design LSI lters to restore the image within certain limits A few wellknown

LSI lters for image restoration are described in the following subsections

Inverse ltering

Let us consider the degradation model expressed in matrix form see Sec

tion as

g h f

with no noise being present The restoration problem may be stated as follows

Given g and h estimate f

In order to develop a mathematical statement of the problem let us consider

an approximation

f to f In the leastsquares approach the criterion for

obtaining the optimal solution is stated as follows Minimize the squared error

between the observed response g and the response

g had the input been

The error between g and

g is given by

g g h

The squared error is given as

g h

g g

Now we can state the image restoration problem as an optimization prob

lem Find

f that minimizes

Taking the derivative of the squared error

in Equation with respect to

f we get

g h

Setting this expression to zero we get

f h

This is the leastsquares or pseudoinverse solution If h is square and non

singular we get

f h

If h is circulant or blockcirculant we have h

see Sec

tion Then

f WD

which leads to

F u v

Gu v

Hu v

This operation represents the inverse lter which may be expressed as

u v

Hu v

It is evident that the inverse lter requires knowledge of the MTF of the

degradation process see Sections and for discussions on meth

ods to derive this information

The major drawback of the inverse lter is that it fails if Hu v has zeros

or if h is singular Furthermore if noise is present as in Equation we

get

F u v F u v

u v

Hu v

Problems arise because Hu v is usually a lowpass function whereas u v

is uniformly distributed over the entire spectrum then the amplied noise at

higher frequencies the second component in the equation above overshadows

the restored image

An approach to address the singularity problem associated with the inverse

lter is the use of the singular value decomposition SVD method A

widely used implementation of this approach is an iterative algorithm based on

Bialys theorem to solve the normal equation the algorithm is also known

as the Landweber iterative method McGlamery demonstrated

the application of the inverse lter to restore images blurred by atmospheric

turbulence

In an interesting extension of the inverse lter to compensate for distor

tions or aberrations caused by abnormalities in the human eye Alonso and

Barreto applied a predistortion or precompensation inverse lter to test

images prior to being displayed on a computer monitor The PSF of the af

fected eye was estimated using the wavefront aberration function measured

using a wavefront analyzer In order to overcome the limitations of the in

verse lter a weighting function similar to the parametric Wiener lter see

Section was applied The subjects participating in the study indicated

improved visual acuity in reading predistorted images of testchart letters

than in reading directly displayed test images

Example The original Shapes test image of size pixels is

shown in Figure a along with its logmagnitude spectrum in part b of

the gure The image was blurred via convolution with an isotropic Gaussian

PSF having a radial standard deviation of two pixels The PSF and the

related MTF are shown in parts c and d of the gure respectively The

blurred image is shown in part e of the gure Gaussiandistributed noise of

variance was added to the blurred image after normalizing the image to

the range the degraded noisy image is shown in part f of the gure

The results of application of the inverse lter to the noisefree and noisy

blurred versions of the Shapes image are shown in Figure in both the

space and frequency domains The result of inverse ltering of the noisefree

blurred image for radial frequencies up to the maximum frequency in u v

shown in part a of the gure demonstrates eective deblurring A small

amount of ringing artifact may be observed upon close inspection due to the

removal of frequency components beyond a circular region see the spectrum

in part b of the gure Inverse ltering of the noisy degraded image even

when limited to radial frequencies less than times the maximum frequency

in u v resulted in signicant amplication of noise that led to the complete

loss of the restored image information as shown in part c of the gure Lim

iting the inverse lter to radial frequencies less than times the maximum

frequency in u v prevented noise amplication but also severely curtailed

the restoration process as shown by the result in part e of the gure The

results illustrate a severe practical limitation of the inverse lter See also

Figures and

Power spectrum equalization

Considering the degradation model in Equation the method of power

spectrum equalization PSE takes the following approach Find a linear

transform L so as to obtain an estimate

fx y Lgx y subject to the

constraint

u v

u v that is the PSD of the restored image be

equal to the PSD of the original image Applying the linear transform L to

a b

c d

e f

FIGURE

a Shapes test image size pixels b Logmagnitude spectrum

of the test image c PSF with Gaussian shape radial standard deviation

pixels d MTF related to the PSF in c e Test image blurred with

the PSF in c f Blurred image in e after normalization to and the

addition of Gaussian noise with variance

a b

c d

e f

FIGURE

a Result of inverse ltering the blurred Shapes image in Figure e

Result of inverse ltering the noisy blurred Shapes image in Figure f

using the inverse lter up to the radial frequency equal to c times and

e times the maximum frequency in u v The logmagnitude spectra of

the images in a c and e are shown in b d and f respectively

Equation as well as the constraint mentioned above to the result we get

u v jLu vj

jHu vj

u v

where Lu v represents the MTF of the lter L Rearranging the expression

above we get

PSE

u v jLu vj

u v

jHu vj

u v

jHu vj

A detailed inspection of the equation above indicates the following proper

ties of the PSE lter

The PSE lter requires knowledge of the PSDs of the original image and

noise processes or models thereof

The PSE lter tends toward the inverse lter in magnitude when the

noise PSD tends toward zero This property may be viewed in terms of

the entire noise PSD or at individual frequency samples

The PSE lter performs restoration in spectral magnitude only Phase

correction if required may be applied in a separate step In most

practical cases the degrading PSF and MTF are isotropic and Hu v

has no phase

The gain of the PSE lter is not aected by zeros in Hu v as long

u v is also not zero at the same frequencies In most cases the

noise PSD is nonzero at all frequencies

The gain of the PSE lter reduces to zero wherever the original im

age PSD is zero The noisetosignal PSD ratio in the denominator of

Equation controls the gain of the lter in the presence of noise

Models of the PSDs of the original image and noise processes may be esti

mated from practical measurements or experiments see Section See

Section for a discussion on extending the PSE lter to blind deblurring

Examples of application of the PSE lter are provided in Sections and

The Wiener lter

Wiener lter theory provides for optimal ltering by taking into account the

statistical characteristics of the image and noise processes

The lter characteristics are optimized with reference to a performance cri

terion The output is guaranteed to be the best achievable result under the

conditions imposed and the information provided The Wiener lter is a pow

erful conceptual tool that changed traditional approaches to signal processing

The Wiener lter performs probabilistic stochastic restoration with the

leastsquares error criterion The basic degradation model used is

g h f

where f and are realvalued secondorderstationary random processes that

are statistically independent with known rstorder and secondorder mo

ments Observe that this equation is the matrix form of Equation The

approach taken to estimate the original image is to determine a linear esti

mate

f Lg to f from the given image g where L is the lter to be derived

The criterion used is to minimize the MSE

Expressing the MSE as the trace of the outer product matrix of the error

vector we have

f f

In expanding the expression above we could make use of the following rela

tionships

f f

Lh f f

L f

h f f

h f

Because the trace of a sum of matrices is equal to sum of their traces the E

and Tr operators may be interchanged in order We then obtain the following

expressions and relationships

f f

the autocorrelation matrix of the original image

with the observation that

because f and are statistically independent processes and

f f

is the autocorrelation matrix of the noise process

Now the MSE may be written as

Note

and

because the autocorrelation matrices are

symmetric and the trace of a matrix is equal to the trace of its transpose

At this point the MSE is no longer a function of f g or but depends

only on the statistical characteristics of f and as well as on h and L

In order to derive the optimal lter L we could set the derivative of

with

respect to L to zero

which leads to the optimal Wiener lter function

Note that this solution does not depend upon the inverses of the individual

ACF matrices or h but upon the inverse of their combination The combined

matrix could be expected to be invertible even if the individual matrices are

singular

Considerations in implementation of the Wiener lter Consider

the matrix to be inverted in Equation

Z h

This matrix would be of size N

for N N images making inversion

practically impossible Inversion becomes easier if the matrix can be written

as a product of diagonal and unitary matrices A condition that reduces the

complexity of the problem is that h

and

each be circulant or block

circulant We can now make the following observations

We know from Section that h is blockcirculant for D LSI opera

tions expressed using circular convolution

is a diagonal matrix if is white uncorrelated noise

In most real images the correlation between pixels reduces as the dis

tance spatial separation increases

is then banded and may be

approximated by a blockcirculant matrix

Based upon the observations listed above we can write see Section

h WD

and

where the matrices D are diagonal matrices resulting from the application of

the DFT to the corresponding blockcirculant matrices with the subscripts

indicating the related entity f original image h degrading system or

noise Then we have

Z WD

W D

The Wiener lter is then given by

The optimal MMSE estimate is given by

f L

Interpretation of the Wiener lter With reference to Equation

we can make the following observations that help in interpreting the nature

of the Wiener lter

g is related to Gu v the Fourier transform of the given degraded

image gx y

is related to the PSD

u v of the original image fx y

is related to the PSD

u v of the noise process

is related to the transfer function Hu v of the degrading system

Then the output of the Wiener lter before the nal inverse Fourier trans

form is given by

F u v

u vH

u vGu v

Hu v

u vH

u v

jHu vj

Gu v

jHu vj

Gu v

Hu v

The Wiener lter itself is given by the expression

u v

jHu vj

We can now note the following characteristics of the Wiener lter

In the absence of noise we have

u v and the Wiener lter

reduces to the inverse lter This is also applicable at any frequency

u v where

u v

The gain of the Wiener lter is modulated by the noisetosignal spec

tral ratio

If the SNR is high the Wiener lter is close to the

inverse lter

If the SNR is poor and both the signal and noise PSDs are white

the ratio

is large and could be assumed to be a constant Then the

Wiener lter is close to the matched lter H

u v

If the noisetosignal PSD ratio is not available as a function of u v

it may be set equal to a constant K

SNR

The lter is then known

as the parametric Wiener lter

The Wiener lter is not often singular or illconditioned Wherever

Hu v the output

F u v

If hx y x y then Hu v that is there is no blurring and

the degradation is due to additive noise only Then

F u v

u v

Gu v

which is the original Wiener lter for the degradation model g f

see Section and Rangayyan The Wiener lter as in Equa

tion was proposed by Helstrom

Comparative analysis of the inverse PSE and Wiener lters

When the noise PSD is zero or at all frequencies where the noise PSD

is equal to zero the PSE lter is equivalent to the inverse lter in mag

nitude

When the noise PSD is zero or at all frequencies where the noise PSD

is equal to zero the Wiener lter is equivalent to the inverse lter

The gains of the inverse PSE and Wiener lter are related as

u vj jL

PSE

u vj jL

u vj

The PSE lter is the geometric mean of the inverse and Wiener lters

that is

PSE

u v L

u v

Because jL

PSE

u vj jL

u vj the PSE lter admits more high

frequency components with larger gain than the Wiener lter There

fore the result will be a sharper but more noisy image

The PSE lter does not have a phase component Phase correction if

required may be applied in a separate step

Examples The original Shapes test image and a degraded version of the

image with blurring via convolution using an isotropic Gaussian PSF having

a radial standard deviation of two pixels as well as the addition of Gaussian

distributed noise of variance to the blurred image after normalizing the

image to the range are shown in Figure a and b respectively

See also Figures and In order to design the appropriate Wiener

and PSE lters a model PSD of the image was prepared by computing the

PSD of an image of a square object of size pixels values of the PSD less

than a limit equal to of its maximum were replaced by the limit The

noise PSD was prepared as an array of constant value equal to the known

variance of the noise times N

where N N is the size the PSD array The

results of the Wiener and PSE lters are shown in Figure c and d

respectively It is evident that both lters have removed the blur to some

extent however as expected the result of the PSE lter is noisier than that

of the Wiener lter

Figures and illustrate the application of the inverse Wiener and

PSE lters to degraded versions of the Lenna and Cameraman images re

spectively The Lenna image was blurred using a Gaussianshaped blur func

tion with a radial standard deviation of three pixels and additive noise to

SNR dB The Cameraman image was blurred with a straight line of

length nine pixels simulating blurring due to motion in the horizontal direc

tion and then degraded further with additive noise to SNR dB Inverse

ltering even to limited ranges of frequency has resulted in distorted images

due to the amplication of noise as well as ringing artifacts The Wiener and

PSE lters have removed the blur to a good extent with control over the

noise Further examples of application of the Wiener lter are provided in

Section

Aghdasi et al applied the parametric Wiener lter to deblur mammo

grams after removing noise using the local LMMSE lter see Section or

a Bayesian lter King et al and Honda et al applied the Wiener

lter to the restoration of nuclear medicine images A comparative analysis

of the performance of the PSE Wiener and Metz lters in the restoration of

nuclear medicine images is provided in Section

a b

c d

FIGURE

a Shapes test image size pixels b Test image blurred with

a Gaussian PSF of radial standard deviation pixels and degraded with

additive Gaussian noise with variance after normalization to

c Result of Wiener restoration d Result of restoration using the PSE

lter

a b

c d

e f

FIGURE

a Lenna test image of size pixels b Blurred image with a

Gaussianshaped blur function and noise to SNR dB Deblurred im

ages c Inverse ltering up to the radial frequency times the maximum

frequency in u v d Inverse ltering up to the radial frequency times

the maximum frequency in u v e Wiener ltering f Power spectrum

equalization Reproduced with permission from TF Rabie RM Rangayyan

and RB Paranjape Adaptiveneighborhood image deblurring Journal of

Electronic Imaging

SPIE

a b

c d

e f

FIGURE

a Cameraman test image of size pixels and graylevel range

b Image blurred by pixel horizontal motion and degraded by additive

Gaussian noise to SNR dB Deblurred images c Inverse ltering up

to the radial frequency times the maximum frequency in u v d Inverse

ltering up to the radial frequency times the maximum frequency in u v

e Wiener ltering f Power spectrum equalization Reproduced with per

mission from TF Rabie RM Rangayyan and RB Paranjape Adaptive

neighborhood image deblurring Journal of Electronic Imaging

SPIE

Constrained leastsquares restoration

The Wiener lter is derived using a statistical procedure based on the cor

relation matrices of the image and noise processes The lter is optimal in

an average sense for the class of images represented by the statistical entities

used It is possible that the result provided by the Wiener lter for a specic

image on hand is less than satisfactory Gonzalez and Wintz describe

a restoration procedure that is optimal for the specic image given under

particular constraints that are imposed as follows

Let L be a linear lter operator Using the degradation model as in Equa

tion the restoration problem may be posed as the following constrained

optimization problem Minimize kL

subject to kgh

where

is the estimate of f being derived Using the method of Lagrange multipliers

we now seek

f that minimizes the function

f kL

kg h

where is the Lagrange multiplier Taking the derivative of J

f with respect

f and setting it equal to zero we get

f h

g h

which leads to

f h

h L

where

Due to the illconditioned nature of restoration procedures the results are

often obscured by noise and highfrequency artifacts Artifacts and noise may

be minimized by formulating a criterion of optimality based on smoothness

such as minimizing the Laplacian of the output In order to agree with the

matrixvector formulation as above let us construct a blockcirculant matrix

L using the Laplacian operator in Equation Now L is diagonalized by

the D DFT as D

LW where D

is a diagonal matrix Then we

have

h L

The estimate before the nal inverse Fourier transform is given by

f D

which is equivalent to

F u v

u v

jHu vj

jLu vj

Gu v

where Lu v is the transfer function related to the constraint operator L

The expression in Equation resembles the result of the Wiener lter

in Equation but does not require the PSDs of the image and noise

processes However estimates of the mean and variance of the noise process

are required to determine the optimal value for which must be adjusted to

satisfy the constraint kg h

It is also worth observing that if

the lter reduces to the inverse lter

Gonzalez and Wintz give an iterative procedure to estimate as fol

lows Dene a residual vector as

r g h

g h

h L

We wish to adjust such that

krk

where is a factor of accuracy Iterative trialanderror methods or the

NewtonRaphson procedure may be used to determine the optimal value for

A simple iterative procedure to nd the optimal value for is as follows

Choose an initial value for

Compute

F u v and

Form the residual vector r and compute krk

Increment if krk

or decrement if krk

and return to Step

Stop if krk

Note kk

is the total squared value or total energy of the noise process and

is given by

multiplied with the image area

Other approaches to constrained restoration Several constrained op

timization methods have been developed such as the constrained least squares

method the maximumentropy method and the LeahyGoutis method

The constrained least squares method requires only information about the

mean and variance of the noise function and the estimate of the image is

dened as that which minimizes the output of a linear operator applied to the

image the result may be further subjected to an a priori constraint or an a

posteriori measurable feature such as the smoothness of the image

The maximumentropy approach maximizes the entropy of the image data

subject to the constraint that the solution should reproduce the ideal im

age exactly or within a tolerable uncertainty dened by the observation noise

statistics The LeahyGoutis method optimizes a convex crite

rion function such as minimum energy or cross entropy over the intersection

of two convex constraint sets

In constrained approaches the estimate needs to satisfy a priori constraints

on the ideal image The a priori constraints may be obtained from informa

tion about the blur noise or the image An example of constrained

approaches is the PSE lter which estimates the image using a linear trans

form of the image based on the constraint that the PSD of the ltered image

be equal to that of the ideal image see Section

The maximumaposteriori probability MAP technique is designed to nd

the estimate that maximizes the probability of the actual image conditioned

on the observation posterior probability In the case of linear systems and

under Gaussian assumptions the MAP estimate reduces to the LMMSE es

timate

Several other constraints may be imposed upon restoration lters and their

results A simple yet eective constraint is that of nonnegativity of the im

age values which is relevant in several practical applications where the pixels

represent physical parameters that cannot take negative values Biraud

described techniques to increase the resolving power of signals by imposing

nonnegativity Boas and Kac derived inequalities for the Fourier trans

forms of positive functions that may be used as limits on spectral components

of restored signals or images Webb et al developed a constrained decon

volution lter for application to liver SPECT images several parameters of

the lter could be varied to render the lter equivalent to the inverse Wiener

PSE and other wellknown lters

The Metz lter

Metz and Beck proposed a modication to the inverse lter for applica

tion to nuclear medicine images including noise suppression at high frequen

cies The Metz lter is dened as

u v

Hu v

where is a factor that controls the extent in frequency up to which the in

verse lter is predominant after which the noisesuppression feature becomes

stronger It is apparent that if the gain of the Metz lter reduces to

zero Given that most degradation phenomena have a blurring or smoothing

eect Hu v is a lowpass lter As a consequence H

u v is a highpass

lter and the numerator in Equation is a lowpass lter whose response

is controlled by the factor The factor may be selected so as to minimize

the MSE between the ltered and the ideal images

King et al Gilland et al Honda

et al Hon et al and Boulfelfel et al applied

the Wiener and Metz lters for the restoration of nuclear medicine images A

comparative analysis of the performance of the Metz lter with other lters

in the restoration of nuclear medicine images is provided in Section

Information required for image restoration

It is evident from Equations and that image restoration

using techniques such as the inverse PSE and Wiener lters requires spe

cic items of information including the MTF of the degradation phenomenon

Hu v the PSD of the noise process

u v and the PSD of the original

image process

u v Although this requirement may appear to be onerous

methods and approaches may be devised to obtain each item based upon prior

or additional knowledge measurements and derivations

Several methods to obtain the PSF and MTF are described in Sections

and The PSF or related functions may be measured if the imaging sys

tem is accessible and phantoms may be constructed for imaging as described

in Sections and The PSF may also be modeled mathematically if the

blurring phenomenon is known precisely and can be represented as a mathe

matical process an example of this possibility is described in Section

The PSD of an imagegenerating stochastic process may be estimated as

the average of the PSDs of several observations or realizations of the process

of interest However care needs to be exercised in selecting the samples such

that they characterize the same process or type of images For example the

PSDs of a large collection of highquality CT images of the brain of similar

characteristics may be averaged to derive a PSD model to represent such a

population of images It would be inappropriate to combine CT images of the

brain with CT images of the chest or the abdomen or to combine CT images

with MR images In order to overcome the limitations imposed by noise as

well as a specic and small collection of sample images a practical approach

would be to estimate the average ACF of the images t a model such as a

Laplacian and apply the Fourier transform to obtain the PSD model

The PSD of noise may be estimated by using a collection of images taken

of phantoms with uniform internal characteristics or a collection of parts of

images outside the patient or object imaged representing the background

When the required functions are unavailable it will not be possible to apply

the inverse PSE or Wiener lters However using a few approximations it

becomes possible to overcome the requirement of explicit and distinct knowl

edge of some of the functions image restoration may then be performed while

remaining blind to these entities Methods for blind deblurring are de

scribed in Section

Motion deblurring

Images acquired of living organisms are often aected by blurring due to

motion during the period of exposure imaging In some cases it may be

appropriate to assume that the motion is restricted to within the plane of

the image Furthermore it may also be assumed that the velocity is constant

over the usually short period of exposure Under such conditions it becomes

possible to derive an analytical expression to represent the blurring function

PSF

Consider the imaging of an object fx y over the exposure interval T

Let the relative motion between the object and the camera or the imaging

system be represented as a displacement of t t at the instant of time

t Then the recorded image gx y is given by

gx y

f x t y t dt

Observe that due to the integration over the period of exposure the resulting

image gx y is not a function of time t In order to derive the MTF of the

imaging blurring system we may analyze the Fourier transform of gx y

as follows

Gu v

gx y expj ux vy dx dy

f x t y t dt

expj ux vy

dx dy

Reversing the order of integration with respect to t and x y we get

Gu v

f x t y t expj ux vy dx dy

The expression within the braces above represents the Fourier transform of

f x t y t which is expfj u t v tgF u v Then we

have

Gu v F u v

expfj u t v tg dt

Therefore the MTF of the system is given by

Hu v

expfj u t v tg dt

Thus if the exact nature of the motion during the period of exposure is known

the MTF may be derived in an analytical form

Let us consider the simple situation of linear motion within the plane of the

image and with constant velocity such that the total displacement during the

period of exposure T is given by a b in the x y directions respectively

Then we have t

and t

and

Hu v

exp

sinua vb

ua vb

expj ua vb

In the case of linear motion in one direction only we have a or b

and the D function above reduces to a D function It follows that the

corresponding PSF is a rectangle or box function that reduces to a straight

line in the case of motion in one direction only

Given full knowledge of the MTF of the degradation phenomenon it be

comes possible to design an appropriate restoration lter such as the inverse

PSE or Wiener lter However it should be observed that the sinc function

in Equation has several zeros in the u v plane at which points the

inverse lter would pose problems

An example of simulated motion blurring is given in Section along

with its restoration see also Figure See Sondhi Gonzalez and

Wintz and Gonzalez and Woods for several examples of motion

deblurring

Blind Deblurring

In some cases it may not be possible to obtain distinct models of the degra

dation phenomena An inspection of Equation reveals the fact that the

given degraded image contains information regarding the blurring systems

MTF and the noise spectrum albeit in a combined form in particular we

have

u v jHu vj

u v

Several methods have been proposed to exploit this feature for blind de

blurring with the adjective representing the point that the procedures do

not require explicit models of the degrading phenomena Two procedures for

blind deblurring are described in the following paragraphs

Extension of PSE to blind deblurring The blurred image itself may

be used to derive the parameters required for PSE as follows

Suppose that the given N N image gmn is broken into M M segments

mn l Q

where Q

and M is larger than the dimensions

of the blurring PSF Then

mn hmn f

and

u v jHu vj

u v

In the expressions above the combined eect of blurring across the boundaries

of adjacent subimages is ignored If we now average the PSDs

over all of

the Q

available image segments and make the further assumption that the

averaged PSDs tend toward the true signal and noise PSDs we have

u v jHu vj

u v

where indicates that the corresponding PSDs are estimates The expression

in Equation gives the denominator of the PSE lter as in Equation

the numerator

is required to be estimated separately The MTF Hu v

and the noise PSD

u v need not be estimated individually the procedure

is thus blind to these entities

Iterative blind deblurring

Rabie et al presented an iterative technique for blind deconvolution

under the assumptions that the MTF of the LSI system causing the degrada

tion has zero phase and that its magnitude is a smoothly varying function

of frequency These assumptions are valid for several types of degradation

such as motion blur and outoffocus blur It has been demon

strated in several works that the spectral magnitude in the Fourier repre

sentation of a signal is aected by the blur function while many of the im

portant features of the signal such as edge locations are preserved in the

phase For example it has been shown

that the intelligibility of a sentence is retained if the phase of the Fourier

transform of a long segment of the speech signal is combined with unit mag

nitude

The method of Rabie et al makes use of the image characteristics edge

information preserved in the phase of the blurred image and attempts to re

cover the original magnitude spectrum that was altered by the blur function

Their blind deconvolution algorithm described in the following paragraphs

diers from earlier related work in that the averaging of the

PSD is achieved by smoothing in the frequency domain hence the method is

based upon the use of the entire image rather than sections of the image see

Section This feature relaxes the assumption on the region of support

ROS of the PSF Another key feature of the algorithm is that further en

hancement of the initial estimate of the image is achieved through an iterative

approach

Using the degradation model expressed in Equation and neglecting the

presence of noise we have

u v M

u vM

u v

and

u v

where M

u v is the spectral magnitude of the original image M

u v is

the spectral magnitude of the degraded image M

u v is the degradation

MTF with the property that it is a smooth magnitude function with zero

phase

u v is the spectral phase of the original image and

u v is

the spectral phase of the degraded image The blur model is thus dened

in terms of magnitude and phase The spectral magnitude M

u v of the

original image may be recovered from the spectral magnitude M

u v of

the degraded image as follows The initial estimate of M

u v is based on

smoothing the spectral magnitude of the blurred image M

u v and using

the assumption that M

u v is smooth If we let S denote a D linear

separable smoothing operator then a smoothed M

u v is given by

S M

u v S M

u vM

u v

Because M

u v is a smooth function and S is separable S M

u v M

u v may be approximated by S M

u v M

u v Then we

have

S M

u v S M

u v M

u v

Combining Equation and Equation we obtain

u v M

u v

S M

u v

S M

u v

Equation suggests that if we can obtain an initial approximation of

u v we can rewrite Equation in an iterative form and use it to

rene the initial magnitude estimate Equation can thus be written in

an iterative form as

S M

where l is the iteration number indicates an estimate of the variable under

the symbol and the argument u v has been dropped for compact notation

The initial estimate

may be derived from the phase of the blurred image

that is exp j

u v which retains most of the highfrequency information

spatial edges in the image A simple initial estimate of the original image in

the frequency domain may be dened as the sum of the phase of the blurred

image and the Fourier transform of the blurred image itself

exp j

in terms of magnitudes we have

From Equation it is evident that a unit constant is being added to

the spectral magnitude at all frequencies which would only raise the entire

spectral magnitude response by the corresponding amount This has an eect

comparable to that of adding a highpass ltered version of the image to the

blurred image which would be of amplifying the highfrequency components

in the image This step would however produce a noisy initial approximation

of M

Rather than simply adding a unit magnitude to recover high frequen

cies it would be more appropriate to add those highfrequency components in

that were aected by blurring Although we do not have explicit knowl

edge of M

we do have a ratio between M

and its smoothed version derived

from Equation giving

S M

Adding this ratio to the blurred magnitude spectrum gives

S M

The expression above may be expected to be a more accurate approximation of

than that in Equation because the amount being added is the ratio

S M

which contains more information about the original spectral magnitude

than a simple constant unity in Equation

The advantages of adding

S M

to the blurred magnitude spectrum are the

following

At low frequencies M

S M

Therefore

which

would not signicantly aect the magnitude spectrum

At higher frequencies S M

at high amplitudes of M

because

variations in M

would be averaged out in S M

Thus

S M

at highfrequency components that are likely to represent spatial edges

Adding this ratio to the blurred magnitude spectrum would have the

eect of amplifying highfrequency edges more than the highfrequency

noise components Therefore the highfrequency components of the

phase spectrum exp j

are emphasized in

S M

exp j

for this rea

son Rabie et al referred to the corresponding image as the enhanced

phase image Thus the operation in Equation tends to add to

the normalized highfrequency components that were aected by

blurring

The iterative blind restoration procedure of Rabie et al may be summarized

as follows

Use Equation to obtain an initial estimate of M

Update the estimate iteratively using Equation

Stop when the MSE between the l

estimate and the l

estimate

is less than a certain limit

Using Equation combine the best estimate of M

u v with the

phase function

u v

u v The Fourier transform of the restored

image is given as

F u v

u v exp j

u v

Take the inverse Fourier transform of Equation to obtain the de

blurred image

fx y

Although the method described above neglects noise it can still be used

for deblurring images corrupted by blur and additive noise after rst reducing

the noise in the blurred image

Examples Figure a shows the original test image Lenna Part

b of the gure shows the test image blurred by a Gaussianshaped PSF with

the radial standard deviation

pixels The magnitude of the inverse

Fourier transform of the enhanced phase function

S M

exp j

is shown in

part c of the gure It is evident that the enhanced phase image retains

most of the edges in the original image that were made weak or imperceptible

by blurring Part d of the gure shows the initial estimate of

fx y formed

by the addition of the image in part b of the gure to the image in part

c as described in Equation The addition has emphasized the edges

of the blurred image thus giving a sharper image albeit with a larger MSE

than that of the blurred image as listed in the caption of Figure The

dynamic range of the image in part d is dierent from that of the original

image in part a of the gure

The image generated by combining the rst estimate

obtained by using

Equation with the phase of the blurred image exp j

is shown in part

e of the gure Even after only one iteration the deblurred image closely

resembles the original image The restored image after four iterations shown

in part f of the gure demonstrates sharp features with minimal artifacts

The smoothing operator S was dened as the mean lter applied to the

Fourier spectrum of the image

Figures and illustrate the restoration of an image with text

blurred to two dierent levels The enhanced phase images clearly depict

the edge information retained in the phase of the blurred image The restored

images demonstrate signicantly improved quality in terms of sharpened edges

and improved readability of the text as well as reduced MSE

a b

c d

e f

FIGURE

Iterative blind deconvolution with the Lenna test image of size pixels

and gray levels a original image b blurred image with Gaussian

shaped blur function of radial standard deviation

pixels MSE

c enhanced phase image with the range mapped to the display range

d initial estimate image MSE e deblurred image after the

rst iteration MSE f deblurred image after four iterations MSE

Reproduced with permission from TF Rabie RB Paranjape and RM

Rangayyan Iterative method for blind deconvolution Journal of Electronic

Imaging

SPIE

a b

c d e

FIGURE

Iterative blind deconvolution of a slightly blurred text image digitized to

pixels and gray levels a original image b blurred image with

Gaussianshaped blur function of radial standard deviation

pixels

MSE c enhanced phase image d initial estimate image MSE

e nal restored image after eight iterations MSE Reproduced

with permission from TF Rabie RB Paranjape and RM Rangayyan

Iterative method for blind deconvolution Journal of Electronic Imaging

SPIE

a b

c d e

FIGURE

Iterative blind deconvolution of a severely blurred text image digitized to

pixels and gray levels a original image b blurred version

of the test image with Gaussianshaped blur function of radial standard de

viation

pixels MSE c enhanced phase image d initial

estimate image MSE e restored image after one iteration MSE

Reproduced with permission from TF Rabie RB Paranjape and RM

Rangayyan Iterative method for blind deconvolution Journal of Electronic

Imaging

SPIE

Homomorphic Deconvolution

Consider the case where we have an image that is given by the convolution of

two component images as expressed by the relation

gx y hx y fx y

Similar to the case of the multiplicative homomorphic system described in

Section the goal in homomorphic deconvolution is to convert the convo

lution operation to addition From the convolution property of the Fourier

transform we know that Equation translates to

Gu v Hu v F u v

Thus the application of the Fourier transform converts convolution to multi

plication Now it is readily seen that the multiplicative homomorphic system

may be applied to convert the multiplication above to addition Taking the

complex logarithm of Gu v we have

logGu v logHu v logF u v Hu v F u v u v

fNote logF u v

F u v logjF u vjj

F u vg A linear lter may

now be used to separate the transformed components of f and h with the

assumption as before that they are separable in the transform space A series

of the inverses of the transformations applied initially will take us back to the

original domain

Figure gives a block diagram of the steps involved in a homomorphic

lter for convolved signals Observe that the path formed by the rst three

blocks the top row transforms the convolution operation at the input to ad

dition The set of the last three blocks the bottom row performs the reverse

transformation converting addition to convolution The lter in between thus

deals with transformed images that are combined by simple addition Prac

tical application of the homomorphic lter is not simple Figure gives

a detailed block diagram of the procedure

The complex cepstrum

The formal denition of the complex cepstrum states that it is the inverse

ztransform of the complex logarithm of the ztransform of the input signal

The name cepstrum was derived by transposing the syllables of

the word spectrum other transposed terms are less commonly

used In practice the Fourier transform is used in place of the ztransform

Given gx y hx y fx y it follows that

Gu v

Hu v

F u v

and furthermore that the complex cepstra of the signals are related simply

gx y

hx y

fx y

Here the symbol over a function of frequency indicates the complex loga

rithm of the corresponding function whereas the same symbol over a function

of space indicates its complex cepstrum

An important consideration in the evaluation of the complex logarithm of

Gu v relates to the phase of the signal The phase spectrum computed

as its principal value in the range given by tan

imaginaryfGuvg

realfGuvg

will almost always have discontinuities that will conict with the require

ments of the inverse Fourier transform to follow Thus Gu v needs to

be separated into its magnitude and phase components the logarithmic op

eration applied to the magnitude the phase corrected to be continuous by

adding correction factors of at discontinuities larger than and the

two components combined again before the subsequent inverse transforma

tion Correcting the phase spectrum as above is referred to as phase unwrap

ping In addition it has been shown that a linear

phase term if present in the spectrum of the input may cause rapidly decay

ing oscillations in the complex cepstrum It is advisable to remove the

linear phase term if present during the phaseunwrapping step The linear

phase term may be added to the ltered result if necessary

Taxt applied D homomorphic ltering to improve the quality of med

ical ultrasound images The recorded image was modeled as the convolution

of a eld representing the soundreecting structures the anatomical surfaces

of interest with a PSF related to the interrogating ultrasound pulse shape

A Wiener lter was then used to deconvolve the eects of the PSF

An important application of D homomorphic ltering is in the extraction of

the basic wavelet from a composite signal made up of quasiperiodic repetitions

or echoes of the wavelet such as a voicedspeech signal or seismic signal

An extension of this technique to the removal of visual

echoes in images is described in Section

Echo removal by Radondomain cepstral ltering

Homomorphic deconvolution has been used to remove repeated versions of

a basic pattern in an image known as a visual echo and other

applications in image processing Martins and Ran

gayyan performed D homomorphic deconvolution on the Radon

transform projections of the given image instead of D homomorphic l

tering the nal ltered image was obtained by applying the ART method of

reconstruction from projections see Section to the ltered projections

The general scheme of the method is shown in Figure the details of the

method are described in the following paragraphs

D e c o n v o l u t i o n D e b l u r r i n g a n d R e s t o r a t i o n

FIGURE

Operations involved in a homomorphic lter for convolved signals The symbol at the input or output of each block indicates

the operation that combines the signal components at the corresponding step Reproduced with permission from Reproduced

with permission from RM Rangayyan Biomedical Signal Analysis A CaseStudy Approach IEEE Press and Wiley New

York NY

IEEE

B i o m e d i c a l I m a g e A n a l y s i s

FIGURE

Detailed block diagram of the steps involved in deconvolution of signals using the complex cepstrum Reproduced with

permission from ACG Martins and RM Rangayyan Complex cepstral ltering of images and echo removal in the Radon

domain Pattern Recognition

Pattern Recognition Society Published by Elsevier Science Ltd

FIGURE

Homomorphic complex cepstral ltering of an image in the Radon domain

Reproduced with permission from ACG Martins and RM Rangayyan

Complex cepstral ltering of images and echo removal in the Radon do

main Pattern Recognition

Pattern Recognition

Society Published by Elsevier Science Ltd

Let f

x y be an image given by

x y fx y dx y

where

dx y x y a x x

y y

with a being a scalar weighting factor In the context of images with visual

echoes we could label fx y as the basic element image dx y as a eld

of impulses at the positions of the echoes including the original element

and f

x y as the composite image Applying the Radon transform see

Section we have the projection p

t of the composite image f

x y at

angle given by

x y x cos y sin t dx dy

which leads to

x y

x cos y sin t dx dy d d

x x

y y

x cos y sin t dx dy d d

Here t is the displacement between the projection samples rays in the Radon

domain Using the properties of the function see Section we get

f cos sin t d d

f x

cos y

sin t d d

Dening

cos y

sin

and

f cos sin t d d

we get

t f

t a f

t n

Here f

t represents the projection of the basic element image fx y at

angle and n

is a displacement or shift in the Radon domain for the given

angle and position x

Equation indicates that the projection of

the composite image at a given angle is the summation of the projections at

the same angle of the basic element image and its replication or echo with

an appropriate shift factor n

In order to demonstrate the eect of an echo in a signal on its complex

cepstrum we could apply the Fourier transform the log operation and then

the inverse Fourier transform as follows Applying the Fourier transform to

Equation we get

w F

w a expj w n

where P

w and F

w are the Fourier transforms of p

t and f

t re

spectively and w is the Fourierdomain variable corresponding to the space

domain Radondomain variable t Applying the natural log we get

w ln a expj w n

where the represents the logtransformed version of the variable under the

symbol If a the log function may be expanded into a power series as

w a expj w n

a expj w n

Applying the inverse Fourier transform we get the complex cepstrum of p

t a t n

t n

Thus the complex cepstrum p

t of the projection p

t at angle of the

image f

x y is composed of the complex cepstrum

t of the projection at

angle of the basic element image fx y and a train of impulses If

and the impulse train are su ciently separated in the cepstral domain the

use of a lowpass lter followed by the inverse cepstrum operation gives the

ltered projection f

t of the basic element image The impulse train may

also be suppressed by using a notch or a comb function After the ltered

projections are obtained at several angles it will be possible to reconstruct

the basic element image fx y If the eects of the echoes are considered to

be a type of image distortion ltering as above may be considered to be an

operation for image restoration

If a that is if the echoes are of amplitude equal to or greater than that

of the basic element image the impulses in the cepstrum will appear with

negative delays and the ltering operation will lead to the extraction of an

echo image The situation may be modied easily by applying decaying

exponential weighting factors to the original image and!or the projections

such that the weighted echoes and!or their projections are attenuated to lower

amplitudes Then the eect is equivalent to the situation with a

Example A composite image with ve occurrences of a simple circular

object is shown in Figure a If the object at the lowerright corner of

the image is considered to be the original object the remaining four objects

could be viewed as echoes of the original object Although echoes would in

general be of lower amplitude intensity than the original object they have

been maintained at the same intensity as the original in this image The test

image is of size pixels the radius of each circle is pixels and the

intensity of each circle is on an inverted gray scale

The Radondomain homomorphic lter was applied to the test image The

image was multiplied by a weighting function given by y

where y is the

vertical axis of the image and the origin is at the topleft corner of the image

Furthermore another weighting function given by

with was

applied to each projection before computing the complex cepstrum The

weighting functions were used in order to reduce the eect of the echoes and

facilitate ltering of the cepstrum The cepstrum was lowpass ltered

with a window of pixels Eighteen projections in the range

in steps

were computed ltered and used to reconstruct the basic element

image The result shown in Figure b was thresholded at the gray

level of It is seen that a single circular object has been extracted with

minimal residual artifact

The method was extended to the extraction of the basic element in images

with periodic texture known as the texton by Martins and Rangayyan

see Section

Spacevariant Restoration

Several image restoration techniques such as the Wiener and PSE lters are

based upon the assumption that the image can be modeled by a stationary

random eld Restoration is achieved by ltering the degraded image with

an LSI lter the frequency response of which is a function of the PSD of the

uncorrupted image There are at least two di culties in this approach most

FIGURE

Illustration of homomorphic complex cepstral ltering of an image in the

Radon domain to remove echoes in images a Original image b Filtered

image after thresholding Reproduced with permission from ACG Martins

and RM Rangayyan Complex cepstral ltering of images and echo removal

in the Radon domain Pattern Recognition

Pat

tern Recognition Society Published by Elsevier Science Ltd

images are not stationary and at best may be described as locally stationary

furthermore in practice the PSD of the uncorrupted original image is not

given and needs to be estimated

A common procedure to estimate the PSD of a signal involves sectioning

the signal into smaller presumably stationary segments and averaging their

modied PSDs or periodograms

The procedure assumes shiftinvariance of the blur PSF and averages

out the nonstationary frequency content of the original image In order for

the sectioning approach to be valid the blurring PSF must have an ROS

spatial extent that is much smaller than the size of the subimages as a

result the size of the subimages cannot be made arbitrarily small Thus the

number of subimages that may used to form the ensemble average is limited

consequently the variance of the PSD estimate from the subimages could be

high which leads to a poor estimate of the PSD of the original image

Another consequence of the assumption of image stationarity and the use

of spaceinvariant ltering is the fact that the deblurred images suer from

artifacts at the boundaries of the image In a simple mathematical rep

resentation of this problem it is assumed that the image is of innite extent

however practical images are of nite extent and the performance of a lter

may vary depending upon the assumptions made regarding the edges of the

image The eect at the edges from deconvolution with incomplete informa

tion due to the lack of information beyond the image boundary could cause

dierent contributions from outside the image boundary during deblurring

than those that were convolved into the image during the actual degrada

tion process This leads to a layer of boundary pixels taking incorrect values

during deblurring and consequently artifacts at the image boundaries

Attempts to overcome the problems caused by the inherent nonstationar

ity of images have led to several methods Angel and Jain proposed

a technique to solve the general superposition integral iteratively by using a

conjugategradient method The method faces problems with convergence in

the presence of noise Techniques have been proposed to enhance the perfor

mance of nonadaptive lters by using radiometric and geometric transforms

to generate images that are nearly stationary block stationary in the rst

and second moments The radiometric transform generates stationary

mean and variance whereas the geometric transform provides stationary au

tocorrelation

Adaptive techniques for spacevariant restoration have been proposed based

upon sectioning the given image into smaller subsections and assuming dier

ent stationary models for each section Two approaches

to sectioned deblurring are described in Sections and

The Kalman lter is the most popular method for truly spacevariant l

tering and is described in Section

Sectioned image restoration

In the iterative sectioned MAP restoration technique proposed by Trussell

and Hunt the input image is divided into small P P sections

and the MAP estimate of the uncorrupted section is developed and iterated

upon for renement This procedure is carried out on each section using an

overlapsave technique to reduce edge artifacts Because sectioning an image

presumably causes each individual section to be close to a stationary process

a simpler approach to sectioned deblurring could be based upon the use of

a conventional LSI lter such as the Wiener or PSE lter to deblur each

section individually the deblurred sections may then be combined to form

the nal deblurred image A technique to reduce edge eects between the

sections is to center each subimage in a square region of size comparable to

that of the input image and then pad the region surrounding the centered

subimage with its mean value prior to ltering Each meanpadded region

may also be multiplied in the space domain with a smooth window function

such as the Hamming window

Using the above argument and assuming that each section of size P P

pixels is large compared to the ROS of the blur PSF but small compared

to the actual image dimensions M M each section of the image may be

expressed as the convolution of the PSF with an equivalent section from the

original undegraded image fmn Then we have

mn hmn f

In the frequency domain Equation becomes

u v Hu vF

u v

Now applying the D Hamming window of size M M given by

cos

to each region padded by its mean to the size M M we obtain

mnw

mn hmn f

mnw

mn hmn f

where the

notation represents the corresponding windowed sections The

PSD of each section g

mn may be expressed as

u v jHu vj

u v

The Wiener or PSE lter may then be applied to each section The nal

restored image

fmn is obtained by extracting the individual sections from

the center of the corresponding ltered meanpadded images and placing

them at their original locations

Limitations exist in sectioned restoration due to the fact that the assump

tion of stationarity within the square sections may not be satised sectioning

using regularly spaced square blocks of equal size cannot discriminate between

at and busy areas of any given image Furthermore because of the lim

itations on the section size that it must be large compared to the ROS of the

blur PSF sections cannot be made arbitrarily small Thus the mean value

of a section could be signicantly dierent from its pixel values and conse

quently artifacts could arise at section boundaries To partially solve the

problem of edge artifacts the sections could be overlapped for example by

onehalf the section size in each dimension This technique however

will not reduce the eects at the image boundaries An adaptiveneighborhood

approach to address this limitation is described in Section along with

examples of application

Adaptiveneighborhood deblurring

Rabie et al proposed an adaptiveneighborhood deblurring AND algo

rithm based on the use of adaptiveneighborhood regions determined individ

ually for each pixel in the input image See Section for a description of

adaptiveneighborhood ltering In the AND approach the image is treated

as being made up of a collection of regions features or objects of relatively

uniform gray levels An adaptive neighborhood is determined for each pixel

in the image called the seed pixel when it is being processed being dened

as the set of pixels connected to the seed pixel and having a dierence in

gray level with respect to the seed that is within specied limits of tolerance

The tolerance used in AND is an additive factor as given by Equation

Thus the tolerance determines the maximum allowed deviation in the gray

level from the seed pixel value within each adaptive neighborhood and any

deviation less than the tolerance is considered to be an intrinsic property

of the adaptiveneighborhood region The number of pixels in an adaptive

neighborhood region may be limited by a predetermined number Q however

there are no restrictions on the shape of the adaptiveneighborhood regions

Assuming that each adaptiveneighborhood region grown is large compared

to the ROS of the PSF each such region may be expressed as the convolution

of the PSF with an equivalent adaptiveneighborhood region grown in the

original undegraded image fmn Thus similar to Equation we have

p q hp q f

p q

where mn is the seed pixel location for which the adaptiveneighborhood

region g

p q was grown and p q give the locations of the pixels within

the region It is assumed that regions corresponding to g

p q may be

identied in the original image and noise elds

Next each adaptiveneighborhood region is centered within a rectangular

region of the same size as the input image M M the area surrounding the

region is padded with its mean value in order to reduce edge artifacts and to

enable the use of the D FFT Thus in the frequency domain Equation

becomes

u v Hu vF

u v

Applying the D Hamming window w

p q see Equation to each

meanpadded adaptiveneighborhood region we obtain

p qw

p q hp q f

p q w

p q

p qw

p q

p q hp q f

p q

where

represents the corresponding windowed regions The PSD of the

region g

p q may be expressed as

u v jHu vj

u v

In deriving the AND lter the stationarity of the adaptiveneighborhood

regions grown is taken into account An estimate of the spectrum of the noise

u v within the current adaptiveneighborhood region grown from the

seed pixel at mn is obtained as

u v A

u vG

u v

where A

u v is a frequencydomain magnitudeonly scale factor that

depends on the spectral characteristics of the adaptiveneighborhood region

grown An estimate of F

u v is obtained from Equation by using

the spectral estimate of the noise

u v in place of

u v Then we

have

u v

Hu v

which reduces to

u v

Hu v

u v

The spectral noise estimator A

u v may be derived by requiring the

PSD of the estimated noise

u v to be equal to the original noise PSD

u v for the current adaptiveneighborhood region Thus using Equa

tion we can describe the relationship between the noise PSD and the

image PSD as

u v A

u v

From Equation and Equation the spectral noise estimator A

u v is given by

u v

jHu vj

u v

where

u v

for Gaussian white noise The quantity

in Equation gives the denominator in Equation Therefore no

additional information is required about the PSD of the original undegraded

image

The frequencydomain estimate of the uncorrupted adaptiveneighborhood

region is obtained by using the value of A

u v computed from Equa

tion in Equation The spectral estimate of the original unde

graded adaptiveneighborhood region is thus given by

u v

Hu v

u v

jHu vj

u v

Hu v

u v

The spacedomain estimate of the uncorrupted adaptiveneighborhood re

gion

p q is obtained from the inverse Fourier transform of the expres

sion above By replacing the seed pixel at mn with the deblurred pixel

mn and running the algorithm above for every pixel in the input im

age we will obtain a deblurred image based on stationary adaptive regions

A computational disadvantage of the algorithm above is the fact that it

requires two M M D FFT operations per pixel An approach to cir

cumvent this di culty to some extent is provided by the intrinsic nature of

the adaptive neighborhood Because most of the pixels inside an adaptive

neighborhood region will have similar adaptiveneighborhood regions when

they become seed pixels because they lie within similar limits of tolerance

instead of growing an adaptiveneighborhood region for each pixel in the in

put image we could grow adaptiveneighborhood regions only from those

pixels that do not already belong to a previously grown region Thus after

ltering an adaptiveneighborhood region the entire adaptiveneighborhood

region may be placed in the output image at the seed pixel location instead of

replacing only the single restored seed pixel Note that the various adaptive

neighborhood regions grown as above could still overlap This approach is a

compromise to reduce the computational requirements

Examples The test image Lenna of size pixels and

gray levels is shown in Figure a The image after degradation by

a Gaussianshaped blur PSF with a radial standard deviation

pixels

and additive white Gaussian noise to dB SNR is shown in part b of the

gure Parts c e of the gure show three dierent xedneighborhood

sections of the blurred image Each section is of size pixels and is

centered in a square region of the same size as the full image

with the surrounding area padded with the mean value of the section Part

f shows the windowed version of the region in part e after weighting with

a Hamming window It is evident from the images in parts c e that the

assumption of stationarity within a given section is not satised each section

contains a variety of image characteristics The values of the meanpadded

areas are also signicantly dierent from the pixel values of the corresponding

centered sections

Three dierent adaptiveneighborhood regions of the blurred test image

are shown in Figure c e Each adaptiveneighborhood region was

allowed to grow to any size as long as the pixel values were within an adaptive

tolerance given by

where

is an estimate of the standard deviation of

the noisefree blurred image gmn hmn fmn Each adaptive

neighborhood region was centered in a square region of the same size as the

full image and the surrounding area was padded with the mean

value of the region Unlike the xed square sections shown in Figure

c e the adaptiveneighborhood regions in Figure c e do

not contain large spatial uctuations such as highvariance edges but rather

slowvarying and relatively smooth details Thus we may assume that each

adaptiveneighborhood region approximates a stationary region in the input

image From Figure c e it should also be observed that the

adaptiveneighborhood regions overlap

Two restored versions of the blurred Lenna image obtained using the sec

tioned Wiener lter with sections of size and pixels with the

adjacent sections overlapped by onehalf of the section size in both dimen

sions are shown in Figure c and d respectively Overlapping the

sections has eectively suppressed the artifacts associated with the inner sec

tions of the image however it does not reduce the artifacts at the boundaries

of the image because of the lack of information beyond the image boundaries

Figure e shows the test image deconvolved using the full image frame

with the Wiener lter In this result less noise smoothing has been achieved

as compared to the restored images using sections of size or

pixels This could be due to the nonstationarity of the fullframe image used

to calculate the PSD required in the Wiener lter equation

The restored image obtained using the AND lter of Equation is

shown in Figure f The result has almost no edge artifacts which is

due to the overlapping adaptiveneighborhood regions used The AND result

has the lowest MSE of all of the results shown as listed in Table

The Cameraman test image of size pixels and gray levels is

shown in Figure a Part b of the gure shows the image after degra

dation by a PSF representing pixel horizontal motion blurring and additive

noise to SNR dB Two restored images obtained using the sectioned

Wiener lter with sections of size and pixels with the adja

cent sections overlapped by onehalf of the section size in both dimensions

are shown in Figure c and d respectively Edge artifacts are appar

ent in these deblurred images the artifacts are more pronounced around the

vertical edges in the image than around the horizontal edges This is due to

the shape of the blur PSF which is a D function along the horizontal axis

The result of deconvolution using the full image frame and the Wiener

lter is shown in Figure e It is evident that edge artifacts do not

exist within the boundaries of this result

The test image restored by the application of the AND lter of Equa

tion is shown in Figure f Almost no edge artifact is present

in this image due to the use of adaptiveneighborhood regions The image is

sharper and cleaner than the other results shown in Figure the AND

result also has the lowest MSE as listed in Table

The results demonstrate that image stationarity plays an important role in

the overall performance of restoration lters The use of adaptiveneighbor

hood regions can improve the performance of restoration lters

The Kalman lter

The Kalman lter is a popular approach to characterize dynamic systems in

terms of statespace concepts The Kalman lter for

mulation could be used for ltering restoration prediction and interpolation

smoothing

Formulation of the Kalman lter In the Kalman lter the signals or

items of information involved are represented as a state vector fn and an

observation vector gn where n refers to the instant of time or is an index

a b

c d

e f

FIGURE

Sectioning of the Lenna image of size pixels and graylevel range of

a Original image b Blurred image with a Gaussianshaped blur

function and noise to SNR dB MSE c d and e Three

sections meanpadded to pixels f Hammingwindowed version

of the region in e Reproduced with permission from TF Rabie RM

Rangayyan and RB Paranjape Adaptiveneighborhood image deblurring

Journal of Electronic Imaging

SPIE

a b

c d

e f

FIGURE

Adaptiveneighborhood segmentation of the Lenna image of size

pixels and graylevel range of a Original image b Blurred image

with a Gaussianshaped blur function and noise to SNR dB MSE

c d and e Three adaptiveneighborhood meanpadded regions

f Hammingwindowed version of the region in e Reproduced with per

mission from TF Rabie RM Rangayyan and RB Paranjape Adaptive

neighborhood image deblurring Journal of Electronic Imaging

SPIE

a b

c d

e f

FIGURE

a Lenna test image b Blurred image with a Gaussianshaped blur func

tion and noise to SNR dB MSE Sectioned deblurring with

overlapped sections of size c pixels MSE and d

pixels MSE e Fullframe Wiener ltering MSE f Adaptive

neighborhood deblurring MSE Reproduced with permission from TF

Rabie RM Rangayyan and RB Paranjape Adaptiveneighborhood image

deblurring Journal of Electronic Imaging

SPIE

a b

c d

e f

FIGURE

a Cameraman test image of size pixels and graylevel range

b Image blurred by pixel horizontal motion and degraded by additive

Gaussian noise to SNR dB MSE Deblurred images c Sec

tioned deblurring with overlapped sections of size pixels MSE

d Sectioned deblurring with overlapped sections of size pixels MSE

e Fullframe Wiener ltering MSE f Adaptiveneighborhood

deblurring MSE Reproduced with permission from TF Rabie RM

Rangayyan and RB Paranjape Adaptiveneighborhood image deblurring

Journal of Electronic Imaging

SPIE

TABLE

Meansquared Errors of the Results of Sectioned and

Adaptiveneighborhood Deblurring of the Lenna and Cameraman

Images of Size Pixels and Gray Levels for Various

Neighborhood Sizes and Two Dierent Blurring Functions and

Approximate Computer Processing Time Using a SUN!Sparc

Workstation

Filter Section Time Meansquared error

size minutes Lenna Cameraman

Degraded " "

Wiener

PSE

Wiener

PSE

Wiener

PSE

Wiener Full frame

PSE Full frame

Adaptive

neighborhood max

Reproduced with permission from TF Rabie RM Rangayyan and RB

Paranjape Adaptiveneighborhood image deblurring Journal of Electronic

Imaging

SPIE

related to the sequence of the state or observation vectors It is assumed

that the input is generated by a driving or process noise source

n and

that the output is aected by an observation noise source

n A state

transition matrix an n is used to indicate the modication of the state

vector from one instant n to the next instant n Note The argument

in the notation an n is used to indicate the dependence of the variable

matrix a upon the instants of observation n and n and not the indices

of the matrix a This notation also represents the memory of the associated

system An observation matrix hn is used to represent the mapping from

the state vector to the observation vector

Figure illustrates the concepts described above in a schematic form A

state vector could represent a series of the values of a D signal or a collection

of the pixels of an image in a local ROS related to the pixel being processed

see Figure for an illustration of two commonly used types of ROS in

image ltering The state vector should be composed of the minimal amount

of data that would be adequate to describe the behavior of the system As we

have seen in Section images may be represented using vectors and image

ltering or transformation operations may be expressed as multiplication with

matrices The state transition matrix an n could represent a linear

prediction or autoregressive model see Section that characterizes the

input state The state vector would then be composed of a series of the

signal or image samples that would be related to the order of the model and

be adequate to predict the subsequent values of the state and observation

vectors with minimal error The observation matrix hn could represent

a blurring PSF that degrades a given image Observe that the Kalman lter

formulation permits the representation of the signals and their statistics as

well as the operators aecting the signals as functions that are nonstationary

dynamic or varying with time or space

FIGURE

Statespace representation of the basic Kalman lter formulation

The Kalman lter formulation represents a new or updated value of

the state vector f recursively in terms of its previous value and new input as

fn an n fn

This is known as the process equation the corresponding model is also known

as the plant process or message model If we let the state vector f be of

size M then the state transition matrix a is of size M M and the

driving noise vector

is of size M We may interpret the driving noise

vector

as the source of excitation of the process system represented by the

state vector f and the state transition matrix a It is assumed that the noise

process

is a zeromean whitenoise process that is statistically independent

of the stochastic process underlying the state vector f The noise process

is characterized by its M M correlation matrix ACF

n if n k

otherwise

The state transition matrix an n is characterized by the following

properties

alm amn al n product rule

mn anm inverse rule

amm I identity

For a stationary system the state transition matrix would be a constant

matrix a that is stationary or independent of time or space

The output side of the dynamic system see Figure is characterized

by a measurement or observation matrix hn that transforms the state vector

fn The result is corrupted by a source of measurement or observation noise

which is assumed to be a zeromean whitenoise process that is statistically

independent of the processes related to the state vector f and the driving noise

It is assumed that the observation vector gn is of size N dierent

from the size of the state vector f which is M then the observation

matrix hn is required to be of size N M and the observation noise

is required to be of size N We then have the measurement or observation

equation

gn hn fn

Similar to the characterization of the driving noise in Equation the

observation noise

n is characterized by itsNN correlation matrix ACF

Due to the assumption of statistical independence of

and

we have

n k

With the situation formulated as above the Kalman ltering problem may

be stated as follows given a series of the observations G

fg g

FIGURE

Regions of support ROS for image ltering a Illustration of the past

present and future in ltering of an image pixel by pixel in rasterscan

order b Nonsymmetric half plane NSHP of order or size P

c Quarter plane QP of order or size P

gng for each n nd the MMSE estimate of the state vector fl

The application is referred to as ltering if l n prediction if l n and

smoothing or interpolation if l n Given our application of interest in

the area of image restoration we would be concerned with ltering Note

The derivation of the Kalman lter presented in the following paragraphs

closely follows those of Haykin and Sage and Melsa

The innovation process An established approach to obtain the solution

to the Kalman ltering problem is via a recursive estimation procedure using

a onestep prediction process and the resultant dierence referred to as the

innovation process Suppose that based upon the set of obser

vations G

fg g gn g the MMSE estimate of fn

has been obtained let this estimate be denoted as

fn jG

Given a

new observation gn we could update the previous estimate and obtain a

new state vector

fnjG

Because the state vector fn and the observation

gn are related via the observation system we may transfer the estimation

procedure to the observation variable and let

gnjG

denote the MMSE

estimate of gn given G

Then the innovation process is dened as

n gn

gnjG

The innovation process n represents the new information contained in gn

that cannot be estimated from G

Using the observation equation

we get

gnjG

fnjG

njG

fnjG

noting that

njG

because the observation noise is orthogonal to

the past observations Combining Equations and we have

n gn hn

fnjG

Using Equation Equation becomes

n hn

n n

where

n n is the predicted state error vector at n using the information

available up to n given by

n n fn

fnjG

The innovation process has the following properties

n is orthogonal to the past observations g g gn and

hence

m n

The innovation process is a series of vectors composed of random vari

ables that are mutually orthogonal and hence

m n

A onetoone correspondence exists between the observations fg

g gng and the vectors of the innovation process f

ng that is one of the two series may be derived from the other

without any loss of information via linear operations

The ACF matrix of the innovation process n is given by

n E

n n h

where

n n E

n n

is the ACF matrix of the predicted state error and the property that

n n

and

n are mutually orthogonal has been used The ACF matrix of the

predicted state error provides a statistical representation of the error in the

predicted state vector

fnjG

Estimation of the state vector using the innovation process The

aim of the estimation process is to derive the MMSE estimate of the state

vector fl Given that g is related to f via a linear transform and that is

linearly related to g we may formulate the state vector as a linear transform

of the innovation process as

fljG

k k

where L

k l n is a series of transformation matrices Now the

predicted state error vector is orthogonal to the innovation process

l n

ffl

fljG

m n

Using Equations and we get

Consequently we get

m E

Using the expression above for L

m in Equation we have

fljG

k k

from which it follows that

fljG

k k

n n

For l n we obtain

fn jG

k k

n n

Using the process equation we have for k n

fan n fn

an nE

In arriving at the expression above the property that

n and k are

mutually orthogonal for k n has been used Now using Equation

and Equation with l n we get

k k

an n

k k

an n

fnjG

Interpretation of the expression in Equation is made easier by the

following formulations

The Kalman gain Let

Kn E

this is a matrix of size M N whose signicance will be apparent after a few

steps The expectation in the equation above represents the crosscorrelation

matrix between the state vector fn and the innovation process n

Using Equations and Equation may be simplied to

fn jG

an n

fnjG

Kn n

This is an important result indicating that we may obtain the MMSE esti

mate of the state vector

fn jG

by applying the state transition matrix

an n to the previous estimate of the state vector

fnjG

and adding a

correction term The correction term Kn n includes the innovation pro

cess n multiplied with the matrix Kn for this reason and in recognition

of the original developer of the underlying procedures the matrix Kn is

referred to as the Kalman gain

In order to facilitate practical implementation of the steps required to com

pute the Kalman gain matrix we need to examine a few related entities as

follows Using Equation we can write

n n

fnjG

n n

where the last step follows from the property that the estimated state vec

tor

fnjG

and the predicted state error vector

n n are mutually

orthogonal Using Equations and we have

E fn

an nE

fn fhn

n n

an nE

n n

an nE

n n

an n

n n h

In arriving at the result above Equation has been used use has also

been made of the property that f and

are independent processes Using

Equation in Equation we get the Kalman gain matrix as

Kn an n

n n h

which upon the use of the expression in Equation for

n gives

Kn an n

n nh

n n h

This expression for the Kalman gain matrix may be used with Equation

to update the estimate of the state vector

Practical computation of the Kalman gain matrix may be facilitated fur

ther by the following derivations Extending Equation one step

further we have

n n fn

fn jG

Putting together Equations and we have

n n an n

fnjG

gn hn

fnjG

Using the measurement equation and Equation the equation

above may be modied to

n n an n

n n

hn fn

n hn

fnjG

an n

n n

Knhn

fnjG

nKn

an nKnhn

n n

nKn

Putting together Equations and and noting the property that

the processes

and

are mutually uncorrelated we have the ACF

matrix of the predicted state error

n n as

n n E

n n

an nKnhn

n n

an nKnhn

n Kn

an n

n n a

n n

Knhn

n n a

n n

an n

n n h

Knhn

n n h

n Kn

an n

n n a

n n

Knhn

n n a

n n

an n

n n h

an n

n n a

n n

Knhn

n n a

n n

an n

n n h

which results in

n n an n

n a

n n

In arriving at the result above Equations and have been used

A new matrix

n of size M M has been introduced dened as

n n an n Knhn

n n

use has been made of the property a

n n an n which follows

from the inverse rule in Equation Equation is known as the

Riccati equation and assists in the recursive computation of the ACF matrix

of the predicted state error

The procedures developed to this point are referred to as Kalmans onestep

or onestage prediction algorithm The algorithm is represented by

in order Equations and

Application to ltering In ltering the aim is compute the estimate

fnjG

The onestep prediction algorithm developed in the preceding para

graphs may be extended to the ltering application as follows

Because the processes f and

are mutually independent it follows from

Equation that the MMSE estimate of fn given G

fn jG

an n

fnjG

njG

an n

fnjG

Premultiplying both sides by an n and using the inverse rule in Equa

tion we get

an n

fn jG

an n an n

fnjG

an n

fn jG

Thus given the result of the onestep prediction algorithm fn jG

can derive the ltered estimate

fnjG

Let us now consider the ltered estimation error dened as

n gn hn

fnjG

Using Equations and we can modify Equation

as follows

n gn hn

an n

fn jG

gn hn

an n fan n

fnjG

Kn ng

gn hn

fnjG

hn an n Kn n

n hn an n Kn n

I hn an n Kn n

This expression indicates that the ltered estimation error

n is related to

the innovation process n through a conversion factor that is given by the

matrix within the square brackets Using Equations and we

may simplify the expression above as follows

I hn an n an n

n n h

I hn

n n h

n hn

n n h

The dierence between the true state vector fn and the ltered estimate

fnjG

labeled as the ltered state error vector

n is given by

n fn

fnjG

Using Equations and we may modify the equation above as

n fn an n

fn jG

fn an n

an n

fnjG

Kn n

fn an n an n

fnjG

an n Kn n

fnjG

an n Kn n

n n an n Kn n

Equation has been used in the last step above where

n n is the

predicted state error vector at n using the information provided up to n

The ACF matrix of the ltered state error vector

n is obtained as

follows

n n

an n KnE

n a

n n

n a

n n

an n KnE

n n

The expectation in the third term in the expression above may be simplied

n n

ffn

fnjG

because

fnjG

is orthogonal to n Using Equation with k n

and premultiplying both sides with a

n n an n we get

an n E

an n Kn

Equation has been used for the second step above Therefore we have

n n

an n Kn

Using a similar procedure the expectation in the fourth term of Equation

may be modied as

n n

n a

n n

Substituting Equations and in Equation we get

n n an n Kn

n a

n n

From Equation we get

n an n

n n h

using which we can modify Equation as

n n an n an n

n n

n a

n n

n n h

n a

n n

The inverse rule of the state transition matrix see Equation has been

used in the last step above Because ACF matrices are symmetric we may

transpose the expression in Equation and obtain

n n an n Knhn

n n

I an n Knhn

n n

where Equation has been used for the last step This result indicates

that the matrix

n introduced in the Riccati equation is the ACF

matrix of the ltered state error

Initial conditions In practice the initial state of the process equation

will not be known However it may be possible to describe it in a

statistical manner in terms of the mean and ACF of the state vector or

estimates thereof The initial conditions given by

fjG

and

f f

result in an unbiased ltered estimate

Summary of the Kalman lter The following description of the Kalman

lter based upon onestep prediction summarizes the main principles and

procedures involved and assists in implementing the lter

Data available The observation vectors G

fg g gng

System parameters assumed to be known

The state transition matrix an n

The observation system matrix hn

The ACF matrix of the driving noise

The ACF matrix of the observation noise

Initial conditions

fjG

a diagonal matrix with values of the order of

Recursive computational steps For n do the following

Using Equation compute the Kalman gain matrix as

Kn an n

n n h

Obtain the innovation process vector using Equation as

n gn hn

fnjG

Using Equation update the estimate of the state vector as

fn jG

an n

fnjG

Kn n

Compute the ACF matrix of the ltered state error given by Equa

tion as

n n an n Knhn

n n

Using Equation update the ACF matrix of the predicted state

error as

n n an n

n a

n n

The Kalman lter formulation is a general formulation of broad scope rel

evance and application Haykin provides a discussion on the Kalman

lter as the unifying basis for recursive leastsquares RLS lters Sage and

Melsa present the derivation of a stationary Kalman lter and relate it

to the Wiener lter

Extension to of the Kalman lter to image restoration Extending

the Kalman lter to D image ltering poses the following challenges

dening the state vector

determining the size spatial extent of the state vector

deriving the dynamic spacevariant state transition matrix process

phenomenon or model

obtaining the dynamic spacevariant observation matrix system

estimating the driving and observation noise correlation matrices and

dealing with the matrices of large size due to the large number of ele

ments in the state vector

Woods and Radewan and Woods and Ingle proposed a scalar

processor that they called as the reducedupdate Kalman lter RUKF dis

cussing in detail the options for the ROS of the lter in particular the NSHP

ROS see Figure b as well as the problems associated with the bound

ary conditions see also Tekalp et al Boulfelfel et al

modied the RUKF algorithm of Woods and Ingle for application to the

restoration of SPECT images with the following main characteristics

The state transition matrix was derived using a D AR model see Sec

tion

The observation matrix was composed by selecting one of PSFs de

pending upon the distance of the pixel being processed from the center

of the image see Figure and Section for details

Filtering operations were performed in a QP ROS as illustrated in Fig

ure c The image was processed pixelbypixel in rasterscan

order by moving the QP ROS window from one pixel to the next

The size of the ROS was determined as the larger of the AR model order

related to the state transition matrix and the width of the PSF

With the assumption that the observation noise follows the Poisson

PDF the variance of the observation noise was estimated as the total

photon count in a local window

FIGURE

In order to apply the nonstationary Kalman lter to SPECT images the

image plane is divided into zones based upon the distance from the center

of the axis of rotation of the camera A dierent PSF MTF is used to

process the pixels in each zone The width of the PSF increases with distance

from the center

The steps and equations of the scalar RUKF algorithm are as follows

Update index n to n at the end of a row reset n and update m

to m

Project the previous estimate of the state vector one step forward by

using the dynamic system model as

m n

Project the error ACF matrix one step forward as

r s

m r n s

and

mnmn

mnm n

Compute the updated Kalman gain matrix as

p q

m n m p n q

p q

m n m p n q

Update the state vector as

p q

p q K

m p n q

gmn

m n

Update the error ACF matrix as

p q

p q K

m p n q

r s

m r n s

In the notation used above the superscript mn indicates the ltering

step position of the QP ltering window the indices within the argument of

a variable indicate the spatial coordinates of the pixels of the variable being

used or processed the subscript b indicates the corresponding variable before

updating the subscript a indicates the corresponding variable after updating

and all summations are over the chosen ROS for the ltering operations The

subscript

has been dropped from the error ACF

With reference to the basic Kalman algorithm documented on page

the RUKF procedure described above has the following dierences

The sequence of operations in the RUKF algorithm above follows the

Kalman lter algorithm described by Sage and Melsa p and

is dierent from that given by Haykin as well as the algorithm on

page

Equation computes only the rst part of the righthand side of

Equation

Equations and together perform the operations repre

sented by Equation

Equation is equivalent to Equation except for the presence

of the state transition matrix an n in the latter

Equation is similar to Equation with the observation

that the innovation process n is given by the expression within the

brackets on the righthand side of the latter

Equation is equivalent to Equation except for the presence

of the state transition matrix an n in the latter Observe that the

term an n in Equation is cancelled by the term an n

in Equation due to the inverse rule given in Equation

Illustrations of application of the D Kalman lter the RUKF algorithm

to the restoration of SPECT images are provided in Section

Application Restoration of Nuclear Medicine Im

ages

Nuclear medicine images including planar SPECT and PET images are

useful in functional imaging of several organs such as the heart brain thyroid

and liver see Section for an introduction to the principles of nuclear

medicine imaging However nuclear medicine images are severely aected

by several factors that degrade their quality and resolution Some of the

important causes of image degradation in nuclear medicine imaging are briey

described below along with suitable methods for image restoration

Poor quality control SPECT images could be blurred by misalign

ment of the axis of rotation of the system in the image reconstruction

process Data in the projection images could become inaccurate due

to defects and nonuniformities in the detector Patient motion during

image data acquisition leads to blurring

Poor statistics The number of photons acquired in a nuclear medicine

image will be limited due to the constrained dose of the radiopharma

ceutical administered the time that the patient can remain immobile

for the imaging procedure the time over which the distribution and ac

tivity of the radiopharmaceutical within the patient will remain stable

the low e ciency of detection of photons imposed by the collimator

the limited photoncapture e ciency of the detection system and the

attenuation of the photons within the body These factors lead to low

SNR in the images SPECT images have poorer statistics due to the

limited time of acquisition of each projection image See Figure

for an illustration of increased levels of noise in planar images due to

limited imaging periods

Photoncounting Poisson noise The emission of gammaray pho

tons is an inherently random process that follows the Poisson distribu

tion see Section When the photon count is large the Poisson

PDF tends toward the Gaussian PDF see Figure The detection of

gammaray photons is also a random process governed by the probabil

ities of interaction of photons with matter The noise in SPECT images

is further aected by the lters used in the process of image reconstruc

tion from projections the noise is amplied and modied to result in a

mottled appearance of relatively uniform regions in SPECT images

Gammaray attenuation Photons are continuously lost as a gamma

ray passes through an object due to scattering and photoelectric ab

sorption p The eect of these phenomena is an overall at

tenuation of the gamma ray The number of photons that arrive at the

detector will depend upon the attenuating eect of the tissues along

the path of the ray and hence will not be a true representation of the

strength at the source the organ that collected the radiopharmaceutical

in relatively higher proportion

Compton scattering Compton scattering occurs when a gamma

ray photon collides with an outershell electron in the medium through

which it passes the scattered or secondary gammaray photon contin

ues in a dierent direction at a lower energy and the electron that gained

the energy is released and scattered pp Gamma

ray photons aected by Compton scattering cause counts in images at

locations that correspond to no true emitting source and appear as

background noise However because the aected photon arrives at the

detector with a lower energy than that at which it was emitted at the

source which is known it may be rejected based upon this knowledge

Poor spatial resolution The spatial resolution of a gamma camera is

aected by two dierent factors the intrinsic resolution of the detector

and the geometric resolution of the collimator The intrinsic resolution

of the detector is related to the diameter of the PMTs the thickness

of the crystal and the energy of the gamma rays used The intrinsic

resolution may be improved by using more PMT tubes a thinner crystal

and gamma rays of higher energy However a thinner crystal would lead

to lower e ciency in the detection of gammaray photons

The geometric resolution of the collimator is a function of the depth

thickness of the collimator the diameter of the holes and the source

tocollimator distance The geometric resolution of a collimator may

be improved by making it thicker and reducing the size of the holes

however these measures would reduce the number of photons that reach

the crystal and hence reduce the overall e ciency of photon detection

Digital image processing techniques may be applied to correct for some of

the degrading phenomena mentioned above lters may also be designed to

remove some of the eects The following sections provide details of

a few methods to improve the quality of nuclear medicine images a schematic

representation of the application of such techniques is shown in Figure

FIGURE

Schematic representation of the application of several techniques to improve

the quality of nuclear medicine images Attenuation correction may be

achieved via averaging of conjugate projections modications to the recon

struction algorithm or other methods applied to the planar projections or the

reconstructed SPECT data The blocks are shown separately to emphasize

the various procedures for quality improvement Only one of the two blocks

shown in dashed lines " prereconstruction or postreconstruction processing

" would be used

Quality control

The quality of images obtained using a gamma camera is aected by imperfec

tions in the detection system of the camera System misalignment detector

nonlinearity and detector nonuniformity are the major contributors to im

age degradation in this category System misalignment errors arise when an

incorrect axis of rotation is used in the reconstruction process or due to non

linearities in the eld of view of the camera and result in smearing of the

SPECT image System misalignment may be corrected through reg

ular calibration of the camera system Detector nonlinearity arises

due to the compression of events located near the centers of the PMTs and an

expansion between the PMTs Detector nonlinearities are usually

not perceptible however they can have signicant eects on image nonuni

formities and misalignment The nite number of PMTs in the detector is

also a source of image nonuniformity which results in variations in counts in

the acquired image of a uniform object

The common approach to correct nonuniformity is by acquiring an image

of a oodeld a uniform source and then using it as a correction matrix for

other images However this method only corrects for variations

in amplitude A more sophisticated method stores correction matrices for

regional dierences in pulseheight spectra and for positions over the entire

area of the detector The correction matrices are then used on an event

byevent basis to compensate for regional dierences by adjusting either the

system gain or the pulseheight window and to compensate for nonlinearities

by accurately repositioning each event

Scatter compensation

Compton scattering results in a broad symmetric distribution centered about

the primary wavelength or energy level of the gamma rays at the source

One approach to reduce the eect of scatter is by energy discrimination in

the camera system by rejecting photons at all energy levels outside a nar

row window centered at the photopeak of the radioisotope most of the

Comptonscattered photons will be rejected however this rejection is not

complete

Jaszczak et al acquired projection images in two dierent energy win

dows one placed at the photopeak of the radioisotope and the other over the

Compton portion of the energy spectrum A fraction of the second image was

subtracted from the photopeak image The procedure was shown to elimi

nate most of the scatter however di culty was encountered in determining

the fraction of the scatter image to be subtracted Egbert and May

modeled Compton scattering using the integral transport equation Correc

tion was performed in an iterative procedure using an attenuationcorrected

reconstructed image as the initial estimate The next image estimate was

computed as a subtraction of the product of the Chang pointwise attenua

tion correction operator with a scattering operator determined from a

given energy threshold from the previous estimate of the image This tech

nique suers from the necessity of having to determine the scattering operator

for each distribution of the scattering medium photon energy and threshold

Axelsson et al proposed a method in which the scatter distribution is

modeled as the convolution of the projection image with an exponential ker

nel and correction is performed by subtracting the estimate of the scatter

from the acquired projection image A similar technique was described by

Floyd et al in which the correction is performed by deconvolution

rather than by subtraction

See Rosenthal et al Buvat et al and Ljungberg et al for

other methods for scatter compensation

Attenuation correction

When gammaray photons pass through body tissues several photons get

attenuated and scattered The amount of attenuation depends upon the at

tenuation coe cient of the tissues and the depth of the source of the photons

The attenuation in nuclear medicine imaging may be represented as

expx

where I

is the intensity of the gamma ray at the source I

is the intensity

at the detector is the attenuation coe cient of the attenuating medium

or tissue assumed to be uniform in the expression above and x is the dis

tance traveled by the gamma ray through the attenuating medium Sev

eral methods for attenuation correction have been described in the literature

the methods may be divided into three categories preprocessing correction

methods intrinsic correction methods and postprocessing correction meth

ods

The most common preprocessing correction procedure is to estimate the

sourcetocollimator distance at each point in the image and use Equation

to compute the gammaray intensity at the source Other ap

proaches include using the arithmetic mean or the geometric mean of con

jugate projections to correct projections for attenuation before reconstruc

tion see Section These two methods are employed in most SPECT

systems and give acceptable results for head and liver SPECT images where

the attenuation coe cients may be assumed to be constant

In the intrinsic correction methods attenuation correction is incorporated

into the reconstruction algorithm Gullberg and Budinger proposed a

technique where the solution to the attenuated Radon integral for the case of

a uniformly attenuating medium is used This technique is limited to images

with high statistics because it can amplify noise and introduce artifacts in

lowcount images Censor et al proposed a correction method where

both the attenuation coe cient map and the radioisotope distribution are

estimated using discrete estimation theory this approach provides noisier

images but allows correction of nonuniform attenuation

A commonly used postprocessing correction method is the Chang iterative

technique in which each pixel of the SPECT image is multiplied by

the reciprocal of the average attenuation along all rays from the pixel to the

boundaries of the attenuating medium The results are then iteratively rened

by reprojecting the image using the assumed attenuation coe cients

The dierence between the true acquired projections and those obtained by

reprojection is used to correct the image This method does not amplify noise

and performs well in the case of inhomogeneous attenuating media

See Rosenthal et al and King et al for other methods for atten

uation correction

Resolution recovery

Image restoration techniques could be applied to nuclear medicine images to

perform resolution recovery deblurring and noise removal Most restoration

methods assume the presence of additive signalindependent white noise and

a shiftinvariant blurring function However nuclear medicine images are

corrupted by Poisson noise that is correlated with or dependent upon the

image data and the blurring function is shiftvariant Furthermore the noise

present in the planar projection images is amplied in SPECT images by the

reconstruction algorithm The SNR of nuclear medicine images is low which

makes recovery or restoration di cult

Several methods have been proposed for restoring nuclear medicine images

SPECT images may

be either restored after reconstruction postreconstruction restoration

or the projection planar images may be restored

rst and SPECT images reconstructed later prereconstruction restoration

see Figure

Boardman applied a constrained deconvolution lter to restore scinti

graphic images of the brain the lter required only the PSF as a priori infor

mation it was assumed that there were no discontinuities in the object Mad

sen and Park applied a Fourierdomain lter on the projection set for the

enhancement of SPECT images prereconstruction restoration they assumed

the PSF to be a D isotropic Gaussian function King et al

applied the Wiener lter and a countdependent Metz lter to restore both

projection and SPECT images and reported that prereconstruction restora

tion showed better results than postreconstruction restoration How

ever a study of image contrast and percent fractional standard deviation of

counts in regions of uniform activity in the test images used in their studies

does not show a clear dierence between the two methods Webb et al

attempted postreconstruction restoration of liver SPECT images they pro

posed a general formulation that unies a number of lters among which are

the inverse maximum entropy parametric Wiener homomorphic Phylips

and hybrid lters A study of image contrast in their work indicated that

properly tuned maximumentropy or homomorphic lters provide good re

sults Honda et al also attempted postreconstruction restoration of

myocardial SPECT images using a combination of Wiener and Butterworth

lters they assumed the SNR to be a constant and the system transfer func

tion to be Gaussian

Boulfelfel et al performed a comparison of prereconstruction

and postreconstruction restoration lters applied to myocardial images The

results obtained showed that prereconstruction restored images present a

signicant decrease in RMS error and an increase in contrast over post

reconstruction restored images Examples from these studies are presented

in Section

An important consideration in the restoration of SPECT images is the

stationarity of the PSF Hon et al and Boulfelfel

et al conducted several experiments to derive the FWHM values

of the PSF of several imaging congurations and tabulated the variation in

the parameter with sourcetocollimator distance and attenuating media

Larsson investigated the eects of using the arithmetic and geometric

mean of LSFs taken at opposing conjugate angular positions in air and water

to improve stationarity It was found that the arithmetic mean of opposing

LSFs resulted in an LSF of approximately the same FWHM as the LSF at

the center of rotation of the camera although there were signicant amplitude

dierences However the geometric mean of opposing LSFs provided compa

rable FWHMs and amplitudes Furthermore SPECT images reconstructed

by using the geometric means of opposing LSFs did not show any signicant

distortion due to the spatial variation in the detector response with distance

from the collimator as compared to the results with the arithmetic means

However Larsson found that when more than one source is employed the use

of the geometric mean resulted in interference artifacts between the sources

Axelsson et al investigated the shape of the LSF resulting from the geo

metric mean of opposing planar projections and concluded that the LSF was

approximately stationary in the central twothirds of their phantom Msaki

et al also found that the stationarity of the D MTF improved in their

study of the variation in the PSF of the geometric means of opposing views for

source positions both orthogonal and parallel to the collimator axis Boulfelfel

et al evaluated the use of the geometric mean of opposing projections

in prereconstruction restoration and found that this preprocessing step could

lead to improved results the details of their methods and results are presented

in Section

Coleman et al used the arithmetic and geometric means of opposing

projections to calculate the D MTF and scatter fraction as a function of the

camera angle and source location Their results showed that the geometric

mean provided an approximately stationary D MTF and scatter fraction

except near the edges of the phantom However it was observed by Coleman

et al that the arithmetic mean provided more nonstationary results than the

geometric mean which were both less nonstationary than the MTF related

to the planar projections It was concluded that the arithmetic mean may be

preferred to the geometric mean in the restoration of SPECT images because

the latter is nonlinear Although Coleman et al recommended the use of

averaging in restoration their study did not perform any restoration experi

ments King et al continuing the study of Coleman et al used the

means of conjugate views both arithmetic and geometric for prereconstruc

tion attenuation correction in their study on the use of the scatter degradation

factor and Metz ltering for the improvement of activity quantitation how

ever the Metz lter they used was predened and did not explicitly make use

of an MTF related to the averaging procedure Glick et al and

Coleman et al also investigated the eect of averaging arithmetic and

geometric of opposing views on the stationarity of the PSF Whereas these

studies concluded that the D PSF is not stationary and is shiftvariant the

D PSF has been shown to be more stationary and only slightly aected by

the sourcetocollimator distance

The PSF of SPECT images has been shown to have a D spread

regardless most of the restoration methods proposed in the literature

assume a D PSF The use of a D PSF results in only a partial restoration

of SPECT images because the interslice blur is ignored Boulfelfel et al

and Rangayyan et al studied the D nature of the PSF and

proposed D lters to address this problem examples from these works are

presented in Section

It has been shown that discretization of the ltered backprojection pro

cess can cause the MTF related to the blurring of SPECT images to be

anisotropic and nonstationary especially near the edges of the cameras eld

of view Furthermore the Poisson noise present in nuclear

medicine images is nonstationary Shiftinvariant restoration techniques will

fail in the restoration of large images because they do not account for such vari

ations in the MTF and the noise Therefore restoration methods for SPECT

images should support the inclusion of a shiftvariant MTF and nonstationary

noise characterization Boulfelfel et al investigated the use of

a shiftvariant D Kalman lter for the restoration of SPECT images The

Kalman lter allows for the use of dierent MTFs and noise parameters at

each pixel and was observed to perform better than shiftinvariant lters in

the restoration of large objects and organs Examples of this application are

presented in Section

Geometric averaging of conjugate projections

Geometric averaging of conjugate projections reviewed in Section in

brief and described in more detail in the following paragraphs may be viewed

as a predistortion mapping technique that consists of a transformation applied

to the degraded image in such a way that the shiftvariant blurring function

that caused the degradation becomes shiftinvariant The widely used coordi

nate transformation method is a predistortion mapping scheme

that eliminates the shiftvariance of a blurring function by changing to a sys

tem of coordinates in which the blur becomes shiftinvariant Shiftinvariant

restoration may then be used with the image changed back to the original

coordinates by an inverse coordinate transformation Coordinate transforma

tion methods have been applied for shiftvariant blurs such as motion blur

where the shiftvariance is linear However this technique is not applicable

in situations where the blur varies with the depth of eld or where a projec

tion integration operation is involved as in planar nuclear medicine images

Geometric averaging of conjugate projections could be interpreted as a pre

distortion mapping scheme that reduces the shiftvariance of the blur but not

eliminating it completely Furthermore because the procedure combines two

nearly identical planar projections images acquired from opposite directions

and averages them only the blur function is aected and the restored image

does not need any processing for coordinate change

In the work of Boulfelfel et al the MTF of the gamma camera was

modeled as

Hu v exp

N f

where N is the width of the MTF in pixels f

is the sampling frequency

and is the standard deviation of the Gaussian in the model of the PSF

The FWHM of the PSF was experimentally measured using a line source see

Section and Figure for various sourcetocollimator distances The

variation of the FWHM with sourcetocollimator distance d was found to be

linear and modeled as

FWHM a b d

where a and b are the parameters of the linear model that were determined

from experimental measurements of FWHM The standard deviation of the

Gaussian model of the PSF and MTF is related to FWHM as

FWHM

a b d

where is a constant of proportionality

Let us assume that the camera is rotating around a point source that is

located away from the center of rotation of the camera as shown in Fig

ure If the radius of rotation is R and the camera is at a distance d

from the point source when at the closest position then rotating the camera

places the point source at a distance Rd from the camera When

the point source is at a distance d from the camera the MTF according to

Equation is

u v exp

a b d

N f

where refers to the angle of the camera The MTF at the conjugate position

u v exp

a b R d

N f

where

refers to the angle of the camera The geometric mean

of the two MTFs given above is given by

u v H

u vH

u v

Therefore we have

u v exp

a b d

a b R d

exp

Rab b

Rd R

which leads to

u v exp

Rab b

Rd R

If the point source is located at the center of rotation of the camera we

have d R and the MTF is reduced to

u v exp

Rab b

Letting

Bu v

we have

u v exp

Bu v a

abd b

u v exp

Bu v

Rab b

Rd R

and

u v exp

Bu v a

Rab b

Equations and are respectively the MTF related to a

point source located at a distance d from a camera rotating around a circle of

radius R the geometric mean of the MTFs related to two opposing projections

of a point source located at distances d and Rd from the camera and the

MTF of a point source located at the center of rotation of the camera with

FIGURE

Gammacamera imaging geometry illustrating conjugate projections being ob

tained for a point source at distances d and R d from the camera in the

two views where R is the radius of rotation of the camera The same prin

ciples apply to imaging with a dualcamera or multicamera imaging system

Reproduced with permission from D Boulfelfel RM Rangayyan LJ Hahn

and R Kloiber Use of the geometric mean of opposing planar projections

in prereconstruction restoration of SPECT images Physics in Medicine and

Biology

IOP Publishing Ltd

d R Equation shows clearly that the MTF is highly dependent on

the sourcetocollimator distance whereas Equation suggests that the

geometric averaging procedure makes the MTF less dependent on the source

tocollimator distance only the last factor in the equation is a function of the

sourcetocollimator distance d and that it is close to the MTF of a point

source located at the center of rotation Equation

Comparing Equations and leads to the comparison

of the following three functions involving the sourcetocollimator distance d

and the other parameters encountered above

d a

abd b

d a

Rab b

Rd R

and

d a

Rab b

The subscripts of the function fd relate to the subscripts of Hu v in

Equations and Figure shows plots of the three

functions f

d f

d and f

d for radius of rotation R cm and

the FWHM model parameters a and b from exper

imental measurements using line sources with the Siemens Rota camera at

the Foothills Hospital Calgary It is seen that f

d a constant

d varies from to and f

d varies between and The

plots and the ranges of the values of the three functions show that geometric

averaging reduces the spacevariance of the MTF

FIGURE

Plots of the distance functions f

d in Equation solid line f

d in

Equation dashed line and f

d in Equation dotted line

Reproduced with permission from D Boulfelfel RM Rangayyan LJ Hahn

and R Kloiber Use of the geometric mean of opposing planar projections

in prereconstruction restoration of SPECT images Physics in Medicine and

Biology

IOP Publishing Ltd

Figure shows proles of the MTFs related to a point source at d

cm and d cm with R cm the averaged MTF with d cm

and d cm and the MTF at the center of rotation of the camera d R

computed using Equations and respectively The gure

shows that the averaged MTF and the MTF at the center of rotation are close

to each other as compared to the MTFs for the point source at d cm

and d cm

Boulfelfel et al performed experimental measurements with a line

source to verify the validity of the theoretical results described above The

line source was constructed with a thin plastic tube of internal radius of mm

and lled with mCi of

Tc No scattering medium was used Figure

FIGURE

Computed proles of MTFs related to point sources in gammacamera imag

ing a averaged MTF with d cm and d cm b MTF at the center

of rotation of the camera d R cm c d cm and d d cm

Reproduced with permission from D Boulfelfel RM Rangayyan LJ Hahn

and R Kloiber Use of the geometric mean of opposing planar projections

in prereconstruction restoration of SPECT images Physics in Medicine and

Biology

IOP Publishing Ltd

shows proles of the PSF derived from the LSF for sourcetocollimator dis

tances d and cm obtained using the ADAC GENESYS camera

with the lowenergy generalpurpose collimator at the Foothills Hospital Cal

gary The averaged PSF for d cm and cm is also plotted in the gure

It is seen that the averaged PSF matches closely the PSF at the central posi

tion of d cm

FIGURE

Experimentally measured proles of PSFs in gammacamera imaging a d

cm b at the center of rotation of the camera d R cm c

d cm and d averaged PSF with d cm and d cm Re

produced with permission from D Boulfelfel RM Rangayyan LJ Hahn

and R Kloiber Use of the geometric mean of opposing planar projections

in prereconstruction restoration of SPECT images Physics in Medicine and

Biology

IOP Publishing Ltd

The preceding derivations and arguments were based upon a point source

being imaged In order to extend the arguments to a combination of dis

tributed sources we could consider a planar source image px y that is par

allel to the plane of the camera A D source may then be modeled as a

collection of several planar sources with the individual planes being parallel

to the plane of the camera Then the acquired image of each plane could be

modeled as being convolved with the blur PSF for the corresponding distance

to the camera The net projection of the D source would be the sum of the

blurred images of the constituent planes

For a planar source px y placed away from the center of the axis of rota

tion of the camera let us consider two planar images acquired at an angle

and its conjugate

In the frequency domain we may represent the planar

images as

u v H

u vP u v

and

u v H

u vP u v

The geometric mean of the pair of conjugate planar images is given as

u v H

u vP u vH

u vP u v

u vH

u v

P u v

u vP u v

Therefore the geometric mean of a pair of conjugate planar images of a planar

source is equal to the original source distribution blurred by an MTF that

is given by the geometric mean of the individual MTFs This implies that

geometric means of pairs of conjugate planar images may be deconvolved by

using the geometric mean of the corresponding MTFs In practice it may be

appropriate to assume that the MTF is independent of the angular position

of the cameras then the same averaged MTF may be used for all angles

Pixelbypixel geometric averaging of opposing projections before prerecon

struction restoration can improve the performance of the restoration lter

because it reduces the spacevariance of the blur furthermore the averag

ing procedure reduces the eects of scatter and attenuation The averaging

technique is applicable when the object to be restored is of medium size over

which the averaged PSF or MTF may be assumed to be spaceinvariant

see Figure Prereconstruction restoration requires large computing

resources because each projection image needs to be restored Averaging re

duces the ltering time for restoration by because each opposing pair of

projections is replaced by a single image Geometric averaging of opposing

projections performs well in applications where the object is not located in a

corner of the eld of view of the camera It is required that projections be

acquired through the full range of

Geometric averaging reduces the shiftvariance of the blur function but

does not completely eliminate the variance Therefore artifacts due to the

shiftvariance of the blur function may remain in regions situated close to the

edges of the eld of view of the camera Prereconstruction restoration ltering

assumes that the blur function for all averaged projections is the same as for

points located at the axis of rotation of the camera Geometric averaging and

prereconstruction restoration procedures are wellsuited to SPECT imaging

of the brain where the image is centered and does not occupy the full eld of

view of the camera

Examples of the application of geometric averaging as a preprocessing step

prior to the restoration of SPECT images are presented in Section

Examples of restoration of SPECT images

Boulfelfel et al conducted several stud

ies on the restoration of SPECT images using the Wiener PSE Metz and

Kalman lters including the options of prereconstruction restoration post

reconstruction restoration and geometric averaging of the projections as well

as the application of the lters in D or D The MTFs of the imaging sys

tems used were experimentally measured using a line source to obtain the

LSF for various imaging parameters and congurations see Sections and

Some of their experiments and results are described in the following

paragraphs

Images of a tubular phantom A tubular phantom was constructed

with an acrylic tube of length cm and internal diameter mm within

which was introduced a solid acrylic rod of diameter mm The tube was

lled with a solution containing

Tc Sixty planar projections each of

size pixels were acquired spanning the angular range of

with the sourcetocollimator distances of cm and cm using a Siemens

Rota camera with a lowenergy allpurpose collimator No scattering medium

was used around the phantom other than the ambient atmosphere SPECT

images of size pixels each were reconstructed using the FBP algorithm

Models of the PSD of the phantom were mathematically derived for both

planar and SPECT imaging from the known geometry of the phantom Blob

like artifacts could arise in ltered images due to the discontinuities that are

present in discrete mathematical models as above In order to prevent

such artifacts the model PSDs were smoothed with a Gaussian window

of the form

W r exp

where r

is the index of radial frequency in the D Fourier space

N is the width of the PSD array in pixels and is a scale factor Under the

condition of MMSE of the restored images of known test images the optimal

value of was determined to be The PSD of the noise was derived from

the known total count in the image being restored

Due to the small size of the phantom in crosssection and due to its

placement at the center of the eld of imaging with respect to the axis of

rotation of the gamma camera it is reasonable to assume that the PSF of the

imaging system is stationary Hence shiftinvariant lters may be applied for

restoration

The Wiener and PSE lters were used for prereconstruction restoration of

the planar images in postreconstruction restoration the same lters were

also applied to the SPECT image Two sample projection images of the

phantom are shown in Figure for sourcetocollimator distances of cm

and cm it is evident that the latter image is blurred to a greater extent

than the former The restored version of the planar image with the source

tocollimator distance of cm using the Wiener lter is shown in part

c of Figure proles of the original image and the restored image are

shown in Figure The restored image clearly demonstrates the expected

reduction in counts at the center of the image due to the solid rod

Figure shows the original and restored SPECT images for various

cases of ltering Prereconstruction restoration gave better results than post

reconstruction restoration The Gaussian window applied to the object model

eectively reduced the hot spots blobs seen in the restored images without

the window Observe that the blobs appear with signicant amplitude in the

postreconstruction restored images but are reduced in the prereconstruction

restored images RMS error values were computed for the original and restored

images with respect to the ideal known crosssection of the phantom shown

in Figure b It was found that whereas the acquired image had an RMS

error of the error for the postreconstruction Wienerrestored image was

and that for the prereconstruction Wienerrestored image was

application of the Gaussian window to the model PSD further decreased the

error to The results of the PSE lter and Metz lter not shown had

comparable RMS errors

Images of a cardiac phantom A cardiac phantom Data Spectrum

Corporation was used to simulate myocardial perfusion mages with

defect inserts to simulate myocardial ischemia and infarction A sectional

view showing the dimensions of the phantom is illustrated in Figure The

phantom consists of two Ushaped barrels of diameter cm and cm The

space between the barrels was lled with MBq mCi of

T lchloride

Several types of defect arrangements were used Figure shows in cross

section a case with a

solid defect and a

defect with activity The

phantom was held inclined at

with respect to the axis of rotation of the

Rota camera with the lowenergy allpurpose collimator and projections of

size pixels each were acquired over

The acquisition time for each

projection was s and the radius of rotation was cm The SPECT images

of dierent transaxial slices were reconstructed and the D information was

used to derive oblique sectional images at

using the Siemens MicroDelta

software The total count for each image was between and

which corresponds to lowstatistics images in clinical studies

In performing restoration of the projection images of the cardiac phantom

an object PSD model needs to be derived If the model is derived directly from

the physical shape of the phantom it will match the actual phantom images

but will not be applicable for the restoration of real myocardial images with

unknown defects A model that matches closely the general case is a hollow

FIGURE

Projection planar images of the tubular phantom with the sourceto

collimator distance of a cm and b cm c Result of Wiener restoration

of the image in b Figures courtesy of D Boulfelfel See also Figures

and

FIGURE

Proles of planar images of the tubular phantom across the center Solid

line ideal true prole Dotted line prole of the acquired planar image

with the sourcetocollimator distance of cm see Figure b Dashed

line prole of the Wienerltered planar image see Figure c Figure

courtesy of D Boulfelfel

cylinder of appropriate length inclined in the same way as the phantom In

the works of Boulfelfel et al for each projection a rotation was performed on

such a model and the Radon integral was computed to model the projection

The acquired SPECT image of the phantom for the defect arrangement

illustrated in Figure as well as the restored images obtained by pre

reconstruction restoration using the Wiener and PSE lters are shown in

Figure The contrast values of the defect regions were computed as de

ned in Equation Groups of pixels representing the defect regions as well

as suitable background regions were selected manually over subimages

using the physical model of the phantom as depicted in Figure The

restoration lters resulted in more than doubling of the contrast values with

respect to the contrast of the same defect in the acquired image

Cardiac SPECT images One of the major applications of nuclear me

dicine imaging is in cardiology With the use of

T l myocardial diseases

and defects such as ischemia and necrosis are detected in SPECT slices as

cold spots of reduced tracer concentration Myocardial perfusion SPECT im

ages in the shortaxis view see Figure a appear similar to the rings

of activity in the images of the tubular phantom in Figure The car

diac phantom described in the preceding paragraphs may also be used to

simulate myocardial perfusion SPECT images related to various pathological

conditions see Figures and

a b c d

e f g h

FIGURE

a SPECT image of a tubular phantom size pixels RMS error

with respect to the known ideal version of the image shown in b Post

reconstruction restoration of the image in a using c the Wiener lter RMS

error and d the PSE lter RMS error SPECT image

after prereconstruction restoration of the planar images using e the Wiener

lter RMS error and f the PSE lter RMS error g

and h correspond to e and f with the inclusion of Gaussian smoothing of

the image PSD model RMS errors and respectively Images

courtesy of D Boulfelfel

FIGURE

Schematic representation of the cardiac phantom Data Spectrum Corpora

tion used to simulate myocardial perfusion images with defect inserts

to simulate myocardial ischemia and infarction The dimensions shown are in

a b c

FIGURE

a Acquired SPECT image of the cardiac phantom with the defect arrange

ment illustrated in Figure Result of prereconstruction restoration us

ing b the Wiener and c PSE lter Reproduced with permission from D

Boulfelfel RM Rangayyan LJ Hahn and R Kloiber Prereconstruction

restoration of single photon emission computed tomography images IEEE

Transactions on Medical Imaging

IEEE

Unfortunately myocardial SPECT images possess poor statistics because

only a small fraction of the injected activity will accumulate in the my

ocardium Furthermore as the peak photon energy of

T l is about keV

scattering has a serious eect on image quality Boulfelfel et al ap

plied prereconstruction restoration and postreconstruction restoration tech

niques using the Wiener PSE and Metz lters to myocardial SPECT images

examples from their works are presented and discussed in the following para

graphs

In the procedure to acquire myocardial SPECT images of human patients

MBq mCi of

T l was injected into the body After accumulation

of the tracer in the myocardium planar projections each of size

pixels spanning the full range of

were acquired The time for the

acquisition of each projection was s Each projection image had a total

count in the range to The projections were acquired in an el

liptic trajectory with the average distance from the heart being about cm

Two energy peaks were used in the acquisition of the projections in order

to perform scatter correction using the dualenergywindow subtraction tech

nique No attenuation correction was performed as the organ is small Given

the imaging protocol as above and the fact that the organ being imaged is

small it is valid to assume that the blurring function is nearly shiftinvariant

hence it becomes possible to apply shiftinvariant lters such as the Wiener

PSE and Metz lters for the restoration of myocardial planar and SPECT

images

Figures and show one projection image each of two patients

Transverse SPECT images were reconstructed and oblique slices perpendic

ular to the long axis of the heart were then computed from the D data

available Figures and show one representative oblique section

image in each case along with several restored versions of the images The

parts of the myocardium with reduced activity cold spots are seen more

clearly in the restored images than in the original images The results of pre

reconstruction restoration applied to the planar images are better than those

of postreconstruction restoration ltering in terms of noise content as well as

improvement in sharpness and clarity

SPECT images of the brain Radionuclide brain scanning has been used

extensively in the study of neurological and psychiatric diseases The main

area of application is the detection of pathology in the cerebral hemispheres

and the cerebellum A number of radiopharmaceuticals are used for brain

scanning however

Tcbased materials are most widely used

An advantage in brain imaging is that the patients head can be positioned

at the center of rotation of the camera which allows imaging over

rather

than over only

as in the case of myocardial imaging Although the brain

may be considered to be a large organ the homogeneity of the scattering

medium and the ability to image it from a short distance using a circular orbit

allow the use of geometric averaging of the planar images as a preprocessing

FIGURE

Top A sample planar projection image of a patient a Shortaxis SPECT

image showing the myocardium of the left ventricle in crosssection Re

sults of postreconstruction restoration applied to the SPECT image using

b the Wiener c the PSE and d Metz lters Results of prereconstruc

tion restoration applied to the planar images using e the Wiener f the

PSE and g Metz lters Images courtesy of D Boulfelfel LJ Hahn and

R Kloiber Foothills Hospital Calgary

FIGURE

Top A sample planar projection image of a patient a Shortaxis SPECT

image showing the myocardium of the left ventricle in crosssection Re

sults of postreconstruction restoration applied to the SPECT image using

b the Wiener c the PSE and d Metz lters Results of prereconstruc

tion restoration applied to the planar images using e the Wiener f the

PSE and g Metz lters Images courtesy of D Boulfelfel LJ Hahn and

R Kloiber Foothills Hospital Calgary

step to reduce both attenuation and shiftvariance of the blur before restora

tion Brain images are also not as low in statistics as myocardial images

In the procedure for nuclear medicine imaging of the brain after the

chloride administered to the patient had accumulated in the brain planar

projections each of size pixels were acquired The time for the

acquisition of each projection was s The projections were acquired over the

full range of

in a circular trajectory with the radius of rotation of cm

Two energy peaks were used in the acquisition of the projections to perform

scatter correction using the dualenergywindow subtraction technique

Figures and show a set of opposing projection images as well

as their geometric mean for two patients Transverse SPECT images were re

constructed after performing geometric averaging of conjugate projections and

prereconstruction restoration using the Wiener PSE and Metz lters

Figures and show one representative SPECT image in each case

along with several restored versions The results show that averaging of conju

gate projections improves the quality of the restored images which are sharper

than the images restored without averaging

Images of a resolution phantom Boulfelfel et al

conducted several restoration experiments with SPECT images of a resolu

tion phantom The phantom contains nine pairs of hot spots of diameters

and mm in the hot lesion insert Nuclear Asso

ciates with a total diameter of mm see Figure for related illustra

tions The phantom was lled with mCi of

T lchloride centered at the

axis of rotation of the gamma camera at a distance of mm and pro

jections each of size pixels were acquired over

SPECT images

of dierent transaxial slices were reconstructed using the Siemens MicroDelta

software

Given the large size of the phantom it would be inappropriate to assume

that the degradation phenomena are shiftinvariant Boulfelfel et al

applied the Kalman lter for restoration of SPECT images of the resolution

phantom Figure shows a representative SPECT image of the phantom

along with postreconstruction restoration of the image using the Kalman

lter the results of application of the shiftinvariant Wiener and PSE lters

are also shown for comparison It is evident that the shiftvariant Kalman

lter has provided better results than the other lters the Kalmanrestored

image clearly shows seven of the nine pairs of hot spots whereas the results

of the Wiener and PSE lters show only four or ve pairs For the sake

of comparison the results of prereconstruction restoration of the resolution

phantom image obtained by applying the shiftinvariant Wiener PSE and

Metz lters after geometric averaging of conjugate projections are shown in

Figure Observe that the orientation of these results is dierent from

that of the images in Figure due to the alignment procedure required

for averaging Although the results show some of the hot spots with more

clarity than the original image in Figure a they are of lower quality

than the result of Kalman ltering shown in Figure d

FIGURE

Top row A sample pair of conjugate projections of a patient along with their

geometric mean a SPECT image showing the brain in crosssection Results

of prereconstruction restoration applied to the planar images using b the

Wiener c the PSE and d Metz lters Results of geometric averaging

and prereconstruction restoration applied to the planar images using e the

Wiener f the PSE and g Metz lters The orientation of the images in

e g is dierent from that of the images in a d due to the alignment

of conjugate projection images for geometric averaging Images courtesy of

D Boulfelfel LJ Hahn and R Kloiber Foothills Hospital Calgary