Introduction

▪ Principal goal of restoration techniques is to improve

an image in some predefined sense.
▪ Image enhancement is largely a subjective process,
while image restoration is an objective process.
▪ Restoration attempts to recover an image that has
been degraded by using a priori knowledge of the
degradation phenomenon.
▪ Restoration techniques are oriented toward modeling
the degradation and applying the inverse process in
order to recover the original image.
▪ Approach usually involves formulating a criterion of
goodness that will yield an optimal estimate of the
desired result.
▪ Some restoration techniques are best formulated in
the spatial domain, while others are better suited for
the frequency domain.
▪ Spatial processing is applicable when the only
degradation is additive noise.
▪ Degradations such as image blur are difficult to
approach in the spatial domain using small filter
masks. So, we use frequency domain models.
Model of the Image Degradation/
Restoration Process
▪ Degradation process is modeled as a degradation function
that, together with an additive noise term, operates on an
input image f(x, y) to produce a degraded image g(x, y).
▪ Given g(x, y), some knowledge about the degradation
function 𝐻, and some knowledge about the additive noise
term 𝜂(𝑥, 𝑦), the objective of restoration is to obtain an
estimate 𝑓መ (𝑥, 𝑦) of the original image.
▪ The more we know about 𝐻 and 𝜂, the closer
𝑓መ (𝑥, 𝑦) will be to 𝑓(𝑥, 𝑦).
▪ If H is linear, position-invariant process, degraded
image in spatial domain is:
𝑔(𝑥, 𝑦) = ℎ(𝑥, 𝑦) ∗ 𝑓(𝑥, 𝑦) + 𝜂(𝑥, 𝑦)
▪ Where ℎ(𝑥, 𝑦) is the spatial representation of the
degradation function.
▪ In frequency domain, the representation is:
𝐺(𝑢, 𝑣) = 𝐻(𝑢, 𝑣)𝐹(𝑢, 𝑣) + 𝑁(𝑢, 𝑣)
▪ In future, we will assume H is the identity operator,
and we deal only with degradations due to noise.
Noise Models
Introduction
▪ Principal sources of noise in digital images arise
during image acquisition and/or transmission.
▪ Performance of imaging sensors is affected by a
variety of factors, such as environmental conditions
during image acquisition, and by the quality of the
sensing elements themselves.
▪ Images corrupted during transmission principally due
to interference in the channel used for transmission.
Spatial & Frequency Properties of Noise
▪ Spatial – Is noise correlated with the image in spatial
domain?
▪ Frequency properties refer to the frequency content
of noise in the Fourier sense.
▪ With the exception of spatially periodic noise we
assume that noise is independent of spatial
coordinates, and that it is uncorrelated with respect to
the image itself.
Some Important Noise Probability Density Functions
are:
➢ Gaussian noise

➢ Rayleigh noise

➢ Erlang (gamma) noise

➢ Exponential noise

➢ Uniform noise

➢ Impulse (Salt-and-pepper) noise

Gaussian (Normal) Noise
▪ Used frequently in practice due to mathematical
tractability in both spatial and frequency domains.
▪ Even used whenever such models are even marginally
applicable at best.
▪ PDF of a Gaussian random variable, z, is given by:
1 − 𝑧−𝑧ҧ 2 /2𝜎 2
𝑝 𝑧 = 𝑒
2𝜋𝜎 2
Where 𝑧 represents the intensity, 𝑧ҧ represents the mean
value of 𝑧, and 𝜎 is the standard deviation. 𝜎 2 is the
variance of 𝑧.
▪ Usually arises due to factors like electronic circuit noise
and sensor noise due to poor illumination and/or high
temperature.
▪ When 𝑧 is described by the equation, approximately
70% of its values will be in the range [(𝑧-𝜎),
ҧ (𝑧+𝜎)],
ҧ
and about 95% will be in the range [(𝑧-2𝜎),
ҧ (𝑧ҧ +2𝜎)]
Rayleigh noise
▪ The PDF of Rayleigh noise is given
by:
2 2 /𝑏
(𝑧 − 𝑎)𝑒 − 𝑧−𝑎 𝑓𝑜𝑟 𝑧 ≥ 𝑎
𝑝 𝑧 = 𝑏
𝑓𝑜𝑟 𝑧 < 𝑎
0
▪ Mean and variance of this density are:
𝑏 4−𝜋
𝑧ҧ = a + 𝜋𝑏/4 and 𝜎2 =
4
▪ Note the displacement in origin, and
that basic shape of the density is
skewed to the right.
▪ This can be useful for approximating
skewed histograms.
▪ Helpful in characterizing noise
phenomena in range imaging.
Erlang (Gamma) Noise
▪ PDF of Erlang noise is given by:
𝑎𝑏 𝑧 𝑏−1 −𝑎𝑧 𝑓𝑜𝑟 𝑧 ≥ 0
𝑝 𝑧 = 𝑏 −1 !𝑒
𝑓𝑜𝑟 𝑧 < 0
0
▪ Where parameters are such that
a>0, b is a positive integer.
▪ Mean and variance are given by:
𝑏 2 𝑏
𝑧ҧ = and 𝜎 =
𝑎 𝑎2
▪ This is also referred to as gamma
density or Erlang Density
▪ Common application in laser
imaging.
Exponential Noise
▪ PDF of Exponential noise is
given by:
𝑎𝑒 −𝑎𝑧 𝑓𝑜𝑟 𝑧 ≥ 0
𝑝 𝑧 =
0 𝑓𝑜𝑟 𝑧 < 0
▪ Where 𝑎 > 0.
▪ Mean and variance are given
by:
1 2 1
𝑧ҧ = and 𝜎 =
𝑎 𝑎2
▪ This is a special case of the
Erlang PDF, with 𝑏 = 1.
▪ Common application in laser
imaging.
Uniform Noise
▪ PDF of Uniform noise is given
by:
1
𝑓𝑜𝑟 𝑎 ≤ 𝑧 ≤ 𝑏
𝑝 𝑧 = 𝑏−𝑎
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
0
▪ Where 𝑎 > 0.
▪ Mean and variance are given
by:
𝑎+𝑏 (𝑏−𝑎)2
𝑧ҧ = and 𝜎2 =
2 12
▪ Least descriptive in practical
situations. Used as basis for
random number generators for
simulation.
Impulse (salt-and-pepper) noise
▪ The PDF of a (bipolar) impulse noise is
given by:
𝑃𝑎 𝑓𝑜𝑟 𝑧 = 𝑎
▪ 𝑝 𝑧 = 𝑃𝑏 𝑓𝑜𝑟 𝑧 = 𝑏
0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
▪ If 𝑏 > 𝑎, intensity b will appear as light
dot in the image.
▪ Conversely, level 𝑎 will appear as a dark
dot.
▪ If either 𝑃𝑎 or 𝑃𝑏 is zero, impulse noise
is called unipolar.
▪ If neither probability is zero, and
approximately equal, impulse noise
values will resemble salt and pepper
granules randomly distributed over
image. Hence the name.
▪ Aka Data-drop-out and spike noise.
▪ Noise impulses can be positive or negative.
▪ Scaling is part of image digitizing process.
▪ Since impulse corruption is usually large compared
to image strength, impulse noise is generally
digitized as extreme values (black or white).
▪ Thus, we assume a and b are saturated values (min or
max allowed values in digitized image)
▪ Positive impulses appear as white (salt) noise.
▪ Impulse noise is found in situations where quick
transients, such as faulty switching, take place during
imaging.
Periodic Noise
▪ Periodic noise in an image arises typically from
electrical or electromechanical interference during
image acquisition.
▪ Only type of spatially dependent noise considered.
▪ Can be reduced significantly via frequency domain
filtering.
▪ Fourier Transform of pure sinusoid is a pair of
conjugate impulses located at the conjugate
frequencies of sine wave.
▪ If amplitude of a sine wave in spatial domain is
strong enough, we see in the spectrum, a pair of
impulses for each sine wave.
Estimation of Noise Parameters
▪ The parameters of periodic noise typically are
estimated by inspection of the Fourier spectrum of
the image.
▪ Periodic noise produces frequency spikes that can
often be detected by visual analysis.
▪ Another approach is to attempt to infer periodicity of
noise components directly from image, but not
always possible.
▪ Automated analysis possible when noise spikes are
either exceptionally pronounced or when knowledge
is available about general location of frequency
components of interference.
▪ Parameters of noise PDF may be partially known
from sensor specifications, but often need to estimate
for imaging arrangement.
▪ How do we do this?
▪ In case of imaging system being available, one way
to study noise characteristics is to capture a set of
images of a “flat” environment. Eg. Imaging a solid
gray board that is uniformly illuminated.
▪ Resulting images are good indicators of system
noise.
▪ In case only images are available, PDF parameters
can be estimated from small patches of reasonably
constant background intensity.
• In case only images are available, PDF parameters can be
estimated from small patches of reasonably constant background
intensity.
• Shape of the histogram identifies closest PDF match.
• Parameters can be estimated from the graphs.
• In case of impulse noise, we need actual probability of
occurrence of white and black pixels. So, we choose a mid-gray,
relatively constant area for analysis.
Restoration in the
Presence of Noise
Only
Spatial Filtering
Introduction
▪ When noise is the only degradation in an image:
𝑔(𝑥, 𝑦) = 𝑓(𝑥, 𝑦) + 𝜂(𝑥, 𝑦)
𝐺(𝑢, 𝑣) = 𝐹(𝑢, 𝑣) + 𝑁(𝑢, 𝑣)
▪ Noise terms are unknown. Subtracting them from
g(x, y) or 𝐺(𝑢, 𝑣)is not a realistic option.
▪ In the case of periodic noise, it usually is possible to
estimate N(u, v) from the spectrum of 𝐺(𝑢, 𝑣).
▪ In this case N(u, v) can be subtracted from G(u, v) to
obtain an estimate of the original image. However, it
is very rare.
▪ Spatial filtering is the method of choice in situations
when only additive random noise is present.
Mean Filters
Arithmetic Mean Filters
▪ Simplest of the mean filters.
▪ Let 𝑆𝑥𝑦 represent the set of coefficients in rectangular
subimage window (neighborhood) of size m x n, centered at
point (x, y).
▪ Arithmetic mean filter computes the average value of the
corrupted image g(x, y) in the area defined by 𝑆𝑥𝑦
▪ The value of the restored image 𝑓መ at point (x, y) is simply the
arithmetic mean computed using the pixels in the region
defined by 𝑆𝑥𝑦
1
𝑓መ 𝑥, 𝑦 = ෍ 𝑔(𝑠, 𝑡)
𝑚𝑛
(𝑠,𝑡)∈𝑆𝑥𝑦
▪ This operation can be implemented using a spatial filter of
size m x 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
𝑚𝑛
𝑓መ 𝑥, 𝑦 = ෑ 𝑔(𝑠, 𝑡)
(𝑠,𝑡)∈𝑆𝑥𝑦

▪ Each restored pixel is given by the product of the

pixels in the subimage window, raised to the power
1/mn.
▪ A geometric mean filter achieves smoothing
comparable to the arithmetic mean filter, but it tends
to lose less image detail in the process.
Harmonic mean filter
▪ The harmonic mean filtering operation is given by
the expression:

𝑚𝑛
𝑓መ 𝑥, 𝑦 =
1
σ(𝑠,𝑡)∈𝑆𝑥𝑦
𝑔(𝑠, 𝑡)

▪ 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:
σ(𝑠,𝑡)∈𝑆𝑥𝑦 𝑔(𝑠, 𝑡)𝑄+1
𝑓መ 𝑥, 𝑦 =
σ(𝑠,𝑡)∈𝑆𝑥𝑦 𝑔(𝑠, 𝑡)𝑄
▪ where Q is called the order of the filter.
▪ This filter is well suited for reducing or virtually
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.
▪ In general, the arithmetic and geometric mean filters
(particularly the latter) are well suited for random
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.
Order-Statistic Filters
▪ Order-statistics filters are nonlinear spatial filters
whose response is based on ordering (ranking) the
pixels contained in the image area encompassed by
the filter.
▪ Value of the center pixel is replaced with the value
determined by the ranking result.
▪ The ranking result determines the response of the
filter.
Median filter
▪ Replaces the value of a pixel by the median of the
intensity levels in the neighborhood of that pixel:

𝑓መ 𝑥, 𝑦 = median{𝑔(𝑠, 𝑡)}
(𝑠,𝑡)∈𝑆𝑥𝑦

▪ The value of the pixel at (x, y) is included in the

computation of the median.
▪ Median filters are quite popular . For certain types of
random noise, they provide excellent noise reduction
capabilities, with considerably less blurring than
linear smoothing fillers of similar size.
▪ Median filters are particularly effective in the
presence of both bipolar and unipolar impulse noise.
▪ The median, ξ, of a set of values is such that half the
values in the set are less than or equal to ξ, and half are
greater than or equal to ξ.
▪ In order to perform median filtering at a point in an
image, we first sort the values of the pixel in question
and its neighbors, determine their median, and assign this
value to that pixel.
▪ When several values in a neighborhood are the same, all
equal values are grouped.
▪ The principal function of median filters is to force points
with distinct gray levels to be more like their neighbors.
Max and Min filters
▪ The median represents the 50th percentile of a
ranked set of numbers. but ranking lends itself to
many other possibilities.
▪ Using the 100th percentile results in the so-called
max filter, given by

𝑓መ 𝑥, 𝑦 = max {𝑔(𝑠, 𝑡)}

(𝑠,𝑡)∈𝑆𝑥𝑦

▪ This is useful for finding brightest points in an

image.
▪ Since pepper noise has very low values, it is reduced
by max filter due to max selection process.
▪ Using the 0th percentile results in the so-called min
filler, given by

𝑓መ 𝑥, 𝑦 = min {𝑔(𝑠, 𝑡)}

(𝑠,𝑡)∈𝑆𝑥𝑦

▪ This is useful for finding darkest points in an image.

▪ Since salt noise has very high values, it is reduced by
min filter due to min selection process.
Midpoint Filter
▪ The midpoint filter simply computes the midpoint
between the maximum and minimum values in the
area encompassed by the filter:
1
𝑓መ 𝑥, 𝑦 = max {𝑔(𝑠, 𝑡)} + min {𝑔(𝑠, 𝑡)}
2 (𝑠,𝑡)∈𝑆𝑥𝑦 (𝑠,𝑡)∈𝑆𝑥𝑦

▪ This filter combines order statistics and averaging. It

works best for randomly 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(s, t) in the neighborhood Sxy.
▪ Let gr(s, t) represent the remaining 𝑚𝑛– 𝑑 pixels.
▪ A filter formed by averaging these remaining pixels is
called an alpha trimmed mean filter:
1
𝑓መ 𝑥, 𝑦 = ෍ 𝑔𝑟 (𝑠, 𝑡)
𝑚𝑛 − 𝑑
(𝑠,𝑡)∈𝑆𝑥𝑦
▪ where the value of d can range from 0 to 𝑚𝑛 − 1.
▪ When 𝑑 = 0, the alpha trimmed filter reduces to the
arithmetic mean filter .
▪ If we choose 𝑑 = 𝑚𝑛 − 1. the filter becomes a median
filter.
▪ For other values of 𝑑, the alpha-trimmed filler is useful
in situations involving multiple types of noise, such as a
combination of salt-and-pepper and Gaussian noise.
Adaptive Filters
Introduction
▪ Earlier filters are applied to image without regard for
variance of image characteristics from point to point.
▪ Adaptive Filters - behavior changes based on
statistical characteristics of the image inside the filter
region defined by the m x n rectangular window 𝑆𝑥𝑦
▪ Adaptive Filters are capable of superior performance.
▪ Price paid for improved filtering power – increased
complexity.
▪ These filters will work for additive noise situations
only.
Adaptive, Local Noise Reduction Filter
▪ Simplest statistical measures of a random variable : mean
and variance.
▪ Mean gives a measure of average intensity in the region
over which the mean is computed, Variance gives a
measure of contrast in that region.
▪ Our filter is to operate on a local region. 𝑆𝑥𝑦 .
▪ The response of the filter at any point (x, y) on which the
region is centered is to be based on four quantities:
o 𝒈(𝒙, 𝒚) – Value of noisy image at (𝑥, 𝑦)
o 𝝈𝟐𝒏 - Variance of the noise corrupting 𝑓(𝑥, 𝑦) to make 𝑔(𝑥, 𝑦)
o 𝒎𝑳 - Local mean of pixels in 𝑆𝑥𝑦
o 𝝈𝟐𝑳 - Local variance of the pixels in 𝑆𝑥𝑦
▪ We want filter behaviour to be as follows:
1. If 𝝈𝟐𝒏 is zero, the filter should return simply the
value of 𝑔(𝑥, 𝑦).This is the trivial, zero-noise case
in which 𝑔(𝑥, 𝑦) is equal to 𝑓(𝑥, 𝑦).
2. If the local variance is high relative to 𝝈𝟐𝒏 , the filter
should return a value close to 𝑔(𝑥, 𝑦). A high local
variance typically is associated with edges, and
these should he preserved.
3. If the two variances are equal, the filter should
return the arithmetic mean value of the pixels in
𝑆𝑥𝑦 . This condition occurs when the local area has
the same properties as the overall image, and local
noise is to be reduced simply by averaging.
▪ An adaptive expression for obtaining 𝑓 𝑥, 𝑦 based
on these assumptions may be written as:
𝜎𝑛2
𝑓መ 𝑥, 𝑦 = 𝑔 𝑥, 𝑦 − 2 𝑔 𝑥, 𝑦 − 𝑚𝐿
𝜎𝐿
▪ The only quantity that needs to be known/estimated
is the variance of the overall noise, 𝜎𝑛2 .
▪ The other parameters are computed from the pixels
in 𝑆𝑥𝑦 at each location 𝑥, 𝑦 , on which the filter
window is centered.
▪ We assume 𝜎𝑛2 ≤ 𝜎𝐿2 , since Noise is additive and
position invariant.
▪ We seldom have exact knowledge of 𝜎𝑛2 , so the
condition can be violated in practice.
▪ So, a test should be built into an implementation of
the equation where ratio is set to 1, if 𝜎𝑛2 > 𝜎𝐿2 .
▪ This makes filter nonlinear, but prevents nonsensical
results due to lack of knowledge about variance of
image noise.
▪ Another approach - allow the negative values to
occur, and then rescale the intensity values at the end
– Loss of Image Dynamic Range.
Adaptive median filter
▪ Earlier median filter works well if spatial density of
impulse noise is low (Pa and Pb < 0.2)
▪ Adaptive filters can handle impulse noise with larger
probabilities.
▪ Additional benefit – seeks to preserve detail while
smoothing non-impulse noise.
▪ Adaptive median filter also works in a rectangular
window area 𝑆𝑥𝑦 . However, adaptive median filter
changes (increases) the size of 𝑆𝑥𝑦 during filter
operation, depending on certain conditions.
▪ Output of the filter is a single value used to replace
the value used to replace value of the pixel at (x, y),
the point on which the window 𝑆𝑥𝑦 is centered at a
given time.
▪ 3 major purposes of algorithm:
o Remove salt and pepper (impulse) noise
o Provide smoothing other non-impulse noise
o Reduce distortion like extensive thickening/thinning of object
boundaries.

▪ Consider the following notation:

o Zmin = minimum intensity value in 𝑆𝑥𝑦
o Zmax = maximum intensity value in 𝑆𝑥𝑦
o Zmed = median of intensity values in 𝑆𝑥𝑦
o 𝑍𝑥𝑦 = intensity value at coordinates (x, y)
o Smax = maximum allowed size of 𝑆𝑥𝑦

▪ Adaptive median filtering algorithm works in two stages:

Stage A and stage B as shown in the next slide:
Stage A: 𝐴1 = 𝑍𝑚𝑒𝑑 − 𝑍𝑚𝑖𝑛
𝐴2 = 𝑍𝑚𝑒𝑑 − 𝑍𝑚𝑎𝑥
If A1> 0 AND 𝐴2 < 0, go to stage B
Else increase the window size
If 𝑤𝑖𝑛𝑑𝑜𝑤 𝑠𝑖𝑧𝑒 ≤ 𝑆𝑚𝑎𝑥 repeat stage A
Else output 𝑍𝑚𝑒𝑑

Stage B: 𝐵1 = 𝑍𝑥𝑦 − 𝑍𝑚𝑖𝑛

𝐵2 = 𝑍𝑥𝑦 − 𝑍𝑚𝑎𝑥
If 𝐵1 > 0 AND 𝐵2 < 0, output 𝑍𝑥𝑦
Else output 𝑍𝑚𝑒𝑑
▪ The values 𝑍𝑚𝑖𝑛 and 𝑍𝑚𝑎𝑥 are considered to be
"impulse-like" noise components, even if these are
not the lowest and highest possible pixel values in
the image.
▪ Purpose of Stage A is to determine if median filter
output, 𝑍𝑚𝑒𝑑 , is an impulse (black or white) or not.
▪ If 𝑍𝑚𝑖𝑛 < 𝑍𝑚𝑒𝑑 < 𝑍𝑚𝑎𝑥 holds, then 𝑍𝑚𝑒𝑑 cannot be
impulse. So, we go to Stage B and test to see if the
point in center of window, 𝑍𝑥𝑦 is itself an impulse.
▪ If 𝑍𝑚𝑖𝑛 < 𝑍𝑥𝑦 < 𝑍𝑚𝑎𝑥 , then 𝑍𝑥𝑦 cannot be an
impulse, so algorithm outputs unchanged pixel value
𝑍𝑥𝑦 .
• In (b) Although the noise was effectively removed, the fitter caused significant 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.

• In (c) Noise removal performance was similar to the median filter. However, the
adaptive filter did a better job of preserving 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
Periodic Noise
Reduction by
Frequency Domain
Filtering
▪ Periodic noise can be analyzed and filtered quite
effectively using frequency domain techniques.
▪ Periodic noise appears as concentrated bursts of
energy in the Fourier transform, at locations
corresponding to the frequencies of the periodic
interference.
▪ Use a selective filter to isolate the noise.
▪ 3 Types of Selective Filters:
o Bandreject filter
o Bandpass filter
o Notch filter
Bandreject Filters
▪ Principal application of bandreject filtering is for
noise removal in applications where the general
location of the noise component(s) in the frequency
domain is approximately known.
o Eg: image corrupted by additive periodic noise that can be
approximated as two-dimensional sinusoidal functions.
▪ Fourier transform of a sine consists of two impulses
that are mirror images of each other about the origin
of the transform. Their locations are:

▪ The impulses are both imaginary (the real part of the

Fourier transform of a sine is zero) and are complex
conjugates of each other.
Bandpass Filters
▪ Bandpass filter performs the opposite operation of a
bandreject filter.
▪ Bandpass Filter is obtained from bandreject filter as:
𝐻𝐵𝑃 𝑢, 𝑣 = 1 − 𝐻𝐵𝑅 𝑢, 𝑣
▪ Performing straight bandpass filtering on an image is
not a common procedure because it generally
removes too much image detail.
▪ Bandpass filtering is useful in isolating the effects on
an image caused by selected frequency bands.
▪ This is a useful result because it simplifies analysis
of the noise, reasonably independently of image
content.
Notch Filters
▪ Rejects (or passes) frequencies in predefined
neighborhoods about a center frequency.
▪ Due to the symmetry of the Fourier transform, notch
filters must appear in symmetric pairs about the
origin in order to obtain meaningful results.
▪ One exception - if the notch filter is located at the
origin, in which case it appears by itself.
▪ Number of pairs of notches and shapes are arbitrary.
▪ We can also obtain notch fillers that pass, rather than
suppress, the frequencies contained in the notch
areas.
▪ Their transfer functions are given by:
𝐻𝑁𝑃 𝑢, 𝑣 = 1 − 𝐻𝑁𝑅 𝑢, 𝑣
 Optimum Notch Filtering
▪ Figure below is another example of periodic image
degradation, shows a digital image of the Martian terrain
taken by the Mariner 6 spacecraft.
▪ The interference pattern is somewhat similar to the earlier
one (periodic noise), but this pattern is more subtle and
harder to detect in the frequency plane.

▪ The starlike components were caused by the interference.

and several pairs of components are present, indicating that
the pattern contains more than just one sinusoidal
component.
▪ In such cases, earlier methods cannot always be
used:
o May remove too much image information during filtering
o Interference components generally are not single-
frequency bursts, but have broad skirts that carry
information about interference patterns.
o These skirts are not always easily detectable from the
Fourier transform background.
▪ Alternative methods are useful in many applications.
▪ Optimum method minimizes local variances of
restored estimate 𝑓መ 𝑥, 𝑦 .
▪ The procedure consists of first isolating the principal
contributions of the interference pattern and then
subtracting a variable, weighted portion of the
pattern from the corrupted image.
▪ Step 1:
o Extract the principal frequency components of the
interference pattern.
o This can be done by placing a notch pass filter 𝐻𝑁𝑃 𝑢, 𝑣 ,
at the location of each spike.
o If the filter is constructed to pass only components
associated with the interference pattern, then the Fourier
transform of the interference noise pattern is given by the
expression:
𝑁 𝑢, 𝑣 = 𝐻𝑁𝑃 𝑢, 𝑣 𝐺(𝑢, 𝑣)
o After a particular filter has been selected, in spatial
domain:
𝜂 𝑥, 𝑦 = 𝔍−1 𝐻𝑁𝑃 𝑢, 𝑣 𝐺(𝑢, 𝑣)
▪ Step 2:
o To get the uncorrupted image 𝑓መ 𝑥, 𝑦 , we then subtract the
noise pattern from 𝑔(𝑥, 𝑦).
o Problem is, filtering procedure yields an approximation of
True noise pattern.
o Effect of components not present in 𝜂 𝑥, 𝑦 can be
minimized by subtracting a weighted portion of 𝜂 𝑥, 𝑦
from 𝑔 𝑥, 𝑦 :
𝑓መ 𝑥, 𝑦 = 𝑔 𝑥, 𝑦 − 𝑤 𝑥, 𝑦 𝜂 𝑥, 𝑦 ----(1)
o 𝑤 𝑥, 𝑦 is called weighting or modulation function, and is
chosen so as to optimize result.

▪ One approach is to select a 𝑤 𝑥, 𝑦 so that the

variance of 𝑓መ 𝑥, 𝑦 is minimized over a specified
neighborhood of every point 𝑥, 𝑦 .
▪ Consider a neighborhood of size (2𝑎 + 1) by (2𝑏 + 1)
about a point 𝑥, 𝑦
▪ The "local" variance of 𝑓መ 𝑥, 𝑦 at coordinates 𝑥, 𝑦 can be
estimated from the samples, as follows:
𝑎 𝑏
1 2
𝜎2 𝑥, 𝑦 = ෍ ෍ 𝑓መ 𝑥 + 𝑠, 𝑦 + 𝑡 − 𝑓መҧ 𝑥, 𝑦
(2𝑎 + 1)(2𝑏 + 1)
𝑠=−𝑎 𝑡=−𝑏

----(2)
▪ Where 𝑓መҧ 𝑥, 𝑦 - average value of 𝑓መ in the neighborhood:
𝑎 𝑏
1
𝑓መҧ 𝑥, 𝑦 = ෍ ෍ 𝑓መ 𝑥 + 𝑠, 𝑦 + 𝑡
(2𝑎 + 1)(2𝑏 + 1)
𝑠=−𝑎 𝑡=−𝑏

----(3)
▪ Points on or near the edge of the image can be treated by
considering partial neighborhoods or by padding the border
with 0s.
▪ Eqn (1) in (2) gives:
𝜎 2 𝑥, 𝑦
𝑎 𝑏
1
= ෍ ෍ ቂ[𝑔 𝑥 + 𝑠, 𝑦 + 𝑡
2𝑎 + 1 2𝑏 + 1
𝑠=−𝑎 𝑡=−𝑏
2
− 𝑤 𝑥 + 𝑠, 𝑦 + 𝑡 𝜂 𝑥 + 𝑠, 𝑦 + 𝑡 ሿ − 𝑔ҧ 𝑥, 𝑦 − 𝑤 𝑥, 𝑦 𝜂 𝑥, 𝑦 ቃ

▪ Assuming 𝑤 𝑥, 𝑦 remains constant over a

neighborhood, we can approximate as:
𝑤 𝑥 + 𝑠, 𝑦 + 𝑡 =𝑤 𝑥, 𝑦
For −𝑎 ≤ 𝑠 ≤ 𝑎 and −𝑏 ≤ 𝑡 ≤ 𝑏.
Also, we get:
𝑤 𝑥, 𝑦 𝜂 𝑥, 𝑦 = 𝑤 𝑥, 𝑦 ഥ𝜂 𝑥, 𝑦
In the neighborhood.
▪ With these approximations, we get:
𝜎 2 𝑥, 𝑦
𝑎 𝑏
1
= ෍ ෍ [[𝑔 𝑥 + 𝑠, 𝑦 + 𝑡 − 𝑤 𝑥, 𝑦 𝜂 𝑥 + 𝑠, 𝑦 + 𝑡 ሿ − [𝑔ҧ 𝑥, 𝑦
(2𝑎 + 1)(2𝑏 + 1)
𝑠=−𝑎 𝑡=−𝑏
− 𝑤 𝑥, 𝑦 ഥ𝜂 𝑥, 𝑦 ሿ ሿ2

▪ To minimize 𝜎 2 𝑥, 𝑦 , we solve:
𝜕𝜎 2 𝑥, 𝑦
=0
𝜕𝑤 𝑥, 𝑦
▪ For 𝑤 𝑥, 𝑦 . The result is:
𝑔 𝑥, 𝑦 𝜂 𝑥, 𝑦 − ഥ𝑔 𝑥, 𝑦 ഥ𝜂 𝑥, 𝑦
𝑤 𝑥, 𝑦 =
𝜂 2 𝑥, 𝑦 − ഥ𝜂 2 (𝑥, 𝑦)
▪ This can be used to obtain restored image 𝑓መ 𝑥, 𝑦 .
▪ Since 𝑤 𝑥, 𝑦 is assumed to be constant in a neighborhood,
we do not need to compute it each time in a non-overlapping
neighborhood.
Linear Position-
Invariant
Degradations
Introduction
▪ Degraded image in spatial domain is:
𝑔 𝑥, 𝑦 = 𝐻 𝑓 𝑥, 𝑦 + 𝜂(𝑥, 𝑦)
▪ Assuming 𝜂 𝑥, 𝑦 = 0 (no additive noise):
𝑔 𝑥, 𝑦 = 𝐻 𝑓 𝑥, 𝑦
▪ 𝐻 is said to be linear if:
𝐻 𝑎𝑓1 𝑥, 𝑦 + 𝑏𝑓2 𝑥, 𝑦 = 𝑎𝐻 𝑓1 𝑥, 𝑦 + 𝑏𝐻 𝑓2 𝑥, 𝑦
▪ Where 𝑎 and 𝑏 are scalars, 𝑓1 𝑥, 𝑦 and 𝑓2 𝑥, 𝑦 are
any two input images.
▪ When 𝑎 = 𝑏 = 1,
𝐻 𝑓1 𝑥, 𝑦 + 𝑓2 𝑥, 𝑦 = 𝐻 𝑓1 𝑥, 𝑦 + 𝐻 𝑓2 𝑥, 𝑦
▪ This is the additivity property
▪ When 𝑓2 𝑥, 𝑦 = 0, equation becomes:
𝐻 𝑎𝑓1 𝑥, 𝑦 = 𝑎𝐻 𝑓1 𝑥, 𝑦
▪ This is the property of homogeneity.
▪ Linear operator possesses both properties of
additivity and homogeneity.
▪ An operator having input-output relationship
𝑔 𝑥, 𝑦 = 𝐻 𝑓 𝑥, 𝑦 is said to be position (or space)
invariant if:
𝐻 𝑓 𝑥 − 𝛼, 𝑦 − 𝛽 = 𝑔 𝑥 − 𝛼, 𝑦 − 𝛽
▪ Response at any point in the image depends only on
the value of the input at that point, not its position.
▪ Consider impulse sifting property of a 2-D function:
∞ ∞
𝑓 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 𝛿(𝑥 − 𝛼, 𝑦 − 𝛽)𝑑𝛼 𝑑𝛽
−∞ −∞

▪ Considering the transformation:

∞ ∞
𝑔 𝑥, 𝑦 = 𝐻 𝑓 𝑥, 𝑦 =𝐻 න න 𝑓 𝛼, 𝛽 𝛿(𝑥 − 𝛼, 𝑦 − 𝛽)𝑑𝛼 𝑑𝛽
−∞ −∞

▪ If 𝐻 is a linear operator, and additivity property is

extended to integrals:
∞ ∞
𝑔 𝑥, 𝑦 = න න 𝐻[𝑓 𝛼, 𝛽 𝛿(𝑥 − 𝛼, 𝑦 − 𝛽)ሿ𝑑𝛼 𝑑𝛽
−∞ −∞

Since 𝑓 𝛼, 𝛽 is independent of 𝑥 and 𝑦.

▪ Using Homogeneity property:
∞ ∞
𝑔 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 𝐻[𝛿(𝑥 − 𝛼, 𝑦 − 𝛽)ሿ𝑑𝛼 𝑑𝛽
−∞ −∞
Where
ℎ 𝑥, 𝛼, 𝑦, 𝛽 = 𝐻[𝛿(𝑥 − 𝛼, 𝑦 − 𝛽)ሿ
Is called the impulse response of 𝐻.
▪ Ie. If 𝜂 𝑥, 𝑦 = 0, then ℎ 𝑥, 𝛼, 𝑦, 𝛽 is the response of 𝐻
to an impulse at co-ordinates 𝑥, 𝑦 .
▪ In optics, the impulse is a point of light and ℎ 𝑥, 𝛼, 𝑦, 𝛽
is commonly referred to as the point spread function
(PSF).
▪ This name arises from the fact that all physical optical
systems blur (spread) a point of light to some degree
▪ The amount of blurring is determined by the quality of
the optical components.
▪ On substitution, we get:
∞ ∞
𝑔 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 ℎ 𝑥, 𝛼, 𝑦, 𝛽 𝑑𝛼 𝑑𝛽
−∞ −∞
▪ Which is called the superposition (or Fredholm)
integral of the first kind.
▪ It states that if the response of H to an impulse is
known, the response to any input 𝑓 𝛼, 𝛽 , can be
calculated by means of above equation.
▪ A linear system H is completely characterized by its
impulse response.
▪ If H is position invariant, then:
𝐻 𝛿 𝑥 − 𝛼, 𝑦 − 𝛽 = ℎ 𝑥 − 𝛼, 𝑦 − 𝛽
▪ The earlier equation reduces to:
∞ ∞
𝑔 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 ℎ 𝑥 − 𝛼, 𝑦 − 𝛽 𝑑𝛼 𝑑𝛽
−∞ −∞
▪ This is the convolution integral for 2 variables!
▪ In presence of noise, linear degradation model
becomes: ∞ ∞
𝑔 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 ℎ 𝑥, 𝛼, 𝑦, 𝛽 𝑑𝛼 𝑑𝛽 + 𝜂(𝑥, 𝑦)
−∞ −∞
▪ In case of H being position invariant:
∞ ∞
𝑔 𝑥, 𝑦 = න න 𝑓 𝛼, 𝛽 ℎ 𝑥 − 𝛼, 𝑦 − 𝛽 𝑑𝛼 𝑑𝛽 + 𝜂(𝑥, 𝑦)
−∞ −∞
▪ The values of the noise term 𝜂(𝑥, 𝑦) are random, and are
assumed to be independent of position.
▪ Image in spatial domain is:
𝑔(𝑥, 𝑦) = ℎ(𝑥, 𝑦) ∗ 𝑓(𝑥, 𝑦) + 𝜂(𝑥, 𝑦)
▪ In frequency domain, the representation is:
𝐺(𝑢, 𝑣) = 𝐻(𝑢, 𝑣)𝐹(𝑢, 𝑣) + 𝑁(𝑢, 𝑣)
▪ Many types of degradations can be approximated by
linear, position-invariant processes.
▪ Advantage - extensive tools of linear system theory
can be used for image restoration problems.
▪ Image Restoration is also called image deconvolution
▪ Filters used in the restoration process often are called
deconvolution filters.
Estimating the
Degradation
Function
Introduction
▪ There are three principal ways to estimate the
degradation function for use in image restoration:
o Observation
o Experimentation
o Mathematical Modeling
▪ The process of restoring an image by using a
degradation function that has been estimated in some
way sometimes is called blind deconvolution.
▪ This is due to 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 degradation function H.
▪ Based on the assumption that the image was degraded by
a linear, position-invariant process, we can estimate H by
gathering information from the image itself.
o Eg: if image is blurred, look at a small rectangular section of
the image containing sample structures, like part of an object
and the background
▪ 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 is to process the subimage to arrive at a
result that is as unblurred as possible.
o Eg: Sharpen the subimage with a sharpening filter and even by
processing small areas by hand
▪ Let observed subimage be 𝑔𝑠 (𝑥, 𝑦) and processed
subimage (estimate of original image in that area) be
𝑓෡𝑠 (𝑥, 𝑦).
▪ Assuming effect of noise is negligible due to our choice
of a strong-signal area, it follows:
𝐺𝑠 𝑢, 𝑣
𝐻𝑠 𝑢, 𝑣 =
𝐹෡𝑠 (𝑢, 𝑣)
▪ From function characteristics, we can deduce complete
degradation function 𝐻(𝑢, 𝑣) based on assumption of
position invariance.
o Eg: suppose that a radial plot of 𝐻𝑠 𝑢, 𝑣 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.
▪ 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 available, it is possible to
obtain an accurate estimate of the degradation.
▪ Images similar to the degraded image can be
acquired with various system settings until they are
degraded as closely as possible to the image we wish
to restore.
▪ Then, obtain the impulse response of the degradation
by imaging an impulse (small dot of light) using the
same system settings.
▪ Practically, an impulse is simulated by a bright dot of
light, as bright as possible to reduce the effect of
noise to negligible values.
▪ Since Fourier transform of an impulse is a constant:
𝐺 𝑢, 𝑣
𝐻 𝑢, 𝑣 =
𝐴
▪ Where 𝐺 𝑢, 𝑣 is the Fourier transform of the
observed image and A is a constant describing the
strength of the impulse.
Estimation by Modeling
▪ Used for many years, due to insight we can get into
image restoration problem.
▪ In some cases, the model can even take into account
environmental conditions that cause degradations.
▪ Eg: A degradation model proposed by Hufnagel and
Stanley [1964] is based on the physical
characteristics of atmospheric turbulence. It is of the
form:
2 2 5ൗ
𝐻 𝑢, 𝑣 = 𝑒 −𝑘(𝑢 +𝑣 ) 6

▪ Where 𝑘 is a constant depending on nature of

turbulence.
▪ This is similar to the Gaussian LPF, and can be used
to model mild, uniform blurring.
▪ Another major approach - derive a mathematical
model starting from basic principles.
▪ Eg: Assume a case with image being blurred by
uniform linear motion between image and sensor
during acquisition.
▪ Suppose image 𝑓(𝑥, 𝑦) undergoes planar motion, and
𝑥0 𝑡 and 𝑦0 (𝑡)are time varying components of
motion in 𝑥 − and 𝑦 − directions.
▪ Total exposure at any point of time is obtained 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, isolates the effect of image
motion.
▪ if T is the duration of the exposure, blurred image is:
𝑇

𝑔 𝑥, 𝑦 = න 𝑓[𝑥 − 𝑥0 𝑡 , 𝑦 − 𝑦0 𝑡 ሿ𝑑𝑡
0

▪ Fourier Transform of the equations:

∞ ∞

𝐺 𝑢, 𝑣 = න න 𝑔(𝑥, 𝑦)𝑒 −𝑗2𝜋(𝑢𝑥+𝑣𝑦) 𝑑𝑥𝑑𝑑𝑦

−∞ −∞

▪ This can be written as:

𝑇 ∞ ∞

𝐺 𝑢, 𝑣 = න න න 𝑓[𝑥 − 𝑥0 𝑡 , 𝑦 − 𝑦0 𝑡 ሿ𝑒 −𝑗2𝜋(𝑢𝑥+𝑣𝑦) 𝑑𝑥𝑑𝑦 𝑑𝑡

0 −∞ −∞

▪ The term inside the outer brackets is the Fourier transform of

the displaced function𝑇𝑓[𝑥 − 𝑥0 𝑡 , 𝑦 − 𝑦0 𝑡 ሿ.
𝐺 𝑢, 𝑣 = 𝐹 𝑢, 𝑣 න 𝑒 −𝑗2𝜋 𝑢𝑥0 𝑡 +𝑣𝑦0 𝑡 𝑑𝑡 = 𝐹 𝑢, 𝑣 𝐻(𝑢, 𝑣)
0

𝑇 −𝑗2𝜋 𝑢𝑥 𝑡 +𝑣𝑦 𝑡
Where 𝐻 𝑢, 𝑣 = ‫׬‬0 𝑒 0 0 𝑑𝑡 is the transfer
function obtained from motion variables.
▪ Eg.- if image undergoes uniform linear motion in
only 𝑥 − direction at a rate 𝑥0 𝑡 = 𝑎𝑡/𝑇. When
t=T, image is displaced by a distance a.
𝑇 𝑇
𝑇
𝐻 𝑢, 𝑣 = න 𝑒 −𝑗2𝜋𝑢𝑥0 𝑡 𝑑𝑡 = න 𝑒 −𝑗2𝜋𝑢𝑎𝑡/𝑇 𝑑𝑡 = sin(𝜋𝑢𝑎)𝑒 −𝑗2𝜋𝑢𝑎
𝜋𝑢𝑎
0 0

▪ Observe that H vanishes at values of 𝑢 given by 𝑢 =

𝑛/𝑎, where 𝑛 is an integer.
▪ If we allow the y- component to vary as well. with
the motion given by 𝑥0 𝑡 = 𝑏𝑡/𝑇, then the
degradation function becomes:
𝑇
𝐻 𝑢, 𝑣 = sin[𝜋 𝑢𝑎 + 𝑣𝑏 ሿ𝑒 −𝑗2𝜋(𝑢𝑎+𝑣𝑏)
𝜋(𝑢𝑎 + 𝑣𝑏)
Inverse Filtering
Introduction
▪ Simplest approach to restoration.
෠ 𝑣) of transform of original
▪ Compute estimate, 𝐹(𝑢,
image, by dividing transform of degraded image
𝐺(𝑢, 𝑣) by the degradation function:
𝐺 𝑢, 𝑣
𝐹෠ 𝑢, 𝑣 =
𝐻(𝑢, 𝑣)
▪ Division is an array operation.
▪ Since 𝐺(𝑢, 𝑣) = 𝐻(𝑢, 𝑣)𝐹(𝑢, 𝑣) + 𝑁(𝑢, 𝑣), we get
𝑁 𝑢, 𝑣
𝐹෠ 𝑢, 𝑣 = 𝐹 𝑢, 𝑣 +
𝐻(𝑢, 𝑣)
▪ Even if we know the degradation function we cannot
recover the undegraded image exactly because
𝑁 𝑢, 𝑣 is not known.
▪ If the degradation function has zero or very small
values, then the ratio N(u, v)/ H(u, v) could easily
dominate the estimate F(u, v). This is quite common.
▪ One approach to get around the zero or small-value
problem is to limit the filter frequencies to values
near the origin.
▪ Since 𝐻(0,0) is usually the largest value in 𝐻(𝑢, 𝑣)
in frequency domain, we can limit analysis to
frequencies near the origin, to reduce probability of
zero values.
▪ Direct inverse filtering delivers poor performance in
general.
Minimum Mean
Square Error
(Wiener) Filtering
Introduction
▪ Inverse filtering approach makes no explicit provision for
handling noise.
▪ This method incorporates both degradation function and
statistical characteristics of noise in restoration.
▪ We consider noise as random variables, and objective is
to find an estimate 𝑓መ of uncorrupted image 𝑓, such that
mean square error between them is minimized.
▪ Error measure is given by:
መ 2
𝑒 2 = E (𝑓 − 𝑓)
where 𝐸{. } is the expected value of the argument.
▪ Assumptions:
o Noise and image are uncorrelated
o One or the other of them has zero mean
o Intensity levels in the estimate are linear functions of degraded
image intensity levels.
▪ Based on conditions, minimum of error functions is
given in frequency domain by:
𝐻 ∗ 𝑢, 𝑣 𝑆𝑓 𝑢, 𝑣
𝐹෠ 𝑢, 𝑣 = 𝐺 𝑢, 𝑣
𝑆𝑓 𝑢, 𝑣 𝐻 𝑢, 𝑣 2 + 𝑆𝜂 𝑢, 𝑣

𝐻 ∗ 𝑢, 𝑣
= 𝐺 𝑢, 𝑣
𝐻 𝑢, 𝑣 2 + 𝑆𝜂 𝑢, 𝑣 /𝑆𝑓 𝑢, 𝑣

1 𝐻 𝑢, 𝑣 2
= 2 + 𝑆 𝑢, 𝑣 /𝑆 𝑢, 𝑣
𝐺 𝑢, 𝑣
𝐻 𝑢, 𝑣 𝐻 𝑢, 𝑣 𝜂 𝑓

▪ This result is known as the Wiener filter, after N.

Wiener who first proposed the concept in 1942.
▪ Also called minimum mean square error filter or the
least square error filler.
▪ Note from equation that the Wiener filter does not
have the same problem as the inverse filter with
zeros in the degradation function, unless the entire
denominator is zero for the same value(s) of u and v.
▪ The terms are as follows:
o 𝐻 𝑢, 𝑣 - Transformation of Degradation Function
o G 𝑢, 𝑣 - Transformation of Degraded Image
o 𝐻 ∗ 𝑢, 𝑣 - Complex Conjugate of Degradation Function
o 𝐻 𝑢, 𝑣 2 = 𝐻 ∗ 𝑢, 𝑣 𝐻 𝑢, 𝑣
o 𝑆𝜂 𝑢, 𝑣 = 𝐻 𝑢, 𝑣 2 - Noise Power Spectrum
o 𝑆𝑓 𝑢, 𝑣 = 𝐹 𝑢, 𝑣 2 - Undegraded Image Power
Spectrum
▪ If the noise is zero. then the noise power spectrum
vanishes and the Wiener filter reduces to the inverse
filter.
Important measures
1. Signal to Noise Ratio
σ𝑀−1 𝑁−1
𝑢=0 σ𝑣=0 𝐹(𝑢, 𝑣)
2
𝑆𝑁𝑅 = 𝑀−1 𝑁−1
σ𝑢=0 σ𝑣=0 𝑁(𝑢, 𝑣) 2

▪ Measure of the level of information bearing signal

power (i.e., of the original, undegraded image) to the
level of noise power.
▪ Images with low noise tend to have a high SNR and,
image with a higher level of noise has a lower SNR.
▪ By itself, SNR has limited value, but is an important
metric to characterize performance of restoration
algorithms.
2. Mean Square Error
𝑀−1 𝑁−1
1 2
𝑀𝑆𝐸 = መ 𝑦)
෍ ෍ 𝑓 𝑥, 𝑦 − 𝑓(𝑥,
𝑀𝑁
𝑢=0 𝑣=0
▪ Restored image can be considered as the signal and
difference between restored and original images is
considered noise.
▪ SNR in spatial domain:
σ𝑀−1 σ 𝑁−1 መ 2
𝑥=0 𝑦=0 𝑓(𝑥, 𝑦)
𝑆𝑁𝑅 = 2
𝑀−1 𝑁−1
σ𝑥=0 σ𝑦=0 𝑓 𝑥, 𝑦 − 𝑓(𝑥,መ 𝑦)
▪ The closer 𝑓 and 𝑓መ are, the larger the ratio will be.
▪ If square root of the measures is used, we call them
root-mean-square-signal-to-noise ratio and the root-
mean-square-error, respectively.
▪ Quantitative metrics do not necessarily relate well to
perceived image quality.
▪ When we are dealing with spectrally white noise. the
spectrum 𝑁(𝑢, 𝑣) 2 is a constant, which simplifies
things considerably.
▪ Power spectrum of undegraded image is seldom
known.
▪ Frequently used approach when quantities are not
known or cannot be estimated for approximation:
1 𝐻 𝑢, 𝑣 2
𝐹෠ 𝑢, 𝑣 = 2
𝐺 𝑢, 𝑣
𝐻 𝑢, 𝑣 𝐻 𝑢, 𝑣 + 𝐾

Where 𝐾is a specified constant added to all terms of

𝐻 𝑢, 𝑣 2
Constrained Least
Squares Filtering
Introduction
▪ All restoration methods need some knowledge of
degradation function H.
▪ In case of Wiener filter, both power spectra of
undegraded image AND noise must be known.
▪ Good results can be obtained from earlier models,
but constant estimate of the ratio of the power
spectra is not always a suitable solution.
▪ In this section, we have a method which requires
knowledge of only the mean and variance of the
noise, which can usually be calculated from a given
degraded image.
▪ Unlike Wiener filter, which is optimal in the average
sense, this method yields an optimal result for each
image to which it is applied.
▪ Using convolution definition and expressing
degradation model equation in matrix form:
𝒈=𝑯𝒇+𝜼
o Eg: Suppose 𝑔(𝑥, 𝑦) is of size (MxN). Then we can form
the first N elements of vector 𝒈 by using image elements in
first row of 𝑔(𝑥, 𝑦), next N elements from the second row
and so on.
o Resulting vector has size (MNx1), similar to 𝒇 and 𝜼.
o Matrix 𝑯 has the dimensions (MNxMN)
▪ Thus, restoration problem can now be reduced to
matrix manipulations.
▪ However, due to large sizes of M and N, the problem
becomes complicated.
▪ Also, 𝑯 is highly sensitive to noise, which
complicates the situation
▪ One way to reduce the noise sensitivity problem is to
base optimality of restoration on a measure of
smoothness, such as the second derivative of an image
(Laplacian)
▪ Objective is to find minimum of criterion function C:

𝑀−1 𝑁−1

𝐶 = ෍ ෍ 𝛻 2 𝑓(𝑥, 𝑦) 2

𝑥=0 𝑦=0
Subject to constraint:
2
෠
𝒈 − 𝑯𝒇 = 𝜼 2
Where 𝒘 2 ≜ 𝒘𝑇 𝒘 is the Euclidean norm, and 𝒇෠ is the
estimate of undegraded image.

𝜕 2𝑓 𝜕 2𝑓
𝛻 2 𝑓 𝑥, 𝑦 = 2 + 2
𝜕𝑥 𝜕𝑦
▪ Solution in frequency domain:
𝐻 ∗ 𝑢, 𝑣
𝐹෠ 𝑢, 𝑣 = 𝐺 𝑢, 𝑣
𝐻 𝑢, 𝑣 2 + 𝛾 𝑃(𝑢, 𝑣) 2

Where 𝛾 is a parameter to be adjusted so that

constraint is satisfied.
▪ 𝑃(𝑢, 𝑣) is the Fourier Transform of the function:
0 −1 0
𝑝 𝑥, 𝑦 = −1 4 −1
0 −1 0
▪ This is the Laplacian operator in spatial domain.
▪ 𝑝 𝑥, 𝑦 and other spatial domain functions should be
padded with zeros as required
▪ The equation reduces to inverse filtering if 𝛾 = 0
▪ It is possible to adjust the parameter 𝛾 interactively
until acceptable results are achieved.
▪ If optimality is to be achieved, 𝛾 must be adjusted so
that constraint is satisfied.
▪ Procedure for computation:
o Define “residual” vector 𝒓 as:
𝒓 = 𝒈 − 𝑯𝒇෠
෠ 𝑣) is a function of 𝛾, then 𝒓 is also a function of
o Since 𝐹(𝑢,
this parameter.
o It can be shown that:
𝜙 𝛾 = 𝒓𝑇 𝒓 = 𝒓 𝟐
is a monotonically increasing function of 𝛾
o We have to adjust 𝛾 so that
𝒓 𝟐 = 𝜼 𝟐±𝑎
Where 𝑎 is an ACCURACY FACTOR
▪ Since 𝜙 𝛾 is strictly monotonic, finding the desired
value of 𝛾 is easy.
▪ Approach is as follows:
1. Specify an initial value of 𝛾
2. Compute 𝒓 𝟐
3. Stop if 𝒓 𝟐 = 𝜼 𝟐 ± 𝑎 is satisfied; otherwise return to
step 2 after:
➢ Increasing 𝛾 if 𝒓 𝟐
< 𝜼 𝟐 −𝑎
OR
➢ decreasing 𝛾 if 𝒓 𝟐 > 𝜼 𝟐
−𝑎.
4. Use the new value of 𝛾 in solution equation, to recompute
the optimum estimate 𝐹෠ 𝑢, 𝑣 .

Other procedures. such as a Newton-Raphson algorithm, can

be used to improve the speed of convergence.
▪ In order 𝟐to use the algorithm we need the quantities 𝒓 𝟐
and 𝜼 .
▪ To compute 𝒓 𝟐 we take inverse of 𝑅(𝑢, 𝑣) from
Residual Function: 𝒓 = 𝒈 − 𝑯𝒇: ෠

𝑅 𝑢, 𝑣 = 𝐺 𝑢, 𝑣 − 𝐻(𝑢, 𝑣)𝐹෠ 𝑢, 𝑣

𝑀−1 𝑁−2
2
𝒓 = ෍ ෍ 𝑟 2 (𝑥, 𝑦)
𝑥=0 𝑦=0
▪ Computation of 𝜼 𝟐 is done by considering
2
noise
variance
2
over the entire image (𝜎𝑛 ) and sample mean
(𝑚𝑛 ) as:
𝜼 𝟐 = 𝑀𝑁 𝜎𝑛2 − 𝑚𝑛2
▪ Thus, optimum restoration algorithm can be
implemented with knowledge of only mean and variance
of noise.
▪ Optimum restoration in the sense of constrained least
squares does not necessarily imply "best" in the
visual sense.
▪ Other parameters in the algorithm for iteratively
determining the optimum estimate also play a role in
the final result.
▪ In general, automatically determined restoration
filters yield inferior results to manual adjustment of
filter parameters.
▪ Particularly true of the constrained least squares
filter, which is completely specified by a single
scalar parameter.