GACS-7205-001 Digital Image Processing Page 285

(Winter Term, 2021-22)

5.7 Inverse Filtering

The simplest approach to restoration is direct inverse filtering,

where we compute an estimate, Fˆ(u, v) , of the transform of the
original image by
G(u, v)
Fˆ(u, v) = (5.7-1)
H (u, v)

Substituting the right side of

G(u, v ) = H (u, v )F (u, v ) + N (u, v ) (5.1-2)

in (5.7-1) yields
N (u, v)
Fˆ(u, v) = F (u, v ) +
H (u, v) . (5.7-2)

The bad news is that we cannot recover the undegraded image

exactly because N (u, v ) is not known.

More bad news is that if the degradation function H (u, v) has zero
or very small values, so the second term of (5.7-2) could easily
dominate the estimate of Fˆ(u, v) .

One approach to get around the zero or small-value problem is to

limit the filter frequencies to values near the origin. As discussed
earlier, we know that H (0, 0) is usually the highest value of
H (u, v ) in the frequency domain.
GACS-7205-001 Digital Image Processing Page 286
(Winter Term, 2021-22)

Example 5.11: Inverse filtering

The image in Figure 5.25 (b) was inverse filtered with

G(u, v)
Fˆ(u, v) = (5.7-1)
H (u, v)
using the exact inverse of the degradation function that generated
that image. That is, the degradation function used was
5/6
-k éëê ( u -M /2 )2 +( v -N /2 )2 ùûú
H (u, v ) = e
with k = 0.0025 . In this case, M = N = 480 .

Although a Gaussian-shape function has no zeros and it is not a

concern here, the degradation values become so small that the
result of full inverse filtering shown in Figure 5.27 (a) is useless.
GACS-7205-001 Digital Image Processing Page 287
(Winter Term, 2021-22)

Figure 5.27 (b) through (d) show the results of cutting off
values of the ratio G ( u, v ) / H ( u, v ) outside a radius of 40, 70,
and 85, respectively.

Values above 70 started to produce degraded images, and

further increases in radius values would produce images that
looked more and more like Figure 5.27 (a).

The results in Example 5.11 show the poor performance of direct

inverse filtering in general.
GACS-7205-001 Digital Image Processing Page 288
(Winter Term, 2021-22)

5.8 Minimum Mean Square Error (Wiener) Filtering

Here we discuss an approach that incorporates both the degradation

function and statistical characteristics of noise into the restoration
process.

Considering images and noise as random variables, the objective

is to find an estimate fˆ of the uncorrupted image f such that the
mean square error between them is minimized.

The error measure is given by

e 2 = E { (f - fˆ)2 } (5.8-1)

where E {} is the expected value of the argument.

By assuming that
1. the noise and the image are uncorrelated;
2. one or the other has zero mean;
3. the intensity levels in the estimate are a linear function of
the levels in the degraded image.

Then, the minimum of the error function in (5.8-1) is given in the

frequency domain by the expression
é H * (u, v )S f (u, v) ù
ˆ
F (u, v) = ê ê ú G(u, v )
2 ú
êë S f (u, v ) H (u, v) + S h (u, v) úû
é H * (u, v) ù
=êê ú G(u, v )
2 ú (5.8-2)
êë H (u, v) + S h (u, v)/ S f (u, v) úû
é 1 H (u, v ) 2 ù
=ê ê ú G(u, v )
ëê
H (u , v ) H (u, v) + S h (u, v)/ S f (u, v) úûú
2
GACS-7205-001 Digital Image Processing Page 289
(Winter Term, 2021-22)

The terms in (5.8-2) are as follows:

Fˆ(u, v ) is the frequency domain estimate

G(u, v ) is the transform of the degraded image

H (u, v ) is the transform of the degradation function

H * (u, v ) is complex conjugate of H (u, v )

H (u, v ) 2 = H * (u, v )H (u, v )

S h (u, v ) = N (u, v ) 2 = power spectrum of the noise

S f (u, v ) = F (u, v ) 2 = power spectrum of the undegraded image

This result is known as the Wiener filter, which also is commonly

referred to as the minimum mean square error filter or the least
square error filter.

The Wiener filter does not have the same problem as the inverse
filter with zeros in the degradation function, unless the entire
denominator is zero for the same value(s) of u and v .

If the noise is zero, then the Wiener filter reduces to the inverse
filter.

One of the most important measures is the signal-to-noise ratio,

approximated using frequency domain quantities such as
M -1 N -1

åå F (u, v) 2
u =0 v =0
SNR = M -1 N -1
(5.8-3)
åå N (u, v) 2
u =0 v =0
GACS-7205-001 Digital Image Processing Page 290
(Winter Term, 2021-22)

The mean square error given in statistical form in (5.8-1) can be

approximated also in terms a summation involving the original
and restored images:
M -1 N -1
1 2
MSE =
MN å å éêë f (x, y) - fˆ(x, y) ùúû (5.8-4)
x =0 y =0

If one considers the restored image to be signal and the difference

between this image and the original to be noise, we can define a
signal-to-noise ratio in the spatial domain as
M -1 N -1

å å fˆ(x, y)2
x =0 y =0
SNR = M -1 N -1
2 (5.8-5)
å å ëê
é f (x , y ) - ˆ
f (x , y ) ù
ûú
x =0 y =0

The closer f and fˆ are, the larger this ratio will be.
2
If we are dealing with white noise, the spectrum N (u, v ) is a
constant, which simplifies things considerably. However,
F (u, v ) 2 is usually unknown.

An approach is used frequently when these quantities are not

known or cannot be estimated:
é 1 H (u, v) 2 ù
ˆ
F (u, v ) = ê ú G(u, v )
ê H (u, v ) H (u, v ) 2 + K ú (5.8-6)
ë û

where K is a specified constant that is added to all terms of

H (u, v ) 2 .

Note: White noise is a random signal (or process) with a flat power spectral
density. In other words, the signal contains equal power within a fixed
bandwidth at any center frequency.
GACS-7205-001 Digital Image Processing Page 291
(Winter Term, 2021-22)

Example 5.12: Comparison of inverse and Wiener filtering

Figure 5.28 shows the advantage of Wiener filtering over direct

inverse filtering.

Figure 5.28 (a) is the full inverse-filtered result from Figure 5.27 (a).

Figure 5.28 (b) is the radially limited inverse result of Figure 5.27 (c).

Figure 5.28 (c) shows the result obtained using

é 1 H (u, v) 2 ù
ˆ
F (u, v ) = ê ú G(u, v )
ê H (u, v ) H (u, v ) 2 + K ú (5.8-6)
ë û
with the degradation function
5/6
-k éêë ( u -M /2 )2 +( v -N /2 )2 ùúû
H (u, v ) = e

used in Example 5.11. The value of K was chosen interactively

to yield the best visual result.

By comparing Figure 5.25 (a) and Figure 5.28 (c), we see that
the Wiener filter yielded a result very close in appearance to the
original image.
GACS-7205-001 Digital Image Processing Page 292
(Winter Term, 2021-22)

Example 5.13: Further comparisons of Wiener filtering

GACS-7205-001 Digital Image Processing Page 293
(Winter Term, 2021-22)

5.10 Geometric Mean Filter

It is possible to generalize the Wiener filter slightly to the so-called

geometric mean filter:
é ù1-a
ê ú
é H * (u, v) ù a ê H * (u, v) ú
Fˆ(u, v) = ê ú ê ú G(u, v )
ê H (u, v ) 2 ú ê é S h (u, v ) ù ú , (5.10-1)
ë û 2
ê H (u, v ) + b ê úú
êë êë S f (u, v) úû úû

where a and b are positive real constants.

If a = 1 , this filter reduces to the inverse filter.

If a = 0 , the filter becomes the so-called parametric Wiener

filter, which reduces to the standard Wiener filter when b = 1 .

If a = 1/ 2 , this filter becomes a product of the two quantities

raised to the same power, which is the definition of the geometric
mean. When b = 1 , the filter is also commonly referred to as the
spectrum equalization filter.

With b = 1 , as a decreases below 1/ 2 , the filter performance

will tend more toward to the inverse filter; as a increases above
1/ 2 , the filter will behave more like the Wiener filter.
GACS-7205-001 Digital Image Processing Page 294
(Winter Term, 2021-22)

5.11 Image Reconstruction from Projections

In this section, we will examine the problem of reconstructing

an image from a series of projections, with a focus on X-ray
computed tomography (CT), which is one of the principal
applications of digital image processing in medicine.

Introduction

Consider Figure 5.32 (a), which consists of a single object on a

uniform background.

Suppose that we pass a thin, flat beam of X-rays from left to right,
and assume that the energy of the beam is absorbed more by the
object than by the background. Using a strip of X-ray absorption
detectors on the other side will yield the signal, whose amplitude
(intensity) is proportional to absorption.

The approach is to project the 1-D signal back across the direction
from which the beam came, as Figure 5.32 (b) shows. This approach
is called backprojection.
GACS-7205-001 Digital Image Processing Page 295
(Winter Term, 2021-22)

We certainly cannot determine a single object or a multitude of

objects along the path of the beam by a single project.

If we rotate the position of the source-detector pair by 90o and

repeat the previous procedure, we will get a backprojection image
shown in Figure 5.32 (d). Adding this result to Figure 5.32 (b)
will result an image illustrated in Figure 5.32 (e).

We should be able to learn more about the shape of the object in

question by taking more views in the same manner, as shown in
Figure 5.33.

As the number of projections increases, the strength of non-

intersecting backprojects deceases relative to the strength of
regions in which multiple backprojects intersect.

Figure 5.33 (f) shows the result formed from 32 projections.

The reconstructed image seems to be a reasonably good

approximation to the shape of the original object. However, the
image is blurred by a “halo” effect, which shows a “star” in Figure
5.33 (e). As the number of views increases, the shape of the “halo”
becomes circular, as shown in Figure 5.33 (e).
GACS-7205-001 Digital Image Processing Page 296
(Winter Term, 2021-22)

Blurring in CT reconstruction is an important issue and will be

addressed in later discussion.

Since the projections 180o apart are mirror images of each other,
we only need to consider angle increments halfway around a circle
in order to generate all the projects required for reconstruction.

Example 5.16: Backprojection of a simple planar region containing

two objects
GACS-7205-001 Digital Image Processing Page 297
(Winter Term, 2021-22)

Principles of Computed Tomography (CT)

The theoretical foundation of CT dates back to Johann Radon, a

mathematician from Vienna who derived a method in 1907 for
projecting a 2-D object along parallel rays as part of his work on
line integrals. The method is referred as the Radon transform now.

Allan M. Cormack, a physicist at Tufts University, partially

“rediscovered” these concepts and applied them to CT. Cormack
published his initial findings in 1963 and 1964. He provided the
mathematical formulae needed for the reconstruction and built a
CT prototype to show his ideas.

Working independently, electrical engineer Godfrey N. Hounsfield

and his colleagues at EMI in London formulated a similar solution
and built the first medical CT machine.

Cormack and Hounsfield shared the 1979 Nobel Prize in Medicine

for their contributions to medical tomography.

Figure 5.35 shows the first four generations of CT scanners.

GACS-7205-001 Digital Image Processing Page 298
(Winter Term, 2021-22)

The fifth-generation (G5) CT scanners eliminate all mechanical

motion by employing electron beams controlled electromagnetically.

The sixth-generation (G6) CT scanners rotate the source-detector

pair continuously through 360o , while the patient is moved at a
constant speed along the axis perpendicular to the scan.

The seventh-generation (G7) CT scanners (also called multislice CT

scanners) use parallel banks of detectors to collect volumetric CT
data simultaneously.

Projections and the Radon Transform

A straight line in Cartesian coordinates can be described either by

its slope-intercept form
y = ax + b ,
or, as in Figure 5.36, by its normal representation
x cos q + y sin q = r . (5.11-1)

Figure 5.36 Normal representation of a straight line.

GACS-7205-001 Digital Image Processing Page 299
(Winter Term, 2021-22)

The projection of a parallel-ray beam may be modeled by a set of

such lines, as shown in Figure 5.37.

Figure 5.37 Geometry of a parallel-ray beam.

An arbitrary point in the projection signal is given by the raysum

along the line
x cos qk + y sin qk = r j .

In the case of continuous, the raysum is a line integral, given by

¥ ¥
g(r j , qk ) = ò-¥ ò-¥ f (x, y)d(x cos qk + y sin qk - r j )dxdy (5.11-2)

Recall the properties of the impulse, d , the right side of (5.11-2)

is zero unless the argument of d is zero. It indicates that the
integral is computed only along the line x cos qk + y sin qk = r j .
GACS-7205-001 Digital Image Processing Page 300
(Winter Term, 2021-22)

If we consider all values of r and q , (5.11-2) generalizes

¥ ¥
g(r, q) = ò-¥ ò-¥ f (x, y)d(x cos q + y sin q - r)dxdy . (5.11-3)

The equation (5.11-3) gives the projection of f (x, y ) along an

arbitrary line in the xy -plane , is called the Radon transform.

The Radon transform is the cornerstone of reconstruction from

projections, with CT being its principle application in the field of
image processing.

In the discrete case, (5.11-3) becomes

M -1 N -1
g(r, q) = å å f (x, y)d ( x cos q + y sin q - r ) , (5.11-4)
x =0 y =0

where x , y , r , and q are now discrete variables.

If we fix q and allow r to vary, (5.11-4) simply sums the pixels

of f (x, y ) along the line defined by specified values of these two
parameters.

Incrementing through all values of r required to span the image

(with q fixed) yields one projection. Changing q and repeating
the same procedure will yield another projection.
GACS-7205-001 Digital Image Processing Page 301
(Winter Term, 2021-22)

Example 5.17: Using the Radon transform to obtain the projection

of a circular region.

We want to obtain the Radon transform for the projection of

the circular object
ì
ï A x 2 + y2 £ r 2
ï
f (x , y ) = í
ï 0 otherwise ,
ï
î
where A is a constant and r is the radius of the object. The
circular object is shown in Figure 5.38 (a).

Since the object is circularly symmetric, its projections are the

same for all angles, so all we need is to obtain the projection
for q = 0o . From (5.11-3), we get
¥ ¥
g(r, q) = ò-¥ ò-¥ f (x, y)d(x - r)dxdy
¥
= ò-¥ f (r, y)dy
This is a line integral along the line L ( r, 0 ) .
GACS-7205-001 Digital Image Processing Page 302
(Winter Term, 2021-22)

Note that g ( r, q ) = 0 when r > r . When r £ r , the

2 2 2 2
integral is evaluated from y = - r - r to y = r - r .

Therefore,

r 2 -r2 r 2 -r2
g(r, q) = ò- r 2 -r2
f (r, y )dy = ò- r 2 -r2
Ady

It yields
ìï 2A r 2 - r2 r £r
ï
g ( r, q ) = g ( r ) = í
ïï 0 otherwise
ïî
Figure 5.38 (b) shows the result.
g ( r, q ) = g ( r ) indicates that g is independent of q because
the object is symmetric about the origin.
GACS-7205-001 Digital Image Processing Page 303
(Winter Term, 2021-22)

When the Radon transform, g ( r, q ) , is displayed as an image

with r and q as rectilinear coordinates, the result is called a
sinogram, similar in concept to displaying the Fourier spectrum.
Like the Fourier spectrum, a sinogram contains the data necessary
to reconstruct f ( x, y ) .

Figure 5.39 (b) is the sinogram of the rectangle shown in Figure

5.39 (a).

Figure 5.39 (c) shows an image of the Shepp-Logan phantom, a

widely used synthetic image designed to simulate the absorption
of major areas of the brain. The sinogram of Figure 5.39 (c) is
shown in Figure 5.39 (d).
GACS-7205-001 Digital Image Processing Page 304
(Winter Term, 2021-22)

To obtain a formal expression for a back-projected image from

Radon transform, referring to Figure 5.37, we begin with a single
point, g ( r j , qk ) , of the complete projection, g ( r, qk ) , for a fixed
value of rotation, qk .

Forming part of an image by back-projecting this single point is

simply to copy the line L ( r j , qk ) onto the image, where the value
of each point in that line is g ( r j , qk ) . Repeating this process of
all values of r j in the projected signal results
fqk ( x , y ) = g ( r, qk ) = g ( x cos qk + y sin qk , qk ) .
This equation holds for an arbitrary value of qk , therefore, we
can write in general that the image formed from a single
backprojection obtained at an angle q is given by
fq (x, y ) = g ( x cos q + y sin q, q ) . (5.11-5)
We form the final image by integrating over all the back-projected
images
p
f (x, y ) = ò0 fq ( x , y )d q (5.11-6)

In the discrete case, the integral becomes a sum of all back-

projected images:
p
f (x , y ) = å fq ( x, y ) (5.11-7)
q= 0

For example, if 0.5o increments are being used, the summation is

from 0 to 179.5 .

A back-projected image formed in this manner is referred to as a

laminogram, which is only an approximation to the image from
which the projections were generated.
GACS-7205-001 Digital Image Processing Page 305
(Winter Term, 2021-22)

Example 5.18: Obtaining back-projected images from sinograms

Equation
p
f (x , y ) = å fq ( x, y ) (5.11-7)
q= 0

was used to generate the back-projected images in Figure 5.32

through Figure 5.34 from projections obtained with
M -1 N -1
g(r, q) = å å f (x, y)d ( x cos q + y sin q - r ) . (5.11-4)
x =0 y =0

These equations were also used to generate Figure 5.40 (a) and
Figure 5.40 (b), which show the back-projected images
corresponding to the sinograms in Figure 5.39 (b) and Figure
5.39 (d).

Note that there is a significant amount of blurring shown in

Figure 5.40 (a) and (b). It is obvious that a straight use of
Equations (5.11-4) and (5.11-7) will not yield acceptable results.