318 Chapter 5 Image Restoration and Reconstruction

5.1 A MODEL OF THE IMAGE DEGRADATION/RESTORATION

5.1
PROCESS
In this chapter, we model image degradation as an operator Ᏼ that, together with an
additive noise term, operates on an input image f ( x, y) to produce a degraded image
g( x, y) (see Fig. 5.1). Given g( x, y), some knowledge about Ᏼ, and some knowledge
about the additive noise term h( x, y), the objective of restoration is to obtain an
estimate fˆ ( x, y) of the original image. We want the estimate to be as close as possible
to the original image and, in general, the more we know about Ᏼ and h, the closer
fˆ ( x, y) will be to f ( x, y).
We will show in Section 5.5 that, if Ᏼ is a linear, position-invariant operator, then
the degraded image is given in the spatial domain by

g( x, y) = (h 夹 f )( x, y) + h( x, y) (5-1)

where h( x, y) is the spatial representation of the degradation function. As in Chapters

3 and 4, the symbol “夹 ” indicates convolution. It follows from the convolution theorem
that the equivalent of Eq. (5-1) in the frequency domain is

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

where the terms in capital letters are the Fourier transforms of the corresponding
terms in Eq. (5-1). These two equations are the foundation for most of the restora-
tion material in this chapter.
In the following three sections, we work only with degradations caused by noise.
Beginning in Section 5.5 we look at several methods for image restoration in the
presence of both Ᏼ and h.

5.2 NOISE MODELS

5.2

The principal sources of noise in digital images arise during image acquisition and/or
transmission. The performance of imaging sensors is affected by a variety of environ-
mental factors during image acquisition, and by the quality of the sensing elements
themselves. For instance, in acquiring images with a CCD camera, light levels and
sensor temperature are major factors affecting the amount of noise in the resulting
image. Images are corrupted during transmission principally by interference in the
transmission channel. For example, an image transmitted using a wireless network
might be corrupted by lightning or other atmospheric disturbance.

FIGURE 5.1
A model of the g(x, y)
f (x, y) Degradation Restoration fˆ(x, y)
image
Ᏼ filter(s)
degradation/
restoration
Noise
process.
h(x, y)

DEGRADATION RESTORATION

[Link]

DIP4E_GLOBAL_Print_Ready.indb 318 6/16/2017 2:06:57 PM

 5.2 Noise Models 319

SPATIAL AND FREQUENCY PROPERTIES OF NOISE

Relevant to our discussion are parameters that define the spatial characteristics of
noise, and whether the noise is correlated with the image. Frequency properties refer
to the frequency content of noise in the Fourier (frequency) domain discussed in
detail in Chapter 4. For example, when the Fourier spectrum of noise is constant, the
noise is called white noise. This terminology is a carryover from the physical prop-
erties of white light, which contains all frequencies in the visible spectrum in equal
proportions.
With the exception of spatially periodic noise, we assume in this chapter that
noise is independent of spatial coordinates, and that it is uncorrelated with respect
to the image itself (that is, there is no correlation between pixel values and the values
of noise components). Although these assumptions are at least partially invalid in
some applications (quantum-limited imaging, such as in X-ray and nuclear-medicine
imaging, is a good example), the complexities of dealing with spatially dependent
and correlated noise are beyond the scope of our discussion.

SOME IMPORTANT NOISE PROBABILITY DENSITY FUNCTIONS

You may find it helpful

In the discussion that follows, we shall be concerned with the statistical behavior of
to take a look at the the intensity values in the noise component of the model in Fig. 5.1. These may be
Tutorials section of the
book website for a brief
considered random variables, characterized by a probability density function (PDF),
review of probability. as noted briefly as noted earlier. The noise component of the model in Fig. 5.1 is an
image, h( x, y), of the same size as the input image. We create a noise image for simu-
lation purposes by generating an array whose intensity values are random numbers
with a specified probability density function. This approach is true for all the PDFs
to be discussed shortly, with the exception of salt-and-pepper noise, which is applied
differently. The following are among the most common noise PDFs found in image
processing applications.

Gaussian Noise
Because of its mathematical tractability in both the spatial and frequency domains,
Gaussian noise models are used frequently in practice. In fact, this tractability is so
convenient that it often results in Gaussian models being used in situations in which
they are marginally applicable at best.
The PDF of a Gaussian random variable, z, is defined by the following familiar
expression:
( z − z )2
1 −
p(z) = e 2s2 − < z <  (5-3)
2ps

where z represents intensity, z is the mean (average) value of z, and s is its standard
deviation. Figure 5.2(a) shows a plot of this function. We know that for a Gaussian
random variable, the probability that values of z are in the range z ± s is approxi-
mately 0.68; the probability is about 0.95 that the values of z are in the range z ± 2s.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 319 6/16/2017 2:06:58 PM

 320 Chapter 5 Image Restoration and Reconstruction

p(z) p(z) p(z)

1 2
0.607 K
2ps b
Gaussian Rayleigh Erlang (Gamma)
0.607
2ps a(b  1)b1 (b1)
K e
(b  1)!

_ _ _ z z
zs z z s a z (b  1)/a
b
a
2
p(z) p(z) p(z)

1 1 − (Ps + Pp )
a ba
Exponential Salt-and-
Uniform
Ps pepper
Pp

z a b z 0 V 2k − 1 z
a b c
d e f
FIGURE 5.2 Some important probability density functions.

Rayleigh Noise
The PDF of Rayleigh noise is given by

⎧2 − ( z − a )2
⎪ ( z − a) e z≥ a
b
p(z) = ⎨ b (5-4)
⎪⎩0 z< a

The mean and variance of z when this random variable is characterized by a Ray-
leigh PDF are

z=a + pb 4 (5-5)
and
b (4 − p)
s2 = (5-6)
4
Figure 5.2(b) shows a plot of the Rayleigh density. Note the displacement from the
origin, and the fact that the basic shape of the density is skewed to the right. The
Rayleigh density can be quite useful for modeling the shape of skewed histograms.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 320 6/16/2017 2:06:58 PM

 5.2 Noise Models 321

Erlang (Gamma) Noise

The PDF of Erlang noise is

⎧ ab zb − 1 − az
⎪ e z≥ 0
p(z) = ⎨ (b − 1)! (5-7)
⎪0 z< 0
⎩
where the parameters are such that a > b, b is a positive integer, and “!” indicates
factorial. The mean and variance of z are
b
z= (5-8)
a
and
b
s2 = (5-9)
a2
Figure 5.2(c) shows a plot of this density. Although Eq. (5-9) often is referred to as
the gamma density, strictly speaking this is correct only when the denominator is
the gamma function, (b). When the denominator is as shown, the density is more
appropriately called the Erlang density.

Exponential Noise
The PDF of exponential noise is given by

− az
⎪⎧ae z≥ 0
p(z) = ⎨ (5-10)
⎩⎪0 z< 0

where a > 0. The mean and variance of z are

1
z= (5-11)
a
and
1
s2 = (5-12)
a2
Note that this PDF is a special case of the Erlang PDF with b = 1. Figure 5.2(d)
shows a plot of the exponential density function.

Uniform Noise
The PDF of uniform noise is

⎧ 1
⎪ a≤ z≤b
p(z) = ⎨ b − a (5-13)
⎪0 otherwise
⎩

[Link]

DIP4E_GLOBAL_Print_Ready.indb 321 6/16/2017 2:07:00 PM

 322 Chapter 5 Image Restoration and Reconstruction

The mean and variance of z are

a+ b
z= (5-14)
2
and
(b − a)2
s2 = (5-15)
12
Figure 5.2(e) shows a plot of the uniform density.

Salt-and-Pepper Noise
If k represents the number of bits used to represent the intensity values in a digital
image, then the range of possible intensity values for that image is [0, 2 k − 1] (e.g.,
[0, 255] for an 8-bit image). The PDF of salt-and-pepper noise is given by
When image intensities
are scaled to the range
[0, 1], we replace by 1 the ⎧ Ps for z = 2 k − 1
value of salt in this equa- ⎪
tion. V then becomes a p(z) = ⎨ Pp for z = 0 (5-16)
fractional value in the ⎪
open interval (0, 1).
⎩1 − (Ps + Pp ) for z = V

where V is any integer value in the range 0 < V < 2 k − 1.

Let h( x, y) denote a salt-and-pepper noise image, whose intensity values satisfy
Eq. (5-16). Given an image, f ( x, y), of the same size as h( x, y), we corrupt it with salt-
and-pepper noise by assigning a 0 to all locations in f where a 0 occurs in h. Similarly,
we assign a value of 2 k − 1 to all location in f where that value appears in h. Finally,
we leave unchanged all location in f where V occurs in h.
If neither Ps nor Pp is zero, and especially if they are equal, noise values satisfy-
ing Eq. (5-16) will be white (2 k − 1) or black (0), and will resemble salt and pepper
granules distributed randomly over the image; hence the name of this type of noise.
Other names you will find used in the literature are bipolar impulse noise (unipolar
if either Ps or Pp is 0), data-drop-out noise, and spike noise. We use the terms impulse
and salt-and-pepper noise interchangeably.
The probability, P, that a pixel is corrupted by salt or pepper noise is P = Ps + Pp .
It is common terminology to refer to P as the noise density. If, for example, Ps = 0.02
and Pp = 0.01, then P = 0.03 and we say that approximately 2% of the pixels in an
image are corrupted by salt noise, 1% are corrupted by pepper noise, and the noise
density is 3%, meaning that approximately 3% of the pixels in the image are cor-
rupted by salt-and-pepper noise.
Although, as you have seen, salt-and-pepper noise is specified by the probability
of each, and not by the mean and variance, we include the latter here for complete-
ness. The mean of salt-and-pepper noise is given by

z = (0)Pp + K (1 − Ps − Pp ) + (2 k − 1)Ps (5-17)

and the variance by

[Link]

DIP4E_GLOBAL_Print_Ready.indb 322 6/16/2017 2:07:01 PM

 5.2 Noise Models 323

s 2 = (0 − z )2 Pp + (K − z )2 (1 − Ps − Pp ) + (2 k − 1)2 Ps (5-18)

where we have included 0 as a value explicit in both equations to indicate that the
value of pepper noise is assumed to be zero.
As a group, the preceding PDFs provide useful tools for modeling a broad range
of noise corruption situations found in practice. For example, Gaussian noise arises
in an image due to factors such as electronic circuit noise and sensor noise caused by
poor illumination and/or high temperature. The Rayleigh density is helpful in char-
acterizing noise phenomena in range imaging. The exponential and gamma densities
find application in laser imaging. Impulse noise is found in situations where quick
transients, such as faulty switching, take place during imaging. The uniform density
is perhaps the least descriptive of practical situations. However, the uniform density
is quite useful as the basis for numerous random number generators that are used
extensively in simulations (Gonzalez, Woods, and Eddins [2009]).

EXAMPLE 5.1 : Noisy images and their histograms.

Figure 5.3 shows a test pattern used for illustrating the noise models just discussed. This is a suitable pat-
tern to use because it is composed of simple, constant areas that span the gray scale from black to near
white in only three increments. This facilitates visual analysis of the characteristics of the various noise
components added to an image.
Figure 5.4 shows the test pattern after addition of the six types of noise in Fig. 5.2. Below each image
is the histogram computed directly from that image. The parameters of the noise were chosen in each
case so that the histogram corresponding to the three intensity levels in the test pattern would start to
merge. This made the noise quite visible, without obscuring the basic structure of the underlying image.
We see a close correspondence in comparing the histograms in Fig. 5.4 with the PDFs in Fig. 5.2.
The histogram for the salt-and-pepper example does not contain a specific peak for V because, as you
will recall, V is used only during the creation of the noise image to leave values in the original image
unchanged. Of course, in addition to the salt and pepper peaks, there are peaks for the other intensi-
ties in the image. With the exception of slightly different overall intensity, it is difficult to differentiate

FIGURE 5.3
Test pattern used
to illustrate the
characteristics of
the PDFs from
Fig. 5.2.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 323 6/16/2017 2:07:01 PM

 324 Chapter 5 Image Restoration and Reconstruction

a b c
d e f
FIGURE 5.4 Images and histograms resulting from adding Gaussian, Rayleigh, and Erlanga noise to the image in
Fig. 5.3.

visually between the first five images in Fig. 5.4, even though their histograms are significantly different.
The salt-and-pepper appearance of the image in Fig. 5.4(i) is the only one that is visually indicative of
the type of noise causing the degradation.

PERIODIC NOISE
Periodic noise in images typically arises from electrical or electromechanical inter-
ference during image acquisition. This is the only type of spatially dependent noise
we will consider in this chapter. As we will discuss in Section 5.4, periodic noise can
be reduced significantly via frequency domain filtering. For example, consider the
image in Fig. 5.5(a). This image is corrupted by additive (spatial) sinusoidal noise.
The Fourier transform of a pure sinusoid is a pair of conjugate impulses† located at

† Be careful not to confuse the term impulse in the frequency domain with the use of the same term in impulse
noise discussed earlier, which is in the spatial domain.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 324 6/16/2017 2:07:02 PM

 5.2 Noise Models 325

g h i
j k l
FIGURE 5.4 (continued) Images and histograms resulting from adding exponential, uniform, and salt-and-pepper noise
to the image in Fig. 5.3. In the salt-and-pepper histogram, the peaks in the origin (zero intensity) and at the far end
of the scale are shown displaced slightly so that they do not blend with the page background.

the conjugate frequencies of the sine wave (see Table 4.4). Thus, if the amplitude of
a sine wave in the spatial domain is strong enough, we would expect to see in the
spectrum of the image a pair of impulses for each sine wave in the image. As shown
in Fig. 5.5(b), this is indeed the case. Eliminating or reducing these impulses in the
frequency domain will eliminate or reduce the sinusoidal noise in the spatial domain.
We will have much more to say in Section 5.4 about this and other examples of peri-
odic noise.

ESTIMATING NOISE PARAMETERS

The parameters of periodic noise typically are estimated by inspection of the Fourier
spectrum. Periodic noise tends to produce frequency spikes that often can be detect-
ed even by visual analysis. Another approach is to attempt to infer the periodicity

[Link]

DIP4E_GLOBAL_Print_Ready.indb 325 6/16/2017 2:07:02 PM

 326 Chapter 5 Image Restoration and Reconstruction

a b
FIGURE 5.5
(a) Image
corrupted by
additive
sinusoidal noise.
(b) Spectrum
showing two
conjugate
impulses caused
by the sine wave.
(Original
image courtesy of
NASA.)

of noise components directly from the image, but this is possible only in simplis-
tic cases. Automated analysis is possible in situations in which the noise spikes are
either exceptionally pronounced, or when knowledge is available about the general
location of the frequency components of the interference (see Section 5.4).
The parameters of noise PDFs may be known partially from sensor specifications,
but it is often necessary to estimate them for a particular imaging arrangement. If
the imaging system is available, one simple way to study the characteristics of system
noise is to capture a set of “flat” images. For example, in the case of an optical sen-
sor, this is as simple as imaging a solid gray board that is illuminated uniformly. The
resulting images typically are good indicators of system noise.
When only images already generated by a sensor are available, it is often possible
to estimate the parameters of the PDF from small patches of reasonably constant
background intensity. For example, the vertical strips shown in Fig. 5.6 were cropped
from the Gaussian, Rayleigh, and uniform images in Fig. 5.4. The histograms shown
were calculated using image data from these small strips. The histograms in Fig. 5.4
that correspond to the histograms in Fig. 5.6 are the ones in the middle of the group
of three in Figs. 5.4(d), (e), and (k).We see that the shapes of these histograms cor-
respond quite closely to the shapes of the corresponding histograms in Fig. 5.6. Their
heights are different due to scaling, but the shapes are unmistakably similar.
The simplest use of the data from the image strips is for calculating the mean and
variance of intensity levels. Consider a strip (subimage) denoted by S, and let pS (zi ),
i = 0, 1, 2, … , L − 1, denote the probability estimates (normalized histogram values)
of the intensities of the pixels in S, where L is the number of possible intensities in
the entire image (e.g., 256 for an 8-bit image). As in Eqs. (2-69) and (2-70), we esti-
mate the mean and variance of the pixel values in S as follows:
L −1
z= ∑
i=0
zi pS (zi ) (5-19)

and

[Link]

DIP4E_GLOBAL_Print_Ready.indb 326 6/16/2017 2:07:03 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 327

a b c
FIGURE 5.6 Histograms computed using small strips (shown as inserts) from (a) the Gaussian, (b) the Rayleigh, and
(c) the uniform noisy images in Fig. 5.4.

L −1
s2 = ∑
i=0
(zi − z )2 pS (zi ) (5-20)

The shape of the histogram identifies the closest PDF match. If the shape is approxi-
mately Gaussian, then the mean and variance are all we need because the Gaussian
PDF is specified completely by these two parameters. For the other shapes discussed
earlier, we use the mean and variance to solve for the parameters a and b. Impulse
noise is handled differently because the estimate needed is of the actual probability
of occurrence of white and black pixels. Obtaining this estimate requires that both
black and white pixels be visible, so a mid-gray, relatively constant area is needed in
the image in order to be able to compute a meaningful histogram of the noise. The
heights of the peaks corresponding to black and white pixels are the estimates of Pa
and Pb in Eq. (5-16).

5.3 RESTORATION IN THE PRESENCE OF NOISE ONLY—SPATIAL

5.3
FILTERING
When an image is degraded only by additive noise, Eqs. (5-1) and (5-2) become

g( x, y) = f ( x, y) + h( x, y) (5-21)

and

G(u, v) = F (u, v) + N (u, v) (5-22)

The noise terms generally are unknown, so subtracting them from g( x, y) [G(u, v)]
to obtain f ( x, y) [ F (u, v)] typically is not an option. In the case of periodic noise,

[Link]

DIP4E_GLOBAL_Print_Ready.indb 327 6/16/2017 2:07:03 PM

 328 Chapter 5 Image Restoration and Reconstruction

sometimes it is possible to estimate N(u, v) from the spectrum of G(u, v), as noted
in Section 5.2. In this case N(u, v) can be subtracted from G(u, v) to obtain an esti-
mate of the original image, but this type of knowledge is the exception, rather than
the rule.
Spatial filtering is the method of choice for estimating f ( x, y) [i.e., denoising
image g( x, y)] in situations when only additive random noise is present. Spatial fil-
tering was discussed in detail in Chapter 3. With the exception of the nature of the
computation performed by a specific filter, the mechanics for implementing all the
filters that follow are exactly as discussed in Sections 3.4 through 3.7.

MEAN FILTERS
In this section, we discuss briefly the noise-reduction capabilities of the spatial filters
introduced in Section 3.5 and develop several other filters whose performance is in
many cases superior to the filters discussed in that section.

Arithmetic Mean Filter

The arithmetic mean filter is the simplest of the mean filters (the arithmetic mean
filter is the same as the box filter we discussed in Chapter 3). Let Sxy represent the
set of coordinates in a rectangular subimage window (neighborhood) of size m × n,
We assume that m and
n are odd integers. The centered on point ( x, y). The arithmetic mean filter computes the average value of
size of a mean filter is
the same as the size of
the corrupted image, g( x, y), in the area defined by Sxy . The value of the restored
neighborhood Sxy; that image f̂ at point ( x, y) is the arithmetic mean computed using the pixels in the
is, m  n.
region defined by Sxy .In other words,

1
fˆ ( x, y) = ∑ g( r , c )
mn (r, c )∈Sxy
(5-23)

where, as in Eq. (2-43), r and c are the row and column coordinates of the pixels
contained in the neighborhood Sxy . This operation can be implemented using a spa-
tial kernel of size m × n in which all coefficients have value 1 mn. A mean filter
smooths local variations in an image, and noise is reduced as a result of blurring.

Geometric Mean Filter

An image restored using a geometric mean filter is given by the expression
1
⎡ ⎤ mn
fˆ ( x, y) = ⎢ ∏ g(r, c)⎥ (5-24)
⎢⎣(r, c )∈Sxy ⎥⎦

where indicates multiplication. Here, each restored pixel is given by the product of
all the pixels in the subimage area, raised to the power 1 mn. As Example 5.2 below
illustrates, a geometric mean filter achieves smoothing comparable to an arithmetic
mean filter, but it tends to lose less image detail in the process.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 328 6/16/2017 2:07:05 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 329

Harmonic Mean Filter

The harmonic mean filtering operation is given by the expression

mn
fˆ ( x, y) = (5-25)
1
∑
( r, c )∈Sxy g (r , c )

The harmonic mean filter works well for salt noise, but fails for pepper noise. It does
well also with other types of noise like Gaussian noise.

Contraharmonic Mean Filter

The contraharmonic mean filter yields a restored image based on the expression

∑
( r, c )∈S
g(r, c)Q+ 1
fˆ ( x, y) = xy
(5-26)
∑
( r, c )∈S
g(r, c)Q
xy

where Q is called the order of the filter. This filter is well suited for reducing or vir-
tually eliminating the effects of salt-and-pepper noise. For positive values of Q, the
filter eliminates pepper noise. For negative values of Q, it eliminates salt noise. It
cannot do both simultaneously. Note that the contraharmonic filter reduces to the
arithmetic mean filter if Q = 0, and to the harmonic mean filter if Q = −1.

EXAMPLE 5.2 : Image denoising using spatial mean filters.

Figure 5.7(a) shows an 8-bit X-ray image of a circuit board, and Fig. 5.7(b) shows the same image, but
corrupted with additive Gaussian noise of zero mean and variance of 400. For this type of image, this is
a significant level of noise. Figures 5.7(c) and (d) show, respectively, the result of filtering the noisy image
with an arithmetic mean filter of size 3 × 3 and a geometric mean filter of the same size. Although both
filters did a reasonable job of attenuating the contribution due to noise, the geometric mean filter did
not blur the image as much as the arithmetic filter. For instance, the connector fingers at the top of the
image are sharper in Fig. 5.7(d) than in (c). The same is true in other parts of the image.
Figure 5.8(a) shows the same circuit image, but corrupted now by pepper noise with probability of
0.1. Similarly, Fig. 5.8(b) shows the image corrupted by salt noise with the same probability. Figure 5.8(c)
shows the result of filtering Fig. 5.8(a) using a contraharmonic mean filter with Q = 1.5, and Fig. 5.8(d)
shows the result of filtering Fig. 5.8(b) with Q = −1.5. Both filters did a good job of reducing the effect of
the noise. The positive-order filter did a better job of cleaning the background, at the expense of slightly
thinning and blurring the dark areas. The opposite was true of the negative order filter.
In general, the arithmetic and geometric mean filters (particularly the latter) are well suited for ran-
dom noise like Gaussian or uniform noise. The contraharmonic filter is well suited for impulse noise, but
it has the disadvantage that it must be known whether the noise is dark or light in order to select the
proper sign for Q. The results of choosing the wrong sign for Q can be disastrous, as Fig. 5.9 shows. Some
of the filters discussed in the following sections eliminate this shortcoming.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 329 6/16/2017 2:07:06 PM

 330 Chapter 5 Image Restoration and Reconstruction

a b
c d
FIGURE 5.7
(a) X-ray image
of circuit board.
(b) Image
corrupted by
additive Gaussian
noise. (c) Result
of filtering with
an arithmetic
mean filter of size
3 × 3. (d) Result
of filtering with a
geometric mean
filter of the same
size. (Original
image courtesy of
Mr. Joseph E.
Pascente, Lixi,
Inc.)

ORDER-STATISTIC FILTERS
We introduced order-statistic filters in Section 3.6. We now expand the discussion
in that section and introduce some additional order-statistic filters. As noted in Sec-
tion 3.6, order-statistic filters are spatial filters whose response is based on ordering
(ranking) the values of the pixels contained in the neighborhood encompassed by
the filter. The ranking result determines the response of the filter.

Median Filter
The best-known order-statistic filter in image processing is the median filter, which,
as its name implies, replaces the value of a pixel by the median of the intensity levels
in a predefined neighborhood of that pixel:

fˆ ( x, y) = median { g(r, c)} (5-27)

( r, c )∈Sxy

where, as before, Sxy is a subimage (neighborhood) centered on point ( x, y). The val-
ue of the pixel at ( x, y) is included in the computation of the median. Median filters

[Link]

DIP4E_GLOBAL_Print_Ready.indb 330 6/16/2017 2:07:06 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 331

a b
c d
FIGURE 5.8
(a) Image
corrupted by
pepper noise with
a probability of
0.1. (b) Image
corrupted by salt
noise with the
same
probability.
(c) Result of
filtering (a) with
a 3 × 3 contra-
harmonic filter
Q = 1.5. (d) Result
of filtering (b)
with Q = −1.5.

a b
FIGURE 5.9
Results of
selecting the
wrong sign in
contraharmonic
filtering.
(a) Result of
filtering Fig. 5.8(a)
with a
contraharmonic
filter of size 3 × 3
and Q = −1.5.
(b) Result of
filtering Fig. 5.8(b)
using Q = 1.5.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 331 6/16/2017 2:07:07 PM

 332 Chapter 5 Image Restoration and Reconstruction

are quite popular because, for certain types of random noise, they provide excellent
noise-reduction capabilities, with considerably less blurring than linear smoothing
filters of similar size. Median filters are particularly effective in the presence of both
bipolar and unipolar impulse noise, as Example 5.3 below shows. Computation of
the median and implementation of this filter are discussed in Section 3.6.

Max and Min Filters

Although the median filter is by far the order-statistic filter most used in image pro-
cessing, it is by no means the only one. The median represents the 50th percentile of
a ranked set of numbers, but you will recall from basic statistics that ranking lends
itself to many other possibilities. For example, using the 100th percentile results in
the so-called max filter, given by

fˆ ( x, y) = max
( r , c ) ∈Sxy
{ g(r, c)} (5-28)

This filter is useful for finding the brightest points in an image or for eroding dark
regions adjacent to bright areas. Also, because pepper noise has very low values, it
is reduced by this filter as a result of the max selection process in the subimage area
Sxy .
The 0th percentile filter is the min filter:

fˆ ( x, y) = min
( r, c )∈Sxy
{ g(r, c)} (5-29)

This filter is useful for finding the darkest points in an image or for eroding light
regions adjacent to dark areas. Also, it reduces salt noise as a result of the min opera-
tion.

Midpoint Filter
The midpoint filter computes the midpoint between the maximum and minimum
values in the area encompassed by the filter:

1⎡ ⎤
fˆ ( x, y) = ⎢ max { g(r, c)} + min { g(r, c)}⎥ (5-30)
2 ⎣(r, c )∈Sxy ( r, c )∈Sxy ⎦

Note that this filter combines order statistics and averaging. It works best for ran-
domly distributed noise, like Gaussian or uniform noise.

Alpha-Trimmed Mean Filter

Suppose that we delete the d 2 lowest and the d 2 highest intensity values of g(r, c)
in the neighborhood Sxy . Let gR (r, c) represent the remaining mn − d pixels in Sxy .
A filter formed by averaging these remaining pixels is called an alpha-trimmed mean
filter. The form of this filter is

[Link]

DIP4E_GLOBAL_Print_Ready.indb 332 6/16/2017 2:07:08 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 333

1
fˆ ( x, y) = ∑ g R (r , c )
mn − d (r, c )∈Sxy
(5-31)

where the value of d can range from 0 to mn − 1. When d = 0 the alpha-trimmed fil-
ter reduces to the arithmetic mean filter discussed earlier. If we choose d = mn − 1,
the filter becomes a median filter. For other values of d, the alpha-trimmed filter is
useful in situations involving multiple types of noise, such as a combination of salt-
and-pepper and Gaussian noise.

EXAMPLE 5.3 : Image denoising using order-statistic filters.

Figure 5.10(a) shows the circuit board image corrupted by salt-and-pepper noise with probabilities
Ps = Pp = 0.1. Figure 5.10(b) shows the result of median filtering with a filter of size 3 × 3. The improve-
ment over Fig. 5.10(a) is significant, but several noise points still are visible. A second pass [on the im-
age in Fig. 5.10(b)] with the median filter removed most of these points, leaving only few, barely visible
noise points. These were removed with a third pass of the filter. These results are good examples of the
power of median filtering in handling impulse-like additive noise. Keep in mind that repeated passes
of a median filter will blur the image, so it is desirable to keep the number of passes as low as possible.
Figure 5.11(a) shows the result of applying the max filter to the pepper noise image of Fig. 5.8(a). The
filter did a reasonable job of removing the pepper noise, but we note that it also removed (set to a light
intensity level) some dark pixels from the borders of the dark objects. Figure 5.11(b) shows the result
of applying the min filter to the image in Fig. 5.8(b). In this case, the min filter did a better job than the
max filter on noise removal, but it removed some white points around the border of light objects. These
made the light objects smaller and some of the dark objects larger (like the connector fingers in the top
of the image) because white points around these objects were set to a dark level.
The alpha-trimmed filter is illustrated next. Figure 5.12(a) shows the circuit board image corrupted
this time by additive, uniform noise of variance 800 and zero mean. This is a high level of noise corrup-
tion that is made worse by further addition of salt-and-pepper noise with Ps = Pp = 0.1, as Fig. 5.12(b)
shows. The high level of noise in this image warrants use of larger filters. Figures 5.12(c) through (f) show
the results, respectively, obtained using arithmetic mean, geometric mean, median, and alpha-trimmed
mean (with d = 6) filters of size 5 × 5. As expected, the arithmetic and geometric mean filters (especially
the latter) did not do well because of the presence of impulse noise. The median and alpha-trimmed
filters performed much better, with the alpha-trimmed filter giving slightly better noise reduction. For
example, note in Fig. 5.12(f) that the fourth connector finger from the top left is slightly smoother in
the alpha-trimmed result. This is not unexpected because, for a high value of d, the alpha-trimmed filter
approaches the performance of the median filter, but still retains some smoothing capabilities.

ADAPTIVE FILTERS
Once selected, the filters discussed thus far are applied to an image without regard
for how image characteristics vary from one point to another. In this section, we
take a look at two adaptive filters whose behavior changes based on statistical char-
acteristics of the image inside the filter region defined by the m × n rectangular
neighborhood Sxy . As the following discussion shows, adaptive filters are capable
of performance superior to that of the filters discussed thus far. The price paid for

[Link]

DIP4E_GLOBAL_Print_Ready.indb 333 6/16/2017 2:07:09 PM

 334 Chapter 5 Image Restoration and Reconstruction

a b
c d
FIGURE 5.10
(a) Image
corrupted by salt-
and- pepper noise
with probabilities
Ps = Pp = 0.1.
(b) Result of one
pass with a medi-
an filter of size
3 × 3. (c) Result
of processing (b)
with this filter.
(d) Result of
processing (c)
with the same
filter.

a b
FIGURE 5.11
(a) Result of
filtering Fig. 5.8(a)
with a max filter
of size 3 × 3.
(b) Result of
filtering Fig. 5.8(b)
with a min filter of
the same size.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 334 6/16/2017 2:07:09 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 335

a b
c d
e f
FIGURE 5.12
(a) Image
corrupted by
additive uniform
noise. (b) Image
additionally
corrupted by
additive salt-and-
pepper noise.
(c)-(f) Image (b)
filtered with a
5 × 5:
(c) arithmetic
mean filter;
(d) geometric
mean filter;
(e) median filter;
(f) alpha-trimmed
mean filter, with
d = 6.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 335 6/16/2017 2:07:10 PM

 336 Chapter 5 Image Restoration and Reconstruction

improved filtering power is an increase in filter complexity. Keep in mind that we

still are dealing with the case in which the degraded image is equal to the original
image plus noise. No other types of degradations are being considered yet.

Adaptive, Local Noise Reduction Filter

The simplest statistical measures of a random variable are its mean and variance.
These are reasonable parameters on which to base an adaptive filter because they
are quantities closely related to the appearance of an image. The mean gives a mea-
sure of average intensity in the region over which the mean is computed, and the
variance gives a measure of image contrast in that region.
Our filter is to operate on a neighborhood, Sxy , centered on coordinates ( x, y).
The response of the filter at ( x, y) is to be based on the following quantities: g( x, y),
the value of the noisy image at ( x, y); sh2 , the variance of the noise; zSxy , the local
average intensity of the pixels in Sxy ; and sS2xy, the local variance of the intensities of
pixels in Sxy . We want the behavior of the filter to be as follows:

1. If sh2 is zero, the filter should return simply the value of g at ( x, y). This is the
trivial, zero-noise case in which g is equal to f at ( x, y).
2. If the local variance sS2xy is high relative to sh2 , the filter should return a value
close to g at ( x, y). A high local variance typically is associated with edges, and
these should be preserved.
3. If the two variances are equal, we want the filter to return the arithmetic mean
value of the pixels in Sxy . This condition occurs when the local area has the same
properties as the overall image, and local noise is to be reduced by averaging.

An adaptive expression for obtaining fˆ ( x, y) based on these assumptions may be

written as
sh2
fˆ ( x, y) = g( x, y) − 2 ⎡ g( x, y) − zSxy ⎤ (5-32)
s ⎣ Sxy
⎦

The only quantity that needs to be known a priori is sh2 , the variance of the noise
corrupting image f ( x, y). This is a constant that can be estimated from sample noisy
images using Eq. (3-26). The other parameters are computed from the pixels in
neighborhood Sxy using Eqs. (3-27) and (3-28).
An assumption in Eq. (5-32) is that the ratio of the two variances does not exceed 1,
which implies that sh2 ≤ sS2xy . The noise in our model is additive and position indepen-
dent, so this is a reasonable assumption to make because Sxy is a subset of g( x, y).
However, we seldom have exact knowledge of sh2 . Therefore, it is possible for this
condition to be violated in practice. For that reason, a test should be built into an
implementation of Eq. (5-32) so that the ratio is set to 1 if the condition sh2 > sS2xy
occurs. This makes this filter nonlinear. However, it prevents nonsensical results (i.e.,
negative intensity levels, depending on the value of zSxy ) due to a potential lack of
knowledge about the variance of the image noise. Another approach is to allow the
negative values to occur, and then rescale the intensity values at the end. The result
then would be a loss of dynamic range in the image.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 336 6/16/2017 2:07:12 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 337

EXAMPLE 5.4 : Image denoising using adaptive, local noise-reduction filtering.

Figure 5.13(a) shows the circuit-board image, corrupted this time by additive Gaussian noise of zero
mean and a variance of 1000. This is a significant level of noise corruption, but it makes an ideal test bed
on which to compare relative filter performance. Figure 5.13(b) is the result of processing the noisy im-
age with an arithmetic mean filter of size 7 × 7. The noise was smoothed out, but at the cost of significant
blurring. Similar comments apply to Fig. 5.13(c), which shows the result of processing the noisy image
with a geometric mean filter, also of size 7 × 7. The differences between these two filtered images are
analogous to those we discussed in Example 5.2; only the degree of blurring is different.
Figure 5.13(d) shows the result of using the adaptive filter of Eq. (5-32) with sh2 = 1000. The improve-
ments in this result compared with the two previous filters are significant. In terms of overall noise
reduction, the adaptive filter achieved results similar to the arithmetic and geometric mean filters. How-
ever, the image filtered with the adaptive filter is much sharper. For example, the connector fingers at the
top of the image are significantly sharper in Fig. 5.13(d). Other features, such as holes and the eight legs
of the dark component on the lower left-hand side of the image, are much clearer in Fig. 5.13(d).These
results are typical of what can be achieved with an adaptive filter. As mentioned earlier, the price paid
for the improved performance is additional filter complexity.

a b
c d
FIGURE 5.13
(a) Image
corrupted by
additive
Gaussian noise of
zero mean and a
variance of 1000.
(b) Result of
arithmetic mean
filtering.
(c) Result of
geometric mean
filtering.
(d) Result of
adaptive noise-
reduction filtering.
All filters used
were of size 7 × 7.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 337 6/16/2017 2:07:12 PM

 338 Chapter 5 Image Restoration and Reconstruction

The preceding results used a value for sh2 that matched the variance of the noise exactly. If this
quantity is not known, and the estimate used is too low, the algorithm will return an image that closely
resembles the original because the corrections will be smaller than they should be. Estimates that are
too high will cause the ratio of the variances to be clipped at 1.0, and the algorithm will subtract the
mean from the image more frequently than it would normally. If negative values are allowed and the
image is rescaled at the end, the result will be a loss of dynamic range, as mentioned previously.

Adaptive Median Filter

The median filter in Eq. (5-27) performs well if the spatial density of the salt-and-
pepper noise is low (as a rule of thumb, Ps and Pp less than 0.2). We show in the fol-
lowing discussion that adaptive median filtering can handle noise with probabilities
larger than these. An additional benefit of the adaptive median filter is that it seeks
to preserve detail while simultaneously smoothing non-impulse noise, something
that the “traditional” median filter does not do. As in all the filters discussed in the
preceding sections, the adaptive median filter also works in a rectangular neighbor-
hood, Sxy . Unlike those filters, however, the adaptive median filter changes (increas-
es) the size of Sxy during filtering, depending on certain conditions to be listed short-
ly. Keep in mind that the output of the filter is a single value used to replace the
value of the pixel at ( x, y), the point on which region Sxy is centered at a given time.
We use the following notation:

zmin = minimum intensity value in Sxy

zmax = maximum intensity value in Sxy
zmed = median of intensity values in Sxy
zxy = intensity at coordinates ( x, y)
Smax = maximum allowed size of Sxy

The adaptive median-filtering algorithm uses two processing levels, denoted level A
and level B, at each point ( x, y) :

Level A : If zmin < zmed < zmax , go to Level B

Else, increase the size of Sxy
If Sxy ≤ Smax , repeat level A
Else, output zmed .
Level B : If zmin < zxy < zmax , output zxy
Else output zmed .

where Sxy and Smax are odd, positive integers greater than 1. Another option in the
last step of level A is to output zxy instead of zmed . This produces a slightly less
blurred result, but can fail to detect salt (pepper) noise embedded in a constant
background having the same value as pepper (salt) noise.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 338 6/16/2017 2:07:13 PM

 5.3 Restoration in the Presence of Noise Only—Spatial Filtering 339

This algorithm has three principal objectives: to remove salt-and-pepper (impulse)

noise, to provide smoothing of other noise that may not be impulsive, and to reduce
distortion, such as excessive thinning or thickening of object boundaries. The values
zmin and zmax are considered statistically by the algorithm to be “impulse-like” noise
components in region Sxy , even if these are not the lowest and highest possible pixel
values in the image.
With these observations in mind, we see that the purpose of level A is to deter-
mine if the median filter output, zmed , is an impulse (salt or pepper) or not. If the
condition zmin < zmed < zmax holds, then zmed cannot be an impulse for the reason
mentioned in the previous paragraph. In this case, we go to level B and test to see
if the point in the center of the neighborhood is itself an impulse (recall that ( x, y)
is the location of the point being processed, and zxy is its intensity). If the condition
zmin < zxy < zmax is true, then the pixel at zxy cannot be the intensity of an impulse for
the same reason that zmed was not. In this case, the algorithm outputs the unchanged
pixel value, zxy . By not changing these “intermediate-level” points, distortion is
reduced in the filtered image. If the condition zmin < zxy < zmax is false, then either
zxy = zmin or zxy = zmax . In either case, the value of the pixel is an extreme value and
the algorithm outputs the median value, zmed , which we know from level A is not a
noise impulse. The last step is what the standard median filter does. The problem is
that the standard median filter replaces every point in the image by the median of
the corresponding neighborhood. This causes unnecessary loss of detail.
Continuing with the explanation, suppose that level A does find an impulse (i.e.,
it fails the test that would cause it to branch to level B). The algorithm then increas-
es the size of the neighborhood and repeats level A. This looping continues until
the algorithm either finds a median value that is not an impulse (and branches to
stage B), or the maximum neighborhood size is reached. If the maximum size is
reached, the algorithm returns the value of zmed . Note that there is no guarantee
that this value is not an impulse. The smaller the noise probabilities Pa and/or Pb are,
or the larger Smax is allowed to be, the less likely it is that a premature exit will occur.
This is plausible. As the density of the noise impulses increases, it stands to reason
that we would need a larger window to “clean up” the noise spikes.
Every time the algorithm outputs a value, the center of neighborhood Sxy is
moved to the next location in the image. The algorithm then is reinitialized and
applied to the pixels in the new region encompassed by the neighborhood. As indi-
cated in Problem 3.37, the median value can be updated iteratively from one loca-
tion to the next, thus reducing computational load.

EXAMPLE 5.5 : Image denoising using adaptive median filtering.

Figure 5.14(a) shows the circuit-board image corrupted by salt-and-pepper noise with probabilities
Ps = Pp = 0.25, which is 2.5 times the noise level used in Fig. 5.10(a). Here the noise level is high enough
to obscure most of the detail in the image. As a basis for comparison, the image was filtered first using a
7 × 7 median filter, the smallest filter required to remove most visible traces of impulse noise in this case.
Figure 5.14(b) shows the result. Although the noise was effectively removed, the filter caused significant

[Link]

DIP4E_GLOBAL_Print_Ready.indb 339 6/16/2017 2:07:15 PM

 340 Chapter 5 Image Restoration and Reconstruction

a b c
FIGURE 5.14 (a) Image corrupted by salt-and-pepper noise with probabilities Ps = Pp = 0.25. (b) Result of filtering
with a 7 × 7 median filter. (c) Result of adaptive median filtering with Smax = 7.

loss of detail in the image. For instance, some of the connector fingers at the top of the image appear
distorted or broken. Other image details are similarly distorted.
Figure 5.14(c) shows the result of using the adaptive median filter with Smax = 7. Noise removal
performance was similar to the median filter. However, the adaptive filter did a much better job of pre-
serving sharpness and detail. The connector fingers are less distorted, and some other features that were
either obscured or distorted beyond recognition by the median filter appear sharper and better defined
in Fig. 5.14(c). Two notable examples are the feed-through small white holes throughout the board, and
the dark component with eight legs in the bottom, left quadrant of the image.
Considering the high level of noise in Fig. 5.14(a), the adaptive algorithm performed quite well. The
choice of maximum allowed size for Sxy depends on the application, but a reasonable starting value can
be estimated by experimenting with various sizes of the standard median filter first. This will establish a
visual baseline regarding expectations on the performance of the adaptive algorithm.

5.4 PERIODIC NOISE REDUCTION USING FREQUENCY DOMAIN

5.4

FILTERING
Periodic noise can be analyzed and filtered quite effectively using frequency domain
techniques. The basic idea is that periodic noise appears as concentrated bursts of
energy in the Fourier transform, at locations corresponding to the frequencies of
the periodic interference. The approach is to use a selective filter (see Section 4.10)
to isolate the noise. The three types of selective filters (bandreject, bandpass, and
notch) were discussed in detail in Section 4.10. There is no difference between how
these filters were used in Chapter 4, and the way they are used for image restora-
tion. In restoration of images corrupted by periodic interference, the tool of choice
is a notch filter. In the following discussion we will expand on the notch filtering
approach introduced in Section 4.10, and also develop a more powerful optimum
notch filtering method.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 340 6/16/2017 2:07:16 PM

 5.4 Periodic Noise Reduction Using Frequency Domain Filtering 341

MORE ON NOTCH FILTERING

As explained in Section 4.10, notch reject filter transfer functions are constructed
as products of highpass filter transfer functions whose centers have been translated
to the centers of the notches. The general form of a notch filter transfer function is
Q
H NR (u, v) = ∏ H k (u, v) H − k (u, v) (5-33)
k =1

where H k (u, v) and H − k (u, v) are highpass filter transfer functions whose centers
are at (uk , vk ) and (− uk , − vk ), respectively.† These centers are specified with respect
to the center of the frequency rectangle, [ floor( M 2), floor( N 2)], where, as usual,
M and N are the number of rows and columns in the input image. Thus, the distance
computations for the filter transfer functions are given by
1/ 2
Dk (u, v) = ⎡⎣(u − M 2 − uk )2 + (v − N 2 − vk )2 ⎤⎦ (5-34)

and
1/ 2
D− k (u, v) = ⎡⎣(u − M 2 + uk )2 + (v − N 2 + vk )2 ⎤⎦ (5-35)

For example, the following is a Butterworth notch reject filter transfer function of
order n with three notch pairs:
3 ⎡ ⎤⎡ ⎤
1 1
H NR (u, v) = ∏ ⎢ ⎥ ⎢ ⎥ (5-36)
k =1 ⎢ 1 + [ D0 k Dk (u, v)] ⎥ ⎢ 1 + [ D0 k D− k (u, v)] ⎥
n n
⎣ ⎦⎣ ⎦
Because notches are specified as symmetric pairs, the constant D0 k is the same for
each pair. However, this constant can be different from one pair to another. Other
notch reject filter functions are constructed in the same manner, depending on the
highpass filter function chosen. As explained in Section 4.10, a notch pass filter
transfer function is obtained from a notch reject function using the expression

H NP (u, v) = 1 − H NR (u, v) (5-37)

where H NP (u, v) is the transfer function of the notch pass filter corresponding to
the notch reject filter with transfer function H NR (u, v). Figure 5.15 shows perspec-
tive plots of the transfer functions of ideal, Gaussian, and Butterworth notch reject
filters with one notch pair. As we discussed in Chapter 4, we see again that the shape
of the Butterworth transfer function represents a transition between the sharpness
of the ideal function and the broad, smooth shape of the Gaussian transfer function.
As we show in the second part of the following example, we are not limited to
notch filter transfer functions of the form just discussed. We can construct notch

† Remember, frequency domain transfer functions are symmetric about the center of the frequency rectangle, so
the notches are specified as symmetric pairs. Also, recall from Section 4.10 that we use unpadded images when
working with notch filters in order to simplify the specification of notch locations.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 341 6/16/2017 2:07:17 PM

 342 Chapter 5 Image Restoration and Reconstruction

H ( u, v) H ( u, v) H ( u, v)

u Ideal u Gaussian u Butterworth

v v v

a b c
FIGURE 5.15 Perspective plots of (a) ideal, (b) Gaussian, and (c) Butterworth notch reject filter transfer functions.

filters of arbitrary shapes, provided that they are zero-phase-shift functions, as

defined in Section 4.7.

EXAMPLE 5.6 : Image denoising (interference reduction) using notch filtering.

Figure 5.16(a) is the same as Fig. 2.45(a), which we used in Section 2.6 to introduce the concept of filter-
ing in the frequency domain. We now look in more detail at the process of denoising this image, which is
corrupted by a single, 2-D additive sine wave. You know from Table 4.4 that the Fourier transform of a
pure sine wave is a pair of complex, conjugate impulses, so we would expect the spectrum to have a pair
of bright dots at the frequencies of the sine wave. As Fig. 5.16(b) shows, this is indeed is the case. Because
we can determine the location of these impulses accurately, eliminating them is a simple task, consisting
of using a notch filter transfer function whose notches coincide with the location of the impulses.
Figure 5.16(c) shows an ideal notch reject filter transfer function, which is an array of 1's (shown in
white) and two small circular regions of 0's (shown in black). Figure 5.16(d) shows the result of filtering
the noisy image this transfer function. The sinusoidal noise was virtually eliminated, and a number of
details that were previously obscured by the interference are clearly visible in the filtered image (see, for
example, the thin fiducial marks and the fine detail in the terrain and rock formations). As we showed
in Example 4.25, obtaining an image of the interference pattern is straightforward. We simply turn the
reject filter into a pass filter by subtracting it from 1, and filter the input image with it. Figure 5.17 shows
the result.
Figure 5.18(a) shows the same image as Fig. 4.50(a), but covering a larger area (the interference
pattern is the same). When we discussed lowpass filtering of that image in Chapter 4, we indicated that
there were better ways to reduce the effect of the scan lines. The notch filtering approach that follows
reduces the scan lines significantly, without introducing blurring. Unless blurring is desirable for reasons
we discussed in Section 4.9, notch filtering generally gives much better results.
Just by looking at the nearly horizontal lines of the noise pattern in Fig. 5.18(a), we expect its con-
tribution in the frequency domain to be concentrated along the vertical axis of the DFT. However,
the noise is not dominant enough to have a clear pattern along this axis, as is evident in the spectrum
shown in Fig. 5.18(b). The approach to follow in cases like this is to use a narrow, rectangular notch filter
function that extends along the vertical axis, and thus eliminates all components of the interference
along that axis. We do not filter near the origin to avoid eliminating the dc term and low frequencies,

[Link]

DIP4E_GLOBAL_Print_Ready.indb 342 6/16/2017 2:07:20 PM

 5.4 Periodic Noise Reduction Using Frequency Domain Filtering 343

a b
c d
FIGURE 5.16
(a) Image cor-
rupted by sinusoi-
dal interference.
(b) Spectrum
showing the
bursts of energy
caused by the
interference. (The
bursts were
enlarged for
display purposes.)
(c) Notch filter
(the radius of the
circles is 2 pixels)
used to eliminate
the energy bursts.
(The thin borders
are not part of the
data.)
(d) Result of
notch reject
filtering.
(Original
image courtesy of
NASA.)

which, as you know from Chapter 4, are responsible for the intensity differences between smooth areas.
Figure 5.18(c) shows the filter transfer function we used, and Fig. 5.18(d) shows the filtered result. Most
of the fine scan lines were eliminated or significantly attenuated. In order to get an image of the noise
pattern, we proceed as before by converting the reject filter into a pass filter, and then filtering the input
image with it. Figure 5.19 shows the result.

FIGURE 5.17
Sinusoidal
pattern extracted
from the DFT
of Fig. 5.16(a)
using a notch pass
filter.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 343 6/16/2017 2:07:22 PM

 344 Chapter 5 Image Restoration and Reconstruction

a b
c d
FIGURE 5.18
(a) Satellite image
of Florida and the
Gulf of Mexico.
(Note horizontal
sensor scan lines.)
(b) Spectrum of
(a). (c) Notch
reject filter
transfer
function. (The
thin black border
is not part of the
data.) (d) Filtered
image. (Original
image courtesy of
NOAA.)

FIGURE 5.19
Noise pattern
extracted from
Fig. 5.18(a) by
notch pass
filtering.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 344 6/16/2017 2:07:22 PM

 5.4 Periodic Noise Reduction Using Frequency Domain Filtering 345

OPTIMUM NOTCH FILTERING

In the examples of notch filtering given thus far, the interference patterns have been
simple to identify and characterize in the frequency domain, leading to the specifica-
tion of notch filter transfer functions that also are simple to define heuristically.
When several interference components are present, heuristic specifications of
filter transfer functions are not always acceptable because they may remove too
much image information in the filtering process (a highly undesirable feature when
images are unique and/or expensive to acquire). In addition, the interference com-
ponents generally are not single-frequency bursts. Instead, they tend to have broad
skirts that carry information about the interference pattern. These skirts are not
always easily detectable from the normal transform background. Alternative filter-
ing methods that reduce the effect of these degradations are quite useful in practice.
The method discussed next is optimum, in the sense that it minimizes local variances
of the restored estimate fˆ ( x, y).
The procedure consists of first isolating the principal contributions of the interfer-
ence pattern and then subtracting a variable, weighted portion of the pattern from
the corrupted image. Although we develop the procedure in the context of a specific
application, the basic approach is general and can be applied to other restoration
tasks in which multiple periodic interference is a problem.
We begin by extracting the principal frequency components of the interfer-
ence pattern. As before, we do this by placing a notch pass filter transfer function,
H NP (u, v), at the location of each spike. If the filter is constructed to pass only com-
ponents associated with the interference pattern, then the Fourier transform of the
interference noise pattern is given by the expression

N (u, v) = H NP (u, v) G(u, v) (5-38)

where, as usual, G(u, v) is the DFT of the corrupted image.

Specifying H NP (u, v) requires considerable judgment about what is or is not an
interference spike. For this reason, the notch pass filter generally is constructed inter-
actively by observing the spectrum of G(u, v) on a display. After a particular filter
function has been selected, the corresponding noise pattern in the spatial domain is
obtained using the familiar expression

h( x, y) = ᑣ−1 {H NP (u, v) G(u, v)} (5-39)

Because the corrupted image is assumed to be formed by the addition of the uncor-
rupted image f ( x, y) and the interference, h( x, y), if the latter were known com-
pletely, subtracting the pattern from g( x, y) to obtain f ( x, y) would be a simple mat-
ter. The problem, of course, is that this filtering procedure usually yields only an
approximation of the true noise pattern. The effect of incomplete components not
present in the estimate of h( x, y) can be minimized by subtracting from g( x, y) a
weighted portion of h( x, y) to obtain an estimate of f ( x, y) :

fˆ ( x, y) = g( x, y) − w( x, y)h( x, y) (5-40)

[Link]

DIP4E_GLOBAL_Print_Ready.indb 345 6/16/2017 2:07:23 PM

 346 Chapter 5 Image Restoration and Reconstruction

where, as before, fˆ ( x, y) is the estimate of f ( x, y) and w( x, y) is to be determined.

This function is called a weighting or modulation function, and the objective of the
procedure is to select w( x, y) so that the result is optimized in some meaningful way.
One approach is to select w( x, y) so that the variance of fˆ ( x, y) is minimized over a
specified neighborhood of every point ( x, y).
Consider a neighborhood Sxy of (odd) size m × n, centered on ( x, y). The “local”
variance of fˆ ( x, y) at point ( x, y) can be estimated using the samples in Sxy , as fol-
lows:
_
1
∑ [ f (r , c ) − f ]
^ ^ 2
s 2 ( x, y) = (5-41)
mn (r , c )∈Sxy

where f̂ is the average value of f̂ in neighborhood Sxy ; that is,

1
fˆ = ∑ fˆ (r, c)
mn (r, c )∈Sxy
(5-42)

Points on or near the edge of the image can be treated by considering partial neigh-
borhoods or by padding the border with 0's.
Substituting Eq. (5-40) into Eq. (5-41) we obtain

2
1 ⎧ ____
⎫
s 2 ( x, y) = ∑ ⎨[ g(r, c) − w(r, c)h(r, c)] − [ g − wh]⎬
mn (r, c )∈Sxy ⎩ ⎭
(5-43)

____
where g and wh denote the average values of g and of the product wh in neighbor-
hood Sxy , respectively.
If we assume that w is approximately constant in Sxy we can replace w(r, c) by
the value of w at the center of the neighborhood:

w(r, c) = w( x, y) (5-44)
__
Because w( x, y) is assumed to be constant in Sxy , it follows that w=w( x, y) and,
therefore, that
____ __
wh = w( x, y) h (5-45)

in Sxy , where h is the average value of h in the neighborhood. Using these approxi-
mations, Eq. (5-43) becomes

2
1 ⎧ __ ⎫
s 2 ( x, y) = ∑ ⎨[ g(r, c) − w( x, y)h(r, c)] − [ g − w( x, y) h ]⎬
mn (r, c )∈Sxy ⎩ ⎭
(5-46)

[Link]

DIP4E_GLOBAL_Print_Ready.indb 346 6/16/2017 2:07:26 PM

 5.4 Periodic Noise Reduction Using Frequency Domain Filtering 347

To minimize s 2 ( x, y) with respect to w( x, y) we solve

∂s 2 ( x, y)
=0 (5-47)
∂ w( x, y)

for w( x, y). The result is (see Problem 5.17):

____
gh − g h
w( x, y) = ___
__
(5-48)
h − h2
2

To obtain the value of the restored image at point ( x, y) we use this equation to com-
pute w( x, y) and then substitute it into Eq. (5-40). To obtain the complete restored
image, we perform this procedure at every point in the noisy image, g.

EXAMPLE 5.7 : Denoising (interference removal) using optimum notch filtering.

Figure 5.20(a) shows a digital image of the Martian terrain taken by the Mariner 6 spacecraft. The image
is corrupted by a semi-periodic interference pattern that is considerably more complex (and much more
subtle) than those we have studied thus far. The Fourier spectrum of the image, shown in Fig. 5.20(b),
has a number of “starlike” bursts of energy caused by the interference. As expected, these components
are more difficult to detect than those we have seen before. Figure 5.21 shows the spectrum again, but
without centering. This image offers a somewhat clearer view of the interference components because
the more prominent dc term and low frequencies are “out of way,” in the top left of the spectrum.
Figure 5.22(a) shows the spectrum components that, in the judgement of an experienced image ana-
lyst, are associated with the interference. Applying a notch pass filter to these components and using
Eq. (5-39) yielded the spatial noise pattern, h( x, y), shown in Fig. 5.22(b). Note the similarity between
this pattern and the structure of the noise in Fig. 5.20(a).

a b
FIGURE 5.20
(a) Image of the
Martian
terrain taken by
Mariner 6.
(b) Fourier
spectrum showing
periodic
interference.
(Courtesy of
NASA.)

[Link]

DIP4E_GLOBAL_Print_Ready.indb 347 6/16/2017 2:07:27 PM

 348 Chapter 5 Image Restoration and Reconstruction

FIGURE 5.21
Uncentered
Fourier spectrum
of the image
in Fig. 5.20(a).
(Courtesy of
NASA.)

Finally, Fig. 5.23 shows the restored image, obtained using Eq. (5-40) with the interference pattern just
discussed. Function w( x, y) was computed using the procedure explained in the preceding paragraphs.
As you can see, the periodic interference was virtually eliminated from the noisy image in Fig. 5.20(a).

5.5 LINEAR, POSITION-INVARIANT DEGRADATIONS

5.5

The input-output relationship in Fig. 5.1 before the restoration stage is expressed as

g( x, y) = Ᏼ [ f ( x, y)] + h( x, y) (5-49)

For the moment, let us assume that h( x, y) = 0 so that g( x, y) = Ᏼ [ f ( x, y)]. Based on

the discussion in Section 2.6, Ᏼ is linear if

a b
FIGURE 5.22
(a) Fourier spec-
trum of N(u, v),
and
(b) corresponding
spatial noise
interference
pattern, h( x, y).
(Courtesy of
NASA.)

[Link]

DIP4E_GLOBAL_Print_Ready.indb 348 6/16/2017 2:07:28 PM

 5.5 Linear, Position-Invariant Degradations 349
FIGURE 5.23
Restored image.
(Courtesy of
NASA.)

Ᏼ [ af1 ( x, y) + bf2 ( x, y)] = aᏴ [ f1 ( x, y)] + bᏴ [ f2 ( x, y)] (5-50)

where a and b are scalars and f1 ( x, y) and f2 ( x, y) are any two input images.
If a = b = 1, Eq. (5-50) becomes

Ᏼ [ f1 ( x, y) + f2 ( x, y)] = Ᏼ [ f1 ( x, y)] + Ᏼ [ f2 ( x, y)] (5-51)

which is called the property of additivity. This property says that, if Ᏼ is a linear
operator, the response to a sum of two inputs is equal to the sum of the two responses.
With f2 ( x, y) = 0, Eq. (5-50) becomes

Ᏼ [ af1 ( x, y)] = aᏴ [ f1 ( x, y)] (5-52)

which is called the property of homogeneity. It says that the response to a constant
multiple of any input is equal to the response to that input multiplied by the same
constant. Thus, a linear operator possesses both the property of additivity and the
property of homogeneity.
An operator having the input-output relationship g( x, y) = Ᏼ [ f ( x, y)] is said to
be position (or space) invariant if

Ᏼ [ f ( x − a, y − b)] = g( x − a, y − b) (5-53)

for any f ( x, y) and any two scalars a and b. This definition indicates that the
response at any point in the image depends only on the value of the input at that
point, not on its position.
Using the sifting property of the 2-D continuous impulse [see Eq. (4-55)], we can
write f ( x, y) as

[Link]

DIP4E_GLOBAL_Print_Ready.indb 349 6/16/2017 2:07:29 PM

 350 Chapter 5 Image Restoration and Reconstruction

 
f ( x, y) = f (a, b) d( x − a, y − b) da db (5-54)
2- 2-
Assuming again that h( x, y) = 0, substituting this equation into Eq. (5-49) yields

⎡   ⎤
g( x, y) = Ᏼ [ f ( x, y)] = Ᏼ ⎢ f (a, b) d( x − a, y − b) da db ⎥ (5-55)
⎢⎣ 2
- 2
- ⎥⎦

If Ᏼ is a linear operator and we extend the additivity property to integrals, then

 
g( x, y) = Ᏼ [ f (a, b) d( x − a, y − b)] da db (5-56)
2- 2-
Because f (a, b) is independent of x and y, and using the homogeneity property, it
follows that
 
g( x, y) = f (a, b)Ᏼ [ d( x − a, y − b)] da db (5-57)
2-  2- 
The term

h( x, a, y, b) = Ᏼ [ d( x − a, y − b)] (5-58)

is called the impulse response of Ᏼ. In other words, if h( x, y) = 0 in Eq. (5-49), then

h( x, a, y, b) is the response of Ᏼ to an impulse at coordinates ( x, y). In optics, the
impulse becomes a point of light and h( x, a, y, b) is commonly referred to as the
point spread function (PSF). This name is based on the fact that all physical optical
systems blur (spread) a point of light to some degree, with the amount of blurring
being determined by the quality of the optical components.
Substituting Eq. (5-58) into Eq. (5-57) we obtain the expression
 
g( x, y) = f (a, b)h( x, a, y, b) da db (5-59)
2- 2-
which is called the superposition (or Fredholm) integral of the first kind. This expres-
sion is a fundamental result that is at the core of linear system theory. It states that
if the response of Ᏼ to an impulse is known, the response to any input f (a, b) can
be calculated using Eq. (5-59). In other words, a linear system Ᏼ is characterized
completely by its impulse response.
If Ᏼ is position invariant, then it follows from Eq. (5-53) that

Ᏼ [ d( x − a, y − b)] = h( x − a, y − b) (5-60)

In this case, Eq. (5-59) reduces to

 
g( x, y) = f (a, b)h( x − a, y − b) da db (5-61)
2- 2-

[Link]

DIP4E_GLOBAL_Print_Ready.indb 350 6/16/2017 2:07:31 PM

 5.5 Linear, Position-Invariant Degradations 351

This expression is the convolution integral introduced for one variable in Eq. (4-24)
and extended to 2-D in Problem 4.19. Equation (5-61) tells us that the output of a
linear, position invariant system to any input, is obtained by convolving the input
and the system’s impulse response.
In the presence of additive noise, the expression of the linear degradation model
[Eq. (5-59)] becomes
 
g( x, y) = f (a, b)h( x, a, y, b) da db + h( x, y) (5-62)
2- 2-
If Ᏼ is position invariant, then this equation becomes

 
g( x, y) = f (a, b)h( x − a, y − b) da db + h( x, y) (5-63)
2- 2-
The values of the noise term h( x, y) are random, and are assumed to be independent
of position. Using the familiar notation for convolution introduced in Chapters 3
and 4, we can write Eq. (5-63) as

g( x, y) = (h 夹 f )( x, y) + h( x, y) (5-64)

or, using the convolution theorem, we write the equivalent result in the frequency
domain as

G(u, v) = H (u, v)F (u, v) + N (u, v) (5-65)

These two expressions agree with Eqs. (5-1) and (5-2). Keep in mind that, for dis-
crete quantities, all products are elementwise products, as defined in Section 2.6.
In summary, the preceding discussion indicates that a linear, spatially invariant
degradation system with additive noise can be modeled in the spatial domain as
the convolution of an image with the system’s degradation (point spread) function,
followed by the addition of noise. Based on the convolution theorem, the same pro-
cess can be expressed in the frequency domain as the product of the transforms of
the image and degradation, followed by the addition of the transform of the noise.
When working in the frequency domain, we make use of an FFT algorithm. Howev-
er, unlike in Chapter 4, we do not use image padding in the implementation of any of
the frequency domain restoration filters discussed in this chapter. The reason is that
in restoration work we usually have access only to degraded images. For padding
to be effective, it would have to be applied to images before they were degraded, a
condition that obviously cannot be met in practice. If we had access to the original
images, then restoration would be a mute point.
Many types of degradations can be approximated by linear, position-invariant
processes. The advantage of this approach is that the extensive tools of linear sys-
tem theory then become available for the solution of image restoration problems.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 351 6/16/2017 2:07:31 PM

 352 Chapter 5 Image Restoration and Reconstruction

Nonlinear and position-dependent techniques, although more general (and usually

more accurate), introduce difficulties that often have no known solution or are very
difficult to solve computationally. This chapter focuses on linear, space-invariant res-
toration techniques. Because degradations are modeled as being the result of convo-
lution, and restoration seeks to find filters that apply the process in reverse, the term
image deconvolution is used frequently to signify linear image restoration. Similarly,
the filters used in the restoration process often are called deconvolution filters.

5.6 ESTIMATING THE DEGRADATION FUNCTION

5.6

There are three principal ways to estimate the degradation function for use in image
restoration: (1) observation, (2) experimentation, and (3) mathematical modeling.
These methods are discussed in the following sections. The process of restoring
an image by using a degradation function that has been estimated by any of these
approaches sometimes is called blind deconvolution, to emphasize the fact that the
true degradation function is seldom known completely.

ESTIMATION BY IMAGE OBSERVATION

Suppose that we are given a degraded image without any knowledge about the degra-
dation function Ᏼ. Based on the assumption that the image was degraded by a lin-
ear, position-invariant process, one way to estimate Ᏼ is to gather information from
the image itself. For example, if the image is blurred, we can look at a small rectan-
gular section of the image containing sample structures, like part of an object and
the background. In order to reduce the effect of noise, we would look for an area in
which the signal content is strong (e.g., an area of high contrast). The next step would
be to process the subimage to arrive at a result that is as unblurred as possible.
Let the observed subimage be denoted by gs ( x, y), and let the processed subimage
(which in reality is our estimate of the original image in that area) be denoted by
fˆs ( x, y). Then, assuming that the effect of noise is negligible because of our choice of
a strong-signal area, it follows from Eq. (5-65) that

Gs (u, v)
H s (u, v) = ⁄
(5-66)
F s (u, v)
From the characteristics of this function, we then deduce the complete degradation
function H(u, v) based on our assumption of position invariance. For example, sup-
pose that a radial plot of H s (u, v) has the approximate shape of a Gaussian curve. We
can use that information to construct a function H(u, v) on a larger scale, but having
the same basic shape. We then use H(u, v) in one of the restoration approaches to
be discussed in the following sections. Clearly, this is a laborious process used only in
very specific circumstances, such as restoring an old photograph of historical value.

ESTIMATION BY EXPERIMENTATION
If equipment similar to the equipment used to acquire the degraded image is avail-
able, it is possible in principle to obtain an accurate estimate of the degradation.
Images similar to the degraded image can be acquired with various system settings

[Link]

DIP4E_GLOBAL_Print_Ready.indb 352 6/16/2017 2:07:32 PM

 5.6 Estimating the Degradation Function 353

a b
FIGURE 5.24
Estimating a
degradation by
impulse
characterization.
(a) An impulse
of light (shown
magnified).
(b) Imaged
(degraded)
impulse.

until they are degraded as closely as possible to the image we wish to restore. Then
the idea is to obtain the impulse response of the degradation by imaging an impulse
(small dot of light) using the same system settings. As noted in Section 5.5, a linear,
space-invariant system is characterized completely by its impulse response.
An impulse is simulated by a bright dot of light, as bright as possible to reduce the
effect of noise to negligible values. Then, recalling that the Fourier transform of an
impulse is a constant, it follows from Eq. (5-65) that

G(u, v)
H (u, v) = (5-67)
A

where, as before, G(u, v) is the Fourier transform of the observed image, and A is a
constant describing the strength of the impulse. Figure 5.24 shows an example.

ESTIMATION BY MODELING
Degradation modeling has been used for many years because of the insight it affords
into the image restoration problem. In some cases, the model can even take into
account environmental conditions that cause degradations. For example, a degrada-
tion model proposed by Hufnagel and Stanley [1964] is based on the physical char-
acteristics of atmospheric turbulence. This model has a familiar form:
2
+ v 2 )5 / 6
H (u, v) = e − k(u (5-68)

where k is a constant that depends on the nature of the turbulence. With the excep-
tion of the 5 6 power in the exponent, this equation has the same form as the Gauss-
ian lowpass filter transfer function discussed in Section 4.8. In fact, the Gaussian
LPF is used sometimes to model mild, uniform blurring. Figure 5.25 shows examples
obtained by simulating blurring an image using Eq. (5-68) with values k = 0.0025

[Link]

DIP4E_GLOBAL_Print_Ready.indb 353 6/16/2017 2:07:33 PM

 354 Chapter 5 Image Restoration and Reconstruction

a b
c d
FIGURE 5.25
Modeling
turbulence.
(a) No visible
turbulence.
(b) Severe
turbulence,
k = 0.0025.
(c) Mild
turbulence,
k = 0.001.
(d) Low
turbulence,
k = 0.00025.
All images are
of size 480 × 480
pixels.
(Original
image courtesy of
NASA.)

(severe turbulence), k = 0.001 (mild turbulence), and k = 0.00025 (low turbulence).

We restore these images using various methods later in this chapter.
Another approach used frequently in modeling is to derive a mathematical model
starting from basic principles. We illustrate this procedure by treating in some detail
the case in which an image has been blurred by uniform linear motion between
the image and the sensor during image acquisition. Suppose that an image f ( x, y)
undergoes planar motion and that x0 (t ) and y0 (t ) are the time-varying components
of motion in the x- and y-directions, respectively. We obtain the total exposure at
any point of the recording medium (say, film or digital memory) by integrating the
instantaneous exposure over the time interval during which the imaging system
shutter is open.
Assuming that shutter opening and closing takes place instantaneously, and that
the optical imaging process is perfect, lets us isolate the effects due to image motion.
Then, if T is the duration of the exposure, it follows that

[Link]

DIP4E_GLOBAL_Print_Ready.indb 354 6/16/2017 2:07:33 PM

 5.6 Estimating the Degradation Function 355

T
g( x, y) = f [ x − x0 (t ), y − y0 (t )] dt (5-69)
20
where g( x, y) is the blurred image.
The continuous Fourier transform of this expression is
 
G(u, v) = g( x, y) e − j 2p(ux + vy ) dx dy (5-70)
2- 2-
Substituting Eq. (5-69) into Eq. (5-70) yields

⎡ T  ⎤
G(u, v) = ⎢ f [ x − x0 (t ), y − y0 (t )] dt ⎥ e − j 2p(ux + vy ) dx dy (5-71)
2-  2-  ⎢⎣ 20 ⎥⎦
Reversing the order of integration results in the expression
T
⎡   ⎤
G(u, v) = ⎢ f [ x − x0 (t ), y − y0 (t )] e − j 2 p(ux + vy ) dx dy ⎥ dt (5-72)
20 ⎢⎣ 2
- 2
- ⎥⎦

The term inside the outer brackets is the Fourier transform of the displaced function
[ ]
f x − x0 (t ), y − y0 (t ) . Using entry 3 in Table 4.4 then yields the expression
T
− j 2 p[ ux0 ( t ) + vy0 ( t )]
G(u, v) = F (u, v)e dt
20 (5-73)
T
− j 2 p[ ux0 ( t ) + vy0 ( t )]
= F (u, v) e dt
20
By defining
T
− j 2 p[ ux0 ( t ) + vy0 ( t )]
H (u, v) = e dt (5-74)
20
we can express Eq. (5-73) in the familiar form

G(u, v) = H (u, v)F (u, v) (5-75)

If the motion variables x0 (t ) and y0 (t ) are known, the transfer function H(u, v) can
be obtained directly from Eq. (5-74). As an illustration, suppose that the image in
question undergoes uniform linear motion in the x-direction only (i.e., y0 (t ) = 0), at
a rate x0 (t ) = at T . When t = T , the image has been displaced by a total distance a.
With y0 (t ) = 0, Eq. (5-74) yields
T T
H (u, v) = e − j 2 pux0 (t ) dt = e − j 2 puat T dt
20 20 (5-76)
T
= sin(pua)e − jpua
pua

[Link]

DIP4E_GLOBAL_Print_Ready.indb 355 6/16/2017 2:07:35 PM

 356 Chapter 5 Image Restoration and Reconstruction

a b
FIGURE 5.26
(a) Original
image. (b) Result
of blurring using
the function in
Eq. (5-77) with
a = b = 0.1 and
T = 1.

If we allow the y-component to vary as well, with the motion given by y0 (t ) = bt T ,

then the degradation function becomes

T
H (u, v) = sin [ p(ua + vb)] e − jp(ua + vb )
(5-77)
p(ua + vb)
To generate a discrete filter transfer function of size M × N , we sample this equation
for u = 0, 1, 2, … , M − 1 and v = 0, 1, 2, … , N − 1.

EXAMPLE 5.8 : Image blurring caused by motion.

Figure 5.26(b) is an image blurred by computing the Fourier transform of the image in Fig. 5.26(a), mul-
tiplying the transform by H(u, v) from Eq. (5-77), and taking the inverse transform. The images are of
size 688 × 688 pixels, and we used a = b = 0.1 and T = 1 in Eq. (5-77). As we will discuss in Sections 5.8
and 5.9, recovery of the original image from its blurred counterpart presents some interesting challenges,
particularly when noise is present in the degraded image. As mentioned at the end of Section 5.5, we
perform all DFT computations without padding.

5.7 INVERSE FILTERING

5.7

The material in this section is our first step in studying restoration of images degrad-
ed by a degradation function Ᏼ, which is given, or is obtained by a method such
as those discussed in the previous section. 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 dividing the transform of the degraded image, G(u, v), by the
degradation transfer function:

G(u, v)
Fˆ (u, v) = (5-78)
H (u, v)

[Link]

DIP4E_GLOBAL_Print_Ready.indb 356 6/16/2017 2:07:37 PM

 5.7 Inverse Filtering 357

The division is elementwise, as defined in Section 2.6 and in connection with Eq.
(5-65). Substituting the right side of Eq. (5-2) for G(u, v) in Eq. (5-78) yields

N (u, v)
Fˆ (u, v) = F (u, v) + (5-79)
H (u, v)
This is an interesting expression. It tells us that, even if we know the degradation
function, we cannot recover the undegraded image [the inverse Fourier transform
of F(u, v) ] exactly because N(u, v) is not known. There is more bad news. If the deg-
radation function has zero or very small values, then the ratio N (u, v) H (u, v) could
easily dominate the term F(u, v). In fact, this is frequently the case, as you will see
shortly.
One approach to get around the zero or small-value problem is to limit the filter
frequencies to values near the origin. From the discussion of Eq. (4-92), we know
that H(0, 0) is usually the highest value of H(u, v) in the frequency domain. Thus,
by limiting the analysis to frequencies near the origin, we reduce the likelihood of
encountering zero values. The following example illustrates this approach.

EXAMPLE 5.9 : Image deblurring by inverse filtering.

The image in Fig. 5.25(b) was inverse filtered with Eq. (5-78) 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. The M 2 and N 2 constants are offset values; they center the function so that it will
correspond with the centered Fourier transform, as discussed in the previous chapter. (Remember, we
do not use padding with these functions.) In this case, M = N = 480. We know that a Gaussian function
has no zeros, so that will not be a concern here. However, despite this, the degradation values became so
small that the result of full inverse filtering [Fig. 5.27(a)] is useless. The reasons for this poor result are
as discussed in connection with Eq. (5-79).
Figures 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. The cut off was implemented by applying to the ratio a Butter-
worth lowpass function of order 10. This provided a sharp (but smooth) transition at the desired radius.
Radii near 70 yielded the best visual results [Fig. 5.27(c)]. Radii below 70 resulted in blurred images, as
in Fig. 5.27(b), which was obtained using a radius of 40. Values above 70 started to produce degraded
images, as illustrated in Fig. 5.27(d), which was obtained using a radius of 85. The image content is almost
visible in this image behind a “curtain” of noise, but the noise definitely dominates the result. Further
increases in radius values produced images that looked more and more like Fig. 5.27(a).

The results in the preceding example are illustrative of the poor performance of
direct inverse filtering in general. The basic theme of the three sections that follow is
how to improve on direct inverse filtering.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 357 6/16/2017 2:07:38 PM

 358 Chapter 5 Image Restoration and Reconstruction

a b
c d
FIGURE 5.27
Restoring
Fig. 5.25(b)
using Eq. (5-78).
(a) Result of using
the full filter.
(b) Result with H
cut off outside a
radius of 40.
(c) Result with H
cut off outside a
radius of 70.
(d) Result with H
cut off outside a
radius of 85.

5.8 MINIMUM MEAN SQUARE ERROR (WIENER) FILTERING

5.8

The inverse filtering approach discussed in the previous section makes no explicit
provision for handling noise. In this section, we discuss an approach that incorpo-
rates both the degradation function and statistical characteristics of noise into the
restoration process. The method is founded on considering images and noise as ran-
dom variables, and the objective is to find an estimate fˆ of the uncorrupted image f
such that the mean square error between them is minimized. This error measure is
defined as
⎧ ⁄ ⎫
e 2 = E ⎨( f − f )2 ⎬
⎩ ⎭ (5-80)

where E {i} is the expected value of the argument. We assume that the noise and the
image are uncorrelated, that one or the other has zero mean, and that the intensity
levels in the estimate are a linear function of the levels in the degraded image. Based

[Link]

DIP4E_GLOBAL_Print_Ready.indb 358 6/16/2017 2:07:38 PM

 5.8 Minimum Mean Square Error (Wiener) Filtering 359

on these assumptions, the minimum of the error function in Eq. (5-80) 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) + Sh (u, v) ⎥⎦
⎡ H * (u, v) ⎤
=⎢ 2
⎥ G(u, v) (5-81)
⎢⎣ H (u, v) + Sh (u, v) S f (u, v) ⎥⎦
⎡ 1 H (u, v)
2 ⎤
=⎢ ⎥ G(u, v)
⎢⎣ H (u, v) H (u, v) + Sh (u, v) S f (u, v) ⎥⎦
2

where we used the fact that the product of a complex quantity with its conjugate
is equal to the magnitude of the complex quantity squared. This result is known as
the Wiener filter, after N. Wiener [1942], who first proposed the concept in the year
shown. The filter, which consists of the terms inside the brackets, also is commonly
referred to as the minimum mean square error filter or the least square error filter.
We include references at the end of the chapter to sources containing detailed deri-
vations of the Wiener filter. Note from the first line in Eq. (5-81) that the Wiener
filter does not have the same problem as the inverse filter with zeros in the degrada-
tion function, unless the entire denominator is zero for the same value(s) of u and v.
The terms in Eq. (5-81) are as follows:

1. Fˆ (u, v) = Fourier transform of the estimate of the undegraded image.

2. G(u, v) = Fourier transform of the degraded image.
3. H(u, v) = degradation transfer function (Fourier transform of the spatial
degradation).
4. H ∗ (u, v) = complex conjugate of H(u, v) .
5. H (u, v) = H ∗ (u, v)H (u, v).
2

2
6. Sh (u, v) = N (u, v) = power spectrum of the noise [see Eq. (4-89)]†
2
7. S f (u, v) = F (u, v) = power spectrum of the undegraded image.

The restored image in the spatial domain is given by the inverse Fourier transform
of the frequency-domain estimate Fˆ (u, v). Note that if the noise is zero, then the
noise power spectrum vanishes and the Wiener filter reduces to the inverse filter.
Also, keep in mind the discussion at the end of Section 5.5 regarding the fact that all
transform work in this chapter is done without padding.

2
†
The term N(u, v) also is referred to as the autocorrelation of the noise. This term comes from the correlation
theorem (first line of entry 7 in Table 4.4). When the two functions are the same, correlation becomes autocorrela-
tion and the right side of that entry becomes H ∗ (u, v)H (u, v), which is equal to H(u, v) . Similar comments apply
2
2
to F(u, v) , which is the autocorrelation of the image. We will discuss correlation in more detail in Chapter 12.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 359 6/16/2017 2:07:40 PM

 360 Chapter 5 Image Restoration and Reconstruction

A number of useful measures are based on the power spectra of noise and of the
undegraded image. One of the most important is the signal-to-noise ratio, approxi-
mated using frequency domain quantities such as

M −1 N −1 M −1 N −1
∑∑ ∑∑
2 2
SNR = F (u, v) N (u, v) (5-82)
u=0 v=0 u=0 v=0

This ratio gives a measure of the level of information-bearing signal power (i.e., of
the original, undegraded image) to the level of noise power. An image with low
noise would tend to have a high SNR and, conversely, the same image with a higher
level of noise would have a lower SNR. This ratio is an important measure used in
characterizing the performance of restoration algorithms.
The mean square error given in statistical form in Eq. (5-80) can be approximated
also in terms of a summation involving the original and restored images:

1 M −1 N −1 ⁄ 2
MSE = ∑ ∑
MN x = 0 y = 0
[ f ( x, y) − f ( x, y)] (5-83)

In fact, 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 M −1 M −1 2
∑ ∑ fˆ ( x, y ) ∑ ∑ ⎡⎣ f ( x, y) − fˆ ( x, y)⎤⎦
2
SNR = (5-84)
x=0 y=0 x=0 y=0

The closer f and fˆ are, the larger this ratio will be. Sometimes the square root of the
preceding two measures is used instead, in which case they are referred to as the
root-mean-square-error and the root-mean-square-signal-to-noise ratio, respectively.
As we have mentioned before, keep in mind that quantitative measures do not nec-
essarily relate well to perceived image quality.
When dealing with white noise, the spectrum is a constant, which simplifies things
considerably. However, the power spectrum of the undegraded image seldom is
known. An approach frequently used when these quantities are not known, or can-
not be estimated, is to approximate Eq. (5-81) by the expression

⎡ 1 H (u, v)
2
⎤
Fˆ (u, v) = ⎢ 2
⎥ G(u, v) (5-85)
⎢⎣ H (u, v) H (u, v) + K ⎥⎦

2
where K is a specified constant that is added to all terms of H(u, v) . The following
examples illustrate the use of this expression.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 360 6/16/2017 2:07:40 PM

 5.8 Minimum Mean Square Error (Wiener) Filtering 361

a b c
FIGURE 5.28 Comparison of inverse and Wiener filtering. (a) Result of full inverse filtering of Fig. 5.25(b). (b) Radially
limited inverse filter result. (c) Wiener filter result.

EXAMPLE 5.10 : Comparison of deblurring by inverse and Wiener filtering.

Figure 5.28 illustrates the advantage of Wiener filtering over direct inverse filtering. Figure 5.28(a) is the
full inverse-filtered result from Fig. 5.27(a). Similarly, Fig. 5.28(b) is the radially limited inverse filter result
of Fig, 5.27(c). These images are duplicated here for convenience in making comparisons. Figure 5.28(c)
shows the result obtained using Eq. (5-85) with the degradation function used in Example 5.9. The value
of K was chosen interactively to yield the best visual results. The advantage of Wiener filtering over the
direct inverse approach is evident in this example. By comparing Figs. 5.25(a) and 5.28(c), we see that
the Wiener filter yielded a result very close in appearance to the original, undegraded image.

EXAMPLE 5.11 : More deblurring examples using Wiener filtering.

The first row of Fig. 5.29 shows, from left to right, the blurred image of Fig. 5.26(b) heavily corrupted by
additive Gaussian noise of zero mean and variance of 650; the result of direct inverse filtering; and the
result of Wiener filtering. The Wiener filter of Eq. (5-85) was used, with H(u, v) from Example 5.8, and
with K chosen interactively to give the best possible visual result. As expected, direct inverse filtering
produced an unusable image. Note that the noise in the inverse filtered image is so strong that it masks
completely the content of the image. The Wiener filter result is by no means perfect, but it does give us
a hint as to image content. The text can be read with moderate effort.
The second row of Fig. 5.29 shows the same sequence just discussed, but with the level of the noise
variance reduced by one order of magnitude. This reduction had little effect on the inverse filter, but
the Wiener results are considerably improved. For example, the text is much easier to read now. In the
third row of Fig. 5.29, the noise variance was reduced more than five orders of magnitude from the first
row. In fact, image in Fig. 5.29(g) has no visible noise. The inverse filter result is interesting in this case.
The noise is still quite visible, but the text can be seen through a “curtain” of noise (see Problem 5.30).
The Wiener filter result in Fig. 5.29(i) is excellent, being quite close visually to the original image in Fig.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 361 6/16/2017 2:07:41 PM

 362 Chapter 5 Image Restoration and Reconstruction

a b c
d e f
g h i
FIGURE 5.29 (a) 8-bit image corrupted by motion blur and additive noise. (b) Result of inverse filtering. (c) Result of
Wiener filtering. (d)–(f) Same sequence, but with noise variance one order of magnitude less. (g)–(i) Same sequence,
but noise variance reduced by five orders of magnitude from (a). Note in (h) how the deblurred image is quite vis-
ible through a “curtain” of noise.

[Link]

DIP4E_GLOBAL_Print_Ready.indb 362 6/16/2017 2:07:41 PM

 5.9 Constrained Least Squares Filtering 363

5.26(a). In practice, the results of restoration filtering are seldom this close to the original images. This
example, and Example 5.12 in the next section, were idealized slightly to focus on the effects of noise
on restoration algorithms.

5.9 CONSTRAINED LEAST SQUARES FILTERING

5.9

The problem of having to know something about the degradation function H is com-
mon to all methods discussed in this chapter. However, the Wiener filter presents
an additional difficulty: the power spectra of the undegraded image and noise must
be known also. We showed in the previous section that in some cases it is possible
to achieve acceptable results using the approximation in Eq. (5-85), but a constant
value for the ratio of the power spectra is not always a suitable solution.
The method discussed in this section requires knowledge of only the mean and
variance of the noise. As discussed in Section 5.2, these parameters generally can be
calculated from a given degraded image, so this is an important advantage. Another
difference is that the Wiener filter is based on minimizing a statistical criterion and,
as such, it is optimal in an average sense. The algorithm presented in this section
has the notable feature that it yields an optimal result for each image to which it
is applied. Of course, it is important to keep in mind that these optimality criteria,
although they are comforting from a theoretical point of view, are not related to
the dynamics of visual perception. As a result, the choice of one algorithm over the
other will almost always be determined by the perceived visual quality of the result-
ing images.
By using the definition of convolution given in Eq. (4-94), and as explained in
Section 2.6, we can express Eq. (5-64) in vector-matrix form:
g = Hf + H (5-86)

For example, suppose that g( x, y) is of size M × N . We can form the first N elements
of vector g by using the image elements in the first row of g( x, y), the next N ele-
ments from the second row, and so on. The dimensionality of the resulting vector will
be MN × 1. These are also the dimensions of f and H, as these vectors are formed in
the same manner. Matrix H then has dimensions MN × MN . Its elements are given
by the elements of the convolution in Eq. (4-94).
It would be reasonable to arrive at the conclusion that the restoration problem
can now be reduced to simple matrix manipulations. Unfortunately, this is not the
case. For instance, suppose that we are working with images of medium size, say
M = N = 512. Then the vectors in Eq. (5-86) would be of dimension 262, 144 × 1
and matrix H would be of dimension 262, 144 × 262, 144. Manipulating vectors and
matrices of such sizes is not a trivial task. The problem is complicated further by
the fact that H is highly sensitive to noise (after the experiences we had with the
effect of noise in the previous two sections, this should not be a surprise). The key
See Gonzalez and Woods
[1992] for an entire chap- advantage of formulating the restoration problem in matrix form is that it facilitates
ter devoted to the topic derivation of restoration algorithms.
of algebraic techniques
for image restoration. Although we do not fully derive the method of constrained least squares that
we are about to present, this method has its roots in a matrix formulation. We give

[Link]

DIP4E_GLOBAL_Print_Ready.indb 363 6/16/2017 2:07:42 PM