
 %  Program blur
 % Takes an initial image, and then blurs it using a "blurring"
 % operator, which is multiplication by a matrix H. 
 % The blurring operation is to take the average of the 4 nearest
 % neighbors. Then do it again. This acomplished by constucting
 % a blurring matrix H, and then multiplying by H^2.

 %   We then 
 % attempt do deblur the image by using the pseudo inverse of H^2.

 % We use a 16 by 16 mesh.  If this is too slow on your machine, try n = 10.
 n = 16;
 n = n; n2 = n^2;
 % construct the blurring matrix
 sdiag = ones(1,n2);
 s1 = zeros(1, n2-1);
 s2 = zeros(1, n2-n);
 s3 = s1;
 s4 = s2;
 for j = 1:n-2
   for i = j*n+2:(j+1)*n -1
     sdiag(i) = 0;
     s1(i) = .25;
     s2(i) = .25;
     s3(i-1) = .25;
     s4(i-n) = .25;
   end
 end 

 D = diag(sdiag); 
 Eupper = diag(s1,1) + diag(s2, n);
 Elower = diag(s3,-1) + diag(s4, -n);
 % this is the blurring matrix
 H = D +Eupper + Elower;
 %Square of the blurring matrix
 H2 = H*H;

 % level sets of this function form the image
 f = inline('(x.^2 + y.^2 < 3)', 'x', 'y');
 x = linspace(-2,2,n);
 [X,Y] = meshgrid(x);
 % original image
 Z = f(X,Y);
 Z(7:10, 7:10) = zeros(4,4);
 % put all of Z into a column vector z
 z = Z(1,:)';
 for j = 2:n
    z = [z; Z(j,:)'];
 end 
 % apply the bluring operator H2. H takes leaves the boundary points alone,
 % and at each interior point, it takes the average of the 4 nearest 
 % neighbors.


 zz = H2*z;
 %reassamble the blurred image
 ZZ = zz(1:n)';
 for j = 1:n-1
   ZZ = [ZZ; zz(n*j+1:(j+1)*n)'];
 end
 % ZZ is the blurred image.
 subplot(1,2,1)
   image(64*Z) 
   colormap(gray)
   title('Original image' ) 
  
   disp('Hit return to see blurred image ')
   pause
 subplot(1,2,2)
   image(64*ZZ)
   title('Blurred image ')
 r = rank(H2);
 sprintf('Dimension of blurring matrix H2 is %d x %d. Rank of H2 is %d.', n2,n2, r)
 % the peudo inverse of H
 % deblurred vector 
 zzz = pinv(H2)*zz;
  % construction of the deblurred image from the deblurred vector
 ZZZ = zzz(1:n)';
 for j = 1:n-1
    ZZZ = [ZZZ; zzz(n*j+1:(j+1)*n)'];
 end
 disp('Hit return to see deblurred image  ')
 pause 
 subplot(1,2,2)
   image(64*ZZZ)
   title('Deblurred image using pseudo inverse ')
 
    




 



   
 
