ima = imread('flowers.tif');
ima = rgb2gray(ima);
ima = double(ima); 
ima = ima /max(max(ima)); 
sx = size(ima,1); 
sy = size(ima,2); 
buf = ima(1:2:sx,1:2:sy); 
buf1 = zeros(sx,sy); 
ima = buf1; 
ima(sx/2 - sx/4:sx/2+sx/4-1,sy/2-sy/4:sy/2+sy/4-1) = buf; 
clf;
figure(1); 
subplot(1,3,1);
imagesc(ima);
axis image;
colormap(gray); 
xlabel('original image');
disp('Now blur the image using a [7x7] mean kernel');
disp('Press any key to continue');
pause; % Now blur image
h = ones(7,7); 
ima1 = filter2(h,ima); 
subplot(1,3,2);
colormap(gray); 
imagesc(ima1);
axis image;
xlabel('Blurred image');
H = fft2(h,size(ima1,1),size(ima1,2)); 
disp('Now inverse filtering the inage'); 
% inverse filter 
Filter = (abs(H) > 0.0001) .* (1./H); 
ima_out = ifft2(fft2(ima1).*Filter); 
subplot(1,3,3);
imagesc(real(ima_out)); 
axis image;
xlabel('Deconvolved image');
drawnow;
disp('How sensitive is that to noise?');
disp('Press any key to continue');
pause;
ima2 = imnoise(uint8(ima1),'gaussian',0,0.001);
figure(2);
subplot(1,2,1);
imagesc(ima2);
axis image;
xlabel('original image with noise');
Filter = (abs(H) > 10) .* (1./H); 
%Filter = H./abs(H);
ima_out = ifft2(fft2(ima2).*Filter); 
subplot(1,2,2);
imagesc(real(ima_out)); 
drawnow;
axis image;
xlabel('Deconvolved noisy image');