Image Restoration
Jayanta Mukhopadhyay
Dept. of CSE,
IIT Kharagpur
Image smoothing / enhancement:
Gaussian Filtering
"#!(%! &'! ))!
!
l Low pass filtering 𝐻 𝑢, 𝑣 = 𝐴𝑒 *
l Gaussian filter
l Freq. / Spatial domain
(+ &, ! !)
𝐴 !
ℎ 𝑚, 𝑛 = 𝑒 *)!
4𝜋𝜎 *
l More complex filters
l Using difference of Gaussian
l LPF, BPF, HPF "# (%! &'! ))"!
!
"# (%! &'! ))!!
!
! !
𝐻 𝑢, 𝑣 = 𝐴1𝑒 ! - 𝐴 2𝑒 !
l DoG in spatial domain, too!
!
(* &+ !) !
(* &+ !)
-" ! -! !
ℎ 𝑚, 𝑛 = "#)"! 𝑒 !),! - "#)!! 𝑒 !)!!
Filtering in 1-D the transform domain
/
.*# 0#,
• Freq. Shifting 𝐷𝐹𝑇(𝑓(𝑛)𝑒 ) = 𝐹(< 𝑘 − 𝑘0 > 𝑁)
/
!.*#0,#
• Phase Shifting 𝐼𝐷𝐹𝑇(𝐹(𝑘)𝑒 ) = 𝑓(< 𝑛 − 𝑛0 > 𝑁)
F(k) H(k) G(k)=H(k) F(k)
l Use sufficient 0 padding at the both end to make circular
convolution equivalent to linear convolution
l To take care of boundary effect.
l The length of f(n) and h(n) should be the same.
l H(k) usually provided as symmetric about the center.
l 0th freq. at the N/2 th element.
l Center F(k) as Fc(k) by multiplying f(n) with (-1)n ß (k0=N/2)
l Obtain G(k)= H(k) . Fc(k)
l Multiply G(k) by (-1)k and perform IDFT to get g(n).ß (n0=N/2)
Low pass filters
l Ideal LPF
l H(u,v) =1 , if D(u,v) < D0
= 0, Otherwise
l Use centered freq. transform:
l D (u,v) = [(u – M/2)2+ (v –N/2)2]1/2
l Make DFT of image also centered
l Perform filtering
l Not very practical
l Blurring and ringing effect prominent!
l Sharp discontinuity in frequency response!
l Impulse response in the form of a Sinc function.
4
Butterworth Low Pass Filters
l Transfer function of order n: 1
𝐻 𝑢, 𝑣 = *,
l Avoids sharp discontinuity 𝐷(𝑢, 𝑣)
l 50% of maximum as D(u,v)=D0
1+ 𝐷1
• First order: No ringing
• Second order: Ringing imperceptible
• Higher the order ringing effect higher.
Usually BPLF of order 2
recommended.
Similarly we may have
HPF and BPF.
5
Courtesy: R.C. Gonzalez and R.E Woods © 1992-2008
Homomorphic filtering
l Let f(x,y)=i(x,y) r(x,y)
l i(x,y): Illumination variation
l Slow / Low frequency
l r(x,y): reflectance variation
l High / High frequency
l Make additive component in the log domain
l Apply filtering by suppressing low frequencies
and enhancing high frequencies.
l Dynamic range compression coupled with contrast
enhancement
l Back to the original domain by exponentiation. 6
Homomorphic filtering
5(%,')!
!4
l A typical example 𝐻 𝑢, 𝑣 = 𝛾2 − 𝛾3 1−𝑒 5- ! + 𝛾3
l 𝛾L<1, 𝛾H>1
log(.) DFT
f(x,y)
H(u,v)
g(x,y) exp(.) IDFT
7
Homomorphic filtering: Typical
example
l Simultaneous dynamic range compression and contrast
enhancement
5(%,')!
!4 !
𝐻 𝑢, 𝑣 = 𝛾2 − 𝛾3 1−𝑒 5 - + 𝛾3
Courtesy: R.C. Gonzalez and R.E Woods © 1992-2008
Band reject filters BPF=1-BRF
(BRF)
Distance from the center
𝑊 𝑤
l Ideal 𝐻 𝑢, 𝑣 = < 0 − 2 + 𝐷0 ≤ 𝐷 𝑢, 𝑣 ≤ 𝐷0 + 2
1 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Cut-off frequency
1
l Butterworth 𝐻 𝑢, 𝑣 = *,
𝑊. 𝐷(𝑢, 𝑣)
1+ *
𝐷 𝑢, 𝑣 − 𝐷1*
!
5! %,' !5-!
!
Gaussian 75(%,')
l 𝐻 𝑢, 𝑣 = 1 − 𝑒
9
Used for removing periodic noise.
Notch filters
l A notch filter rejects or passes frequencies in a predefined
frequency about the center of the frequency rectangle.
l Zero phase shift filter should be symmetric about the
origin.
l If there exists a notch at center (u0,v0), there must be a notch
at (-u0,-v0). ;
A general form:
l
𝐻08 𝑢, 𝑣 = I 𝐻/ (𝑢, 𝑣) 𝐻!/ (𝑢, 𝑣)
/9:
l e.g. using Butterworth HPFs
HPF with center at (uk,vk) or (u-k,v-k)
1
𝐻! 𝑢, 𝑣 = #$ %
𝐷"! 𝑀 #
𝑁 # #
1+
𝐷! (𝑢, 𝑣) 𝐷! 𝑢, 𝑣 = 𝑢− − 𝑢! + 𝑣− − 𝑣!
2 2
A degradation model and restoration
l Degradation function: l Linear degradation model.
l Identity (No degradation), 𝑔 𝑥, 𝑦 = ℎ 𝑥, 𝑦 ∗ 𝑓 𝑥, 𝑦 + 𝑛 (𝑥, 𝑦)
Linear filter (Motion Blur),
Transformation of a l To design a restoration filter w(x,y)
functional value l such that w(x,y)*g(x,y) = fa(x,y) close
l Noise model to f(x,y).
l White, Gaussian, l To minimize
Rayleigh, 𝐸( 𝑓 𝑥, 𝑦 − 𝑓" (𝑥, 𝑦) # )
f(x,y) Degradation g(x,y) Restoration fa(x,y)
function
H(.)
+ filters
W(.)
If H(.) is the identity Objective function in Freq. domain?
function, the task is 𝐸( 𝐹 𝑢, 𝑣 − 𝐹G (𝑢, 𝑣) *)
simply noise cleaning.
n(x,y)
Restoration in the absence of noise
𝑔 𝑥, 𝑦 = ℎ 𝑥, 𝑦 ∗ 𝑓 𝑥, 𝑦 G(u,v)=H(u,v) F(u,v)
Fa(u,v)
F(u,v) G(u,v) W(u,v)
H(u,v) ?
l Inverse filtering: W(u,v) = 1/H(u,v)
l Problem with zeros and low values in H(u,v)
l Minimize E(||F(u,v) - Fa(u,v) ||2)
= E(||F(u,v) - W(u,v) H(u,v) F(u,v) ||2)
12
Restoration in the absence of noise
Fa(u,v)
F(u,v) G(u,v) W(u,v)
H(u,v) ?
l Power Spectrum of the image: Sf(u,v)=||F(u,v)||2=F*(u,v)F(u,v)
l To minimize Ew=E(||F(u,v) - W(u,v) H(u,v) F(u,v) ||2)
For convenience, F(u,v) written as F and the same for all others.
Ew=||F||2- (W*H*+WH)||F||2+||W||2||H||2||F||2
𝜕𝐸7
=0 -HSf+W*||H||2Sf=0 W=H*/||H||2=1/H
𝜕𝑊(𝑢, 𝑣)
Pretending W* Inverse filtering! 13
constant for W The same problem!
Restoration in the presence of
noise F (u,v) a
G(u,v)=H(u,v)F(u,v)+N(u,v) W(u,v)
?
l Power Spectrum of the original image: Sf= ||F||2 =F*F
l Power Spectrum of the noisy image: Sg= ||G||2 =G*G
l Sg = ||H||2 ||F||2 +H*F*N+HFN*+ ||N||2
l As noise assumes to be uncorrelated: E(F*N)=E(N*F)=0
l Hence, Sg = ||H||2 ||F||2 +||N||2 = ||H||2 Sf + Sn
l To minimize Ew=E(||F – WG||2) W=(H*Sf )/ (||H||2Sf +Sn)
l Ew= E(||F ||2-FW*G*-WGF*+ ||W||2 ||G||2 )
l = E(||F ||2-FW*G*-WHF*F-WF*N+ ||W||2 ||G||2 )
𝜕𝐸7
=0 -HSf+W* (||H||2Sf +Sn)=0
𝜕𝑊(𝑢, 𝑣)
Weiner Filter (Least square error
filter) n(x,y)
f(x,y) g(x,y) fa(x,y)
h(x,y) + w(x,y)
l Solution in frequency domain: W=(H*Sf )/ (||H||2Sf +Sn)
𝐻∗ 𝐻∗
𝑊= 𝑊=
𝐻 #+𝐾 # 𝑆%
𝐻 +
Noise to Signal Ratio 𝑆&
1 𝐻 #
𝑊= Weighted Inverse filter!
𝐻 𝐻 #+𝐾
Tasks of restoration 𝐻∗
𝑊=
# 𝑆%
𝐻 +𝑆
n(x,y) &
f(x,y) g(x,y) fa(x,y)
h(x,y) + w(x,y)
l Model degradation:
l Design / derive h(x,y)
l Model Noise
l Identify PDF and estimate parameters
l Derive W or w(x,y)
l Apply filtering: fa(x,y)= w(x,y)* g(x,y)
Degradation of a defocused
image
l Defocused image: The projected point is not sharp.
l e.g. The projection forms a circle of radius r.
l Spatial resolution along x: ∆x and y: ∆y
N: Total number of
pixels within the circle
h(i,j)
h(i , j)= 1/N, (i. ∆x )2+(j. ∆y )2 < r
= 0 Otherwise
Degradation due to motion blur
l Motion blur: Movement of camera, Movement of
object.
l e.g. The camera moving with a velocity vx in the
direction of x
l Exposure time: t
l Spatial resolution along x: ∆x
l Number of pixels covered in a shot due to movement:
N=(vxt) / ∆x
h(i,j)= 1/N, -(N-1) < i < 0
= 0 Otheriwse
h(I,j)
Degradation model for a general
planar rigid body motion
l Time varying translational motion vector : x0(t) i + y0(t) j
l Duration of exposure: T
l Shutter opening and closing takes place instantaneously
l Blurred image g(x,y) given by:
(
𝑔 𝑥, 𝑦 = 3 𝑓(𝑥 − 𝑥' 𝑡 , 𝑦 − 𝑦' 𝑡 ) 𝑑𝑡
'
Fourier transform of g(x,y):
(
𝐺 𝑢, 𝑣 = 3 𝐹(𝑢, 𝑣)𝑒 )*#+(-.! / 012! / )
𝑑𝑡 19
'
Degradation model for a general
planar rigid body motion
(
𝑔 𝑥, 𝑦 = 3 𝑓(𝑥 − 𝑥' 𝑡 , 𝑦 − 𝑦' 𝑡 ) 𝑑𝑡
'
(
𝐺 𝑢, 𝑣 = 3 𝐹(𝑢, 𝑣)𝑒 )*#+(-.! / 012! / )
𝑑𝑡
'
(
𝐺 𝑢, 𝑣 = 𝐹 𝑢, 𝑣 [3 𝑒 )*#+ -.! / 012! /
𝑑𝑡 ]
'
H(u,v)
Motion blur special cases
l Uniform linear motion in horizontal direction:
l x0(t)=at/T, and y0(t)=0
(
𝐻 𝑢, 𝑣 = 3 𝑒 )*#+ -.! / 012! /
𝑑𝑡
'
( -"/
)*#+
𝐻 𝑢, 𝑣 = 3 𝑒 ( 𝑑𝑡
'
𝑇
= sin(𝜋𝑢𝑎)𝑒 )*+-"
𝜋𝑢𝑎 21
Motion blur special cases
l Uniform linear motion in a plane:
l x0(t)=at/T, and y0(t)= bt/T
(
𝐻 𝑢, 𝑣 = 3 𝑒 )*#+ -.! / 012! /
𝑑𝑡
'
( -"/ 14/
)*#+ 0
𝐻 𝑢, 𝑣 = 3 𝑒 ( ( 𝑑𝑡
'
𝑇
= sin(𝜋(𝑢𝑎 + 𝑣𝑏))𝑒 )*+(-"014)
𝜋(𝑢𝑎 + 𝑣𝑏) 22
Atmospheric degradation model
(Hufnagel & Stanley (1964))
#
)5(- " 01 " )$
𝐻 𝑢, 𝑣 = 𝑒
l k is a constant that depends on the turbulence.
l Almost the form of a Gaussian function
l LPF
How do you relate them in discrete domain?
Use normalized frequency: k/N for sequence of length N.
23
Variance is the measure of noise power.
White noise with a pdf has the same power
Noise models at every freq.
1 6 .)7 "
)
l Gaussian: 𝑝 𝑥 = 𝑒 # 8
2𝜋𝜎
𝑚𝑒𝑎𝑛 = 𝜇 𝑣𝑎𝑟 = 𝜎*
2 .)" "
) &:; .<"
𝑝 𝑥 = H𝑏 (𝑥 − 𝑎)𝑒 4
l Rayleigh:
𝜋𝑏
0 𝑓𝑜𝑟 𝑥 < 𝑎
𝑚𝑒𝑎𝑛 = 𝑎 +
4
𝑏(4 − 𝜋)
𝑣𝑎𝑟 =
4 𝑎4 𝑥 4)6 )".
Erlang 𝑒 𝑓𝑜𝑟 𝑥 ≥ 0
l 𝑝 𝑥 =H 𝑏−1 !
(Gamma): 𝑏 0 𝑓𝑜𝑟 𝑥 < 𝑎
𝑏
𝑚𝑒𝑎𝑛 = 𝑣𝑎𝑟 = 𝑎*
24
𝑎
Noise models
l Exponential: 𝑎𝑒 )". 𝑓𝑜𝑟 𝑥 ≥ 0
1 𝑝 𝑥 = O
1
𝑣𝑎𝑟 = 0 𝑓𝑜𝑟 𝑥 < 0
𝑚𝑒𝑎𝑛 =
𝑎 𝑎*
1
l Uniform: 𝑝 𝑥 = P𝑏 − 𝑎 𝑓𝑜𝑟 𝑎 ≤ 𝑥 ≤ 𝑏
𝑎+𝑏
𝑚𝑒𝑎𝑛 =
2
0 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
𝑏−𝑎 *
𝑣𝑎𝑟 =
12 𝑃" 𝑓𝑜𝑟 𝑥 = 𝑎
l Impulse (Salt 𝑝 𝑥 = P 𝑃4 𝑓𝑜𝑟 𝑥 = 𝑏
and Pepper): 0 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Bipolar impulse 25
Estimation of noise
l Study relatively flat (constant) region and study the
histogram.
l The shape of histogram may indicate appropriate
PDF to be chosen.
l Compute mean and variance of the flat region.
l Relate to them to the parameters of the
distribution. (Optional)
l Make robust estimation of variance (by sampling
and estimation using the central limit theorem)
l Use variance for the noise power at every
frequency (modeling as white noise) 26
Noise removal: linear and nonlinear
filters exploiting local statistics
0!:
l Arithmetic mean. 1 T 𝑥 l Order Statistics.
H
𝑁 l Median
H91
:
0!: l Max
0
l Geometric mean. I 𝑥H l Min
H91 l Mid-point
l Harmonic mean. 1 l Average of Max and Min.
1 0!: 1 l Alpha-trimmed mean
∑
𝑁 H91 𝑥H l Mean excluding top (d/2)
and bottom (d/2) in the
l Contraharmonic ∑ 𝑥H;&:
rank order.
mean. ∑ 𝑥H;
Q: Order of filter 27
Q= 0 (A.M.), Q= -1 (H.M.)
Adaptive filter for restoration
l Exploit local statistics
l Local mean: 𝛍 Local variance: 𝛔2
l Local noise variance: 𝜂2
l Pixel value: g(x,y)
l Desirable
l If 𝜂2=0 return g(x,y)
l If 𝛔2 high return close to g(x,y)
l If 𝜂2 = 𝛔2 return local mean 𝛍
l Adaptive expression
𝜂#
𝑓" 𝑥, 𝑦 = 𝑔 𝑥, 𝑦 − 𝟐 (𝑔 𝑥, 𝑦 − 𝝁) 28
𝝈
𝐻∗
Constrained least 𝑊=
# 𝑆%
𝐻 +𝑆
squares filtering &
l Limitations of Wiener Filtering:
l Power spectra of the undegraded image and
noise required to be known.
l Based on minimizing a statistical criterion.
l Optimization in average sense
l Problem statement reformulated for taking care of
individual image separately.
l Criteria set from filtering in spatial domain
l But solution obtained in the frequency domain
29
Constrained least squares
filtering in Spatial domain
l Consider degraded and undegraded image in the form of
a derasterized vector
l g and f of length N, respectively.
l Let the degradation filter response given by a vector
l h of length N
l The noise at each point expressed by a vector
l 𝜂 of length N
l The convolution operation can be expressed by a matrix
operation as follows:
l g=Hf+ 𝜂
l H is a matrix of dimension NxN 30
Optimization problem
l Minimize a smoothness criteria C (measured as sum of
squares of Laplacians) subject to keeping noise at a
desired level
J!: ;!:
𝐶 = T T 𝛻 *𝑓(𝑥, 𝑦) *
I91 K91
Subject to the constraint A small
positive
𝐠 − 𝐇𝐟 * = 𝛈 *
value
or *
𝑎𝑏𝑠( 𝐠 − 𝐇𝐟 − 𝛈 * ) < 𝑎2
31
Residual difference
Frequency domain solution
𝐻 ∗ (𝑢, 𝑣)
𝐹 𝑢, 𝑣 = * * 𝐺(𝑢, 𝑣)
|𝐻 𝑢, 𝑣 | + 𝛾|𝐿 𝑢, 𝑣 |
l 𝛾: a parameter, may be chosen iteratively such that the
residual difference < a2
l L(u,v) is the Fourier Transform of the Laplacian operator
l(x,y)
0 −1 0
𝑙 𝑥, 𝑦 = −1 4 −1
0 −1 0 32
Summary
l Filtering for noise removal
l Degradation model
l Weiner (LSE) filter
l to model degradation filter and noise, and then obtain the
LSE filter.
l Apply Weiner filter on degraded image to restore it.
l Modeling motion blur and defocusing.
l Noise model
l Use of various local statistics for removing noise.
l Adaptive filter exploiting local statistics and variance of
noise.
l Constrained least squares filters
Thank You
34
