% Matlab program for image deblurring problem.
% Rich Kozick, Spring 1997.

% Read in the blurred and noisy image, and display it.

load blurred
figure(1)
imshow( mat2gray(B), 256 );
title('Blurred and noisy image')

% Now you need to figure out what process has blurred this image.
% A hint is to look at the upper left corner of the image.  You should
% see a faint streak.  From previous photographs, it is known that this
% streak corresponds to a blurred version of a bright star.  If the image
% were not blurred, the star would cover a 5 x 5 square of pixels.

% Define H as a K x K matrix that is the blurring function.  In other
% words, the matrix B is a 2-D convolution between the true image (which
% you want to restore) and H.  Matrix B also contains additive noise.

% Define K and H on the following lines.  (Remove the '%' signs.)

% K = ;
% H = ;

[M, N] = size(B);

fprintf('Computing FFTs of B and H ... ')
Bf = fft2(B, M, N);
Hf = fft2(H, M, N);	% Must zero-pad H before fft
fprintf('Done\n')

% Define the threshold to use in placing zeros in the approximation to
% the inverse of Hf.  If you choose a small threshold, then you will
% do better deblurring, but the noise will be amplified more.  If you
% choose a larger threshold, then less deblurring is done, but the
% noise is also kept under control.  You must find the right tradeoff.
% If done correctly, you can get a very good final image.

% The threshold is adapted to the range of values in Hf by comparing all
% values with the MAXIMUM element in Hf.

% Define your threshold below.  (Remove the '%' sign.)

% thresh = ;

% Now the rest of the program does the processing for you.
% It creates the inverse to Hf, with entries set to zero according to
% your threshold.  Then it multiplies this inverse times the FFT of
% the blurred, noisy image, and finally does an inverse FFT to get
% the final image.

fprintf('Computing Hfinv ... ')
maxHf = max(max(abs(Hf)));
[zr,zc] = find(abs(Hf/maxHf) < thresh); % Find elements below threshold.
for k = 1:size(zr,1)
  Hf(zr(k),zc(k)) = 1;	% Temporarily set equal to 1 to avoid division by 0.
end
Hfinv = 1 ./ Hf;	% Take reciprocals.
for k = 1:size(zr,1)
  Hfinv(zr(k),zc(k)) = 0; % Now set components to 0 in inverse
end
fprintf('Done\n')

fprintf('Computing IFFT for deblurred image ... ')
Ihat = ifft2( Hfinv .* Bf );
fprintf('Done\n')
figure(3)
imshow( mat2gray( abs( Ihat(1:(M-K+1), 1:(N-K+1)) ) ), 256);
title('Deblurred image')
