Conjugate Gradient Iterative Method for Deblurring Images

```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
Ae
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
X 0  0, R0  C , P0  C
For n  1, 2,3,...
RnT1  Rn1
n  T
Pn1   A  Pn1  B 
Initial values for approximation,
Scalar to minimize the norm and
determine step length
X n  X n1  n  Pn1
Approximate restoration matrix
Rn  Rn1  n  A  Pn1
Residual matrix
RnT  Rn
n  T
Rn1  Rn1
Scalar to measure residual ratio
Pn  Rn  n  Pn1
direction




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)