Transcript Document 7570087

Chapter 6
Image Restoration
Digital Image Processing
Instructor: Dr. Cheng-Chien Liu
Department of Earth Sciences
National Cheng Kung University
Last updated: 10 October 2003
Introduction
 Image restoration
• Use objective criteria and prior knowledge
• cf. image enhancement  subjective criteria
 Two cases need image restoration
• Degradation  gray value altered
• Distortion  pixel shifted
Geometric restoration (image registration)
Aerial photographs
Geometric restoration
 Source of geometric distortion
• Lens (Fig 6.1)
• Irregular movement (Fig 6.2)
 Two-stage operation
• Spatial transformation
x^ = Ox(x, y) = c1x + c2y + c3xy + c4
y^ = Oy(x, y) = c5x + c6y + c7xy + c8
Four tie points  c1, …, c8
• Grey level interpolation
Simple way: g(x^, y^) = ax^ + by^ + cx^y^ + d
Fig 6.3
 Example 6.1
Linear degradation
 Output image
• g(a, b) =  - - f(x, y) h(x, y, a, b) dxdy
• The point spread function: h(x, y, a, b)
 Shift invariant
• h(x, y, a, b) = h(x, y, a - x, b - y)
• g(a, b) =  - - f(x, y) h(x, y, a - x, b - y) dxdy
g depends on the relative position rather than actual position
• G(u, v) = F(u, v) H(u, v)
 For discrete images
• g(i, j) = Sk=1NSl=1Nf(k, l)h(k, l, i, j)
• g=Hf
The point spread function H
 Problems of image restoration: g = H f
• Given the degraded image g, recover the original
undegraded image f
• Obtain the information of H
From the knowledge of the physical process
 e.g.diffraction, atmospheric turbulence, motion, …
From some known objects on the image
 Example 6.2
• Expression of blurred image
 Example 6.3
• Derive H for the blurred image
The point spread function H (cont.)
 Example 6.4
• Calculate H for the blurred image
 Example 6.5
• Derive H for the degradation process of
accelerating motion
 Example 6.6
• Asymptotic solution of Example 6.5
 Example 6.7
• Application of Example 6.6
The point spread function H (cont.)
 Example 6.8
• Calculate H from a bright straight line
 Example 6.9
• Calculate H from an edge
 Example 6.10
• Calculate H from an image device
Straightforward solution
 If H is known
• F(u, v) = G(u, v) / H(u, v)
• F(u, v)  f(u, v)
 However
• Straightforward solution  unacceptable poor
results
H(u, v) = 0 at some points  G = 0  0/0  undetermined
 If there is a small amount of noise  G  0, even if H = 0
For additive noise: G(u, v) = F(u, v) H(u, v) + N(u, v)
 F(u, v) = G(u, v) / H(u, v) - N(u, v) / H(u, v)
 If H(u, v)  0  N(u, v) / H(u, v)   (amplified noise)
Straightforward solution (cont.)
 Avoiding the amplification of noise
• Windowed version of the filter 1 / H
F(u, v) = M(u, v) G(u, v) - M(u, v) N(u, v)
where M(u, v) = 1 / H(u, v) for u2 + v2  w02
M(u, v) = 1
for u2 + v2 > w02
 Where w0 is chosen so that all zeroes of H(u, v) are excluded
• Other windowing filters are also valid
 Example 6.11
• Application of inverse filtering to restore a
motion blurred image
Indirect solution – Wiener filter
 Formal expression of the problem of IR
• To identify f(r) which minimizes e2  E{[f(r) - f(r)]2}
Where f(r) is an estimate of the original undegraded image f(r)
• Shift invariant assumption
g(r) =  - - h(r - r΄) f(r΄)dr΄ + v(r)
 Where g(r), f(r) and h(r) are random fields, v(r) is noise field
 Solution  find the Wiener filter
• If no imposed condition  conditional expectation 
simulated annealing  beyond our scope
• Constraint: f(r) is a linear function of g(r)
f(r) =  - - m(r, r΄) g(r΄)dr΄  we decide (B6.1)
f(r) =  - - m(r - r΄) g(r΄)dr΄  if the random fields are homogeneous
Identify the Wiener filter m(r) with which to convolve g(r΄)  f(r)
Fourier transfer of the Wiener filter
 M(u, v) = F{m(r)} = Sfg(u, v) / Sgg(u, v)
• Proof in B6.3
• Sfg(u, v) is the cross-spectral density of f and g
• Sgg(u, v) is the spectral density of g
 Extra assumption:
• f(r) and v(r) are uncorrelated
• E{v(r)} = 0
•  E{f(r)v(r)} = E{f(r)}E{v(r)} = 0
Fourier transfer of the Wiener filter
(cont.)
 Create Sgf
• g(r) =  - - h(r - r΄) f(r΄)dr΄ + v(r)
• Rgf (s) = E{g(r)f(r - s)}
=  - - h(r - r΄) E{f(r΄)f(r - s)}dr΄ + E{f(r - s)v(r)}
=  - - h(r - r΄) Rff(r΄ - r + s)dr΄
• Sgf(u, v) = H*(u, v)Sff(u, v)  (B6.4)
• Sgg(u, v) = Sff(u, v)|H(u, v)|2 + Svv(u, v)  (B6.4)
 M(u, v)
• M = H*Sff / [Sff|H|2 + Svv]
• M = (1/H) |H|2 / [|H|2 + Svv/Sff ]
Fourier transfer of the Wiener filter
(cont.)
 Noise
• If there is no noise  Svv(u, v) = 0  M = 1/H
So the linear least square error approach simply
determines a correction factor with which the inverse
transfer function of the degradation process has to be
multiplied before it is used as a filter, so that the effect of
noise is taken care of.
• Assumption
White noise:
Svv(u, v) = constant = Svv(0, 0) =  - - Rvv(x, y)dxdy
Ergodic noise: Rvv(x, y) can be obtained from a single pure
noise image i.e. when f(x, y) = 0
Fourier transfer of the Wiener filter
(cont.)
 B6.1
• If m(r - r΄) satisfies E{[f(r) -  - - m(r - r΄)
g(r΄)dr΄]g(s)} = 0, then it minimizes e2  E{[f(r) - f(r)]2}
 Example 6.12
• g(r) =  - - h(t - r) f(t)dt  G(u, v) = H*(u, v) F(u, v)
 B6.2
• Wiener-Khinchine theorem: Rff(u, v) = |Ffg(u, v)|2
 B6.3
• M(u, v) = F{m(r)} = Sfg(u, v) / Sgg(u, v)
Fourier transfer of the Wiener filter
(cont.)
 B6.4
• Sgg(u, v) = Sff(u, v)|H(u, v)|2 + Svv(u, v)
 Example 6.13
• Apply Wiener filtering to restore a motion
blurred image
Problems of the straightforward
solution
 Straightforward solution
• g = Hf
• Including noise: g = Hf + v
• Inversion: f = H-1g – H-1v
H is an N2  N2 matrix
f, g and v are N2  1 vectors
• Problems
f is very sensitive to v (Example 6.14)
Formidable task to inverse an N2  N2 matrix
Circulant matrix
 Definition
• The circulant matrix D (Eq. 6.78)
Each column of a matrix can be obtained from the precious one by
shifting all elements one place and putting the last element at the top
• The block circulant matrix (Eq. 6.77)
 Dw(k) = l(k)w(k)
• l(k) are the eigenvalues of D
l(k)  d(0) + d(M-1)exp[2pjk/M] + d(M-2)exp[2pj2k/M] + … +
d(1)exp[2pj(M-1)k/M]
• w(k) are the eigenvectors of D
w(k)  [1, exp[2pjk/M], exp[2pj2k/M], …, exp[2pj(M-1)k/M]]T
Inversion of the circulant matrix
 Inversion of D
• D = WLW-1
W is formed by having the eigenvectors of D as columns
W-1(k, j) = (1/M)exp[-2pj/M ki] (Example 6.15)
 L is a diagonal matrix with the eigenvalues alone its diagonal.
•
•
•
•
D-1 = (WLW-1)-1 = (W-1)-1L-1W-1 = WLW-1
Example 6.16: A case of M = 3
Example 6.17: A case of M = 4
Example 6.18:
W  WN  WN  W-1 = Z  WN-1  WN-1
 WN (k, n) = (N)-1/2 exp[2pj/N kn]
 WN-1 (k, n) = (N)-1/2 exp[-2pj/N kn]
 Kronecker product 
Inverting H – Overcome one problem
of the straightforward solution
 H is block circulant
•
•
•
•
g=Hf
g(i, j) = Sk=0N-1Sl=0N-1h(k, l, i, j) f(k, l)
For a shift invariant point spread function
g(i, j) = Sk=0N-1Sl=0N-1f(k, l) h(i-k, j-l)
 Diagonalize H
• H = WLW-1 (B 6.5)
WN (k, n) = (N)-1/2 exp[2pj/N kn]
WN-1 (k, n) = (N)-1/2 exp[-2pj/N kn]
L(k, i) = NH(kmod N, [k/N]) if i = k
L(k, i) = 0 if i  k
H(m,n) = (1/N) Sx=0N-1Sy=0N-1 h(x,y)e-2pj(mx/N+ny/N)
Inverting H – Overcome one problem
of the straightforward solution (cont.)
 Transpose H
• HT = WL*W-1 (B 6.6)
L* means the complex conjugate of L
 Example 6.19: Laplacian at a pixel position
• D2f(i, j) = f(i-1, j) + f(i, j-1) + f(i+1, j) + f(i, j +1) - 4f(i, j)
 Example 6.20: Identify L to estimate D2f(i, j)
 Example 6.21: Apply the Eq. of D2f(i, j)  L
 Example 6.22:
Constrained matrix inversion filter –
Overcome another problem