Image Restoration

1
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

2
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

3
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 cxy d
• Fig 6.3
• Example 6.1

4
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) Sk1NSl1Nf(k, l)h(k, l, i, j)
• g H f

5
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

6
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

7
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

8
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)

9
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
  gt 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

10
Indirect solution Wiener filter

• Formal expression of the problem of IR
• To identify f(r) which minimizes e2 ? Ef(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)

11
Fourier transfer of the Wiener filter

• M(u, v) Fm(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
• Ev(r) 0
• ? Ef(r)v(r) Ef(r)Ev(r) 0

12
Fourier transfer of the Wiener filter (cont.)

• Create Sgf
• g(r) ? -??? -?? h(r - r?) f(r?)dr? v(r)
• Rgf (s) Eg(r)f(r - s) ? -??? -?? h(r - r?)
  Ef(r?)f(r - s)dr? Ef(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 HSff / SffH2 Svv
• M (1/H) H2 / H2 Svv/Sff

13
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

14
Fourier transfer of the Wiener filter (cont.)

• B6.1
• If m(r - r?) satisfies Ef(r) - ? -??? -?? m(r -
  r?) g(r?)dr?g(s) 0, then it minimizes e2 ?
  Ef(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) Fm(r) Sfg(u, v) / Sgg(u, v)

15
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

16
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

17
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)exp2pjk/M
  d(M-2)exp2pj2k/M d(1)exp2pj(M-1)k/M
• w(k) are the eigenvectors of D
• w(k) ? 1, exp2pjk/M, exp2pj2k/M, ,
  exp2pj(M-1)k/MT

18
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 exp2pj/N kn
• WN-1 (k, n) (N)-1/2 exp-2pj/N kn
• Kronecker product ?

19
Inverting H Overcome one problem of the
straightforward solution

• H is block circulant
• g H f
• g(i, j) Sk0N-1Sl0N-1h(k, l, i, j) f(k, l)
• For a shift invariant point spread function
• g(i, j) Sk0N-1Sl0N-1f(k, l) h(i-k, j-l)
• Diagonalize H
• H WLW-1 (B 6.5)
• WN (k, n) (N)-1/2 exp2pj/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) Sx0N-1Sy0N-1 h(x,y)e-2pj(mx/Nny/
  N)

20
Inverting H Overcome one problem of the
straightforward solution (cont.)

• Transpose H
• HT WLW-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(i1, 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

21
Constrained matrix inversion filter Overcome
another problem
22
(No Transcript)