% EE362 Final Project
% Lucy-Richardson Restoration

clear

% Load Image
image = imread('bibby.jpg');
image = rgb2gray(image); %Convert to Grayscale
image = im2double(image); %Convert to Double Representation

figure
imshow(image)
title('Original')

% Blur Parameters
len = 20;
angle = 45;

% # Iterations Resolution
K = 10;

% Number Iterations in Lucy-Richardson Algorithm * K
N = 1;
% number_it = K*N

% Simulated Blur Function
blur = fspecial('motion',len,angle);

% Blur Image
image_blur = conv2(blur, image);

% Start and End Indices for Restored Image %%%%%%%%%%%%%%%%
% (Restored Image is Larger than Original)

% Get Size of Original Image
im_size_orig = size(image);
% Get Size of Blurred Image
im_size_restore = size(image_blur);
% Image Start and End Indices
rowstart = (im_size_restore(1) - im_size_orig(1))/2 + 1;
rowend = im_size_restore(1) - rowstart + 1;
colstart = (im_size_restore(2) - im_size_orig(2))/2 + 1;
colend = im_size_restore(2) - colstart + 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Crop Blurred Image %%%%%
image_blur_crop = image_blur([rowstart:rowend],[colstart:colend]);

figure
imshow(image_blur_crop)
title('Blurred')
%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate Error from Blurring %%%%%%%%%%

err(1) = mean((image(:)-image_blur_crop(:)).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate Variance of Blurred - Original Image

Var_diff = var(image_blur_crop(:) - image(:));

Npower_orig = 10*log10( Var_diff );

% Calculate SNR of Blurred Image

Var_image = var(image(:));

Spower_orig = 10*log10( Var_image );

SNR_blur = 10*log10( Var_image / Var_diff );

% Apply Lucy Restoration for i*K Iterations

for i= 1:N

    imlucy = LucyMartin(image_blur, len, angle, i*K);
  
    % Crop Image
    imlucy = imlucy([rowstart:rowend],[colstart:colend]);
    
    % SNR Calculations
    Var_restore = var(imlucy(:) - image(:));
    deltaSNR(i) = 10*log10( Var_diff / Var_restore );
    
    % MSE
    err(i+1) = mean((image(:)-imlucy(:)).^2);

    figure
    imshow(abs(imlucy), [])
    title(['Restored, ' num2str(i*K) ' Iterations'])

end

% Plot

figure
plot(K*[1:N],deltaSNR)
title('Delta SNR vs. # Iterations')
ylabel('Delta SNR (dB)')
xlabel('# Iterations')

figure
plot(K*[0:N],err)
title('Mean Squared Error vs # Iterations')
ylabel('MSE')
xlabel('# Iterations')