1 / 22

Chapter 6

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

denali
Download Presentation

Chapter 6

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Chapter 6 Image Restoration Digital Image Processing Instructor: Dr. Cheng-Chien Liu Department of Earth Sciences National Cheng Kung University Last updated: 10October 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 + c4y^ = 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

  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) = Sk=1NSl=1Nf(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 w02M(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

  10. 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) andh(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) = 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) andv(r) are uncorrelated • E{v(r)} = 0 •  E{f(r)v(r)} = E{f(r)}E{v(r)} = 0

  12. 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 ]

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

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

  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 exp[2pj/Nkn] • WN-1(k, n) = (N)-1/2 exp[-2pj/Nkn] • Kronecker product 

  19. Inverting H – Overcome one problem of the straightforward solution • H is block circulant • g = H f • 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/Nkn] • WN-1(k, n) = (N)-1/2 exp[-2pj/Nkn] • L(k, i) = NH(kmodN, [k/N]) if i = k • L(k, i) = 0 if i  k • H(m,n) = (1/N) Sx=0N-1Sy=0N-1h(x,y)e-2pj(mx/N+ny/N)

  20. 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:

  21. Constrained matrix inversion filter – Overcome another problem

More Related