%%%%%%%%%%%% ENEE 631 HOMEWORK #3 %%%%%%%%%%%%%%%%%
% Image restoration using:
% 1. Pseudo-inverse filter
% 2. Wiener filter
%
% Programmed by Daniel Garcia-Romero 3-5-2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load images
 I_motion=imread('.\motion_blur.jpg');
 I_blur=imread('.\out_of_focus.jpg');

% Allocate memory for resulting images
I_motion_inverse=zeros(768,1024,3);
I_blur_inverse=zeros(768,1024,3);

I_motion_wiener=zeros(768,1024,3);
I_blur_wiener=zeros(768,1024,3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process image with out of focus blur: I_blur
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Use a blur model with a psf consisting of a disk with r=4. The value of has
% been experimentally obtained (subjective visual quality) as a trade off between
% a sharp image and ringing effect. 
psf=fspecial('disk',4);


[I_blur_inverse] = my_pseudo_inverse_filter(I_blur, psf, 0.5, 50); % Best experimental parameters for I_blur

[I_blur_wiener] = my_wiener_filter(I_blur, psf, 0); %Noise to Signal ratio set to 0 since the image has no noise

figure;
imshow(I_blur,[]); title('Original out of focus image');

figure;
imshow(I_blur_inverse,[]); title('Pseudo-inverse filter processed');
%imwrite(I_blur_inverse,'.\temp\I_blur_inverse.jpg');

figure;
imshow(I_blur_wiener,[]); title('Wiener filter processed');
%imwrite(I_blur_wiener,'.\temp\I_blur_wiener.jpg');


% Since the ringing is most prominent in the borders, we are going to
% create mask to combine the border of the restored image with the 
% degraded image. A frame of 20 pixels is good enough for a psf='disk' of r=4;
mask=ones(768,1024,3);
%Rows
mask(1:20,:,:)=0;
mask(768-20:768,:,:)=0;
%Columns
mask(:,1:20,:)=0;
mask(:,1024-20:1024,:)=0;

I_mix_inverse=mask.*double(I_blur_inverse)+(1-mask).*double(I_blur);
I_mix_inverse=uint8(I_mix_inverse);
figure;
imshow(I_mix_inverse,[]); title('Mixture');
%imwrite(I_mix_inverse,'.\temp\I_blur_mix.jpg');





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Process image with motion blur: I_motion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The motion blur in the image is modelled as a linear
% movement with length 23 pixels and angle =0. This parameters can be 
% easily infered base on the length of the black borders on both sides of
% the image.

psf = fspecial('motion',23,0);

% In order to reduce the ringing effects introduced by using a DFT based
% filer (due to border effect) we are going to smooth the lateral edges of
% the image using a gaussian filter with sigma =40 and size=[1 80]

I_tap_motion= edgetaper(I_motion,fspecial('gaussian',[1 80],40));

 [I_motion_inverse] = my_pseudo_inverse_filter(I_tap_motion, psf,0.5,200); % Best experimental parameters for I_blur

 [I_motion_wiener] = my_wiener_filter(I_tap_motion, psf, 0); %Noise to Signal ratio set to 0 since the image has no noise


figure;
imshow(I_tap_motion,[]); title('Original motion blur image with smoothed lateral edges');
%imwrite(I_tap_motion,'.\temp\I_motion_smoothed_edges.jpg');

figure;
imshow(I_motion_inverse,[]); title('Pseudo-inverse filter processed');
%imwrite(I_motion_inverse,'.\temp\I_motion_inverse.jpg');

figure;
imshow(I_motion_wiener,[]); title('Wiener filter processed');
%imwrite(I_motion_wiener,'.\temp\I_motion_wiener.jpg');