By Mary Hudachek-Buswell Overview Atmospheric Turbulence Blur k ( x y ) H ( x, y) e 2 2 g ( x, y) H [ f ( x, y)] n( x, y ) G A F E 2-D Separable Gaussian Blur k ( x y ) H ( x, y) e e k ( x y ) 2 Ae 2 k x 2 e 2 k x 2 2 e k y and B e 2 k y This point spread function represents atmospheric turbulence blur. 2 Toeplitz Matrix Blur Toeplitz matrix structure is a matrix with constants along the descending diagonal. These type matrices are spatially invariant and have zero boundary conditions when combined together. Illustration of Spatially Invariant BLur Original image Invariant Blur Variant Blur Illustration of Zero Boundary Conditions 0 0 0 0 X 0 0 0 0 Surrounds the object with a black boundary, or zeros on the outside borders of the image Minimize the 2-Norm of the Residual A X B C R 2 A X B C 2 A is the column blur, B is the row blur, X is the restored image, and C is distorted image. A & B are overdetermined and positive definite. Overdetermined systems have more equations than unknowns and positive definite matrices zT A zhave 0 a unique quality where given any nonzero vector, z, the product Least Squares of an Overdetermined System A B C I X I 0 T A A B T X A I I I C I 0 A B A I X AT C I I B T 2 T A A I X A I C T Least Squares cont. T T B B A C T 2 T A A I X I I 0 B A T I A I X B B I A C B 2 T 2 T T X A A I A C B B B I T 2 1 T T T 2 1 Necessary Properties The system of recurrence formulas that generates a unique sequence of iterates with X n Kn such that n 1 K n X , A X B, A X B ,..., A X B 2 2 With the property that at step n, the norm is minimized n 1 Conjugate Gradient Algorithm X 0 0, R0 C , P0 C For n 1, 2,3,... RnT1 Rn1 n T Pn1 A Pn1 B Initial values for approximation, residual and gradient matrices Scalar to minimize the norm and determine step length X n X n1 n Pn1 Approximate restoration matrix Rn Rn1 n A Pn1 Residual matrix RnT Rn n T Rn1 Rn1 Scalar to measure residual ratio and update gradient Pn Rn n Pn1 Update conjugate gradient search direction Conjugate Gradient MATLAB C = g; ro = g; p = ro; x = zeros(m, n); for i = 1:maxiter alpha = trace (ro'*ro)/trace(p'*(A*p*B)); x = x + alpha * p; rn = ro - alpha * A * p * B; beta = trace (rn'*rn)/trace(ro'*ro); p = rn + beta * p; ro = rn; End The trace function computes the sum of the diagonal elements of the matrix Comparison between CG & Inverse Methods Original Blurred Noisy Image Deblur Inverse Method Conjugate Iterative Deblur Method RMSE = 0.20442 CG iter = 9 RMSE = 0.060552 Visual Comparison of Iterations Original Iter 6 Iter 7 Blurred Noisy Image RMSE = 0.051591 RMSE = 0.051865 Iter 8 Iter 9 Iter 10 Iter 11 RMSE = 0.053561 RMSE = 0.056487 RMSE = 0.066383 RMSE = 0.082311 Iter 12 RMSE = 0.10071 Iter 13 RMSE = 0.12304 Iter 14 RMSE = 0.14584 Iter 15 RMSE = 0.1758 Concluding Thoughts Direct Methods to find the restoration of a blurred image involve inverting matrices which costs O(n3) Conjugate gradient iterations require transposing matrices, then multiplying by scalars and matrices which costs O(n2), considerably less than direct methods The CG computations are faster and can provide better visual results.