260 likes | 386 Views
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
E N D
Chapter 6 Image Restoration Digital Image Processing Instructor: Dr. Cheng-Chien Liu Department of Earth Sciences National Cheng Kung University Last updated: 10October 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 + 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
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
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 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
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)
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
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/Nkn] • WN-1(k, n) = (N)-1/2 exp[-2pj/Nkn] • Kronecker product
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)
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