% Project 2
% Andy Doran
%
% function baran(name,xdim,M,No)
%
% name  = 'input image' (as in 'lenna.256')
% xdim  = x dimension of input image (usually 256 or 512)
% M     = Number of Images to average over.
% No    = Variance of the noise in the filter
%
% This function takes an input image and runs it through a LSI filter
% M times, generating M degraded images.  We assume we know the average
% transfer function of the filter H.  The degraded images are taken to
% the Fourier domain and geometrically averaged, then multiplied by the
% average H.  Our LSI filter H is generated differently each time with
% zero mean Gaussian noise added each time to the filter components.
% This method comes Elec539 class notes).
% 
% The average blurring filter is calculated as h = zeros(4);

function baran(name,xdim,M,No)

% Load and Plot image
pict = freadbin(name,xdim,xdim);

figure(1)
clf
colormap(gray);
imagesc(pict);
txt = [num2str(name) ' before degradation'];
title(txt)
axis square
axis off

% Create LSI degradation model
blur = 4;
h = zeros(blur);

% Initialize degraded frequency spectra
DEG = ones(xdim);

% Initialize average logarithm transfer function
HM = zeros(xdim);
hm = zeros(4);

% Create Gaussian noise, mean = 0, variance comes from input (No)
% To insure ergodicity, make sure each noise filter index has zero
% mean over all the input images.
% To insure zero mean noise, subtract out the mean each time.

h1 = sqrt(No)*randn(1,M);
h1 = h1 - mean(h1);
h2 = sqrt(No)*randn(1,M);
h2 = h2 - mean(h2);
h3 = sqrt(No)*randn(1,M);
h3 = h3 - mean(h3);
h4 = sqrt(No)*randn(1,M);
h4 = h4 - mean(h4);
h5 = sqrt(No)*randn(1,M);
h5 = h5 - mean(h5);
h6 = sqrt(No)*randn(1,M);
h6 = h6 - mean(h6);
h7 = sqrt(No)*randn(1,M);
h7 = h7 - mean(h7);
h8 = sqrt(No)*randn(1,M);
h8 = h8 - mean(h8);
h9 = sqrt(No)*randn(1,M);
h9 = h9 - mean(h9);
h10 = sqrt(No)*randn(1,M);
h10 = h10 - mean(h10);
h11 = sqrt(No)*randn(1,M);
h11 = h11 - mean(h11);
h12 = sqrt(No)*randn(1,M);
h12 = h12 - mean(h12);
h13 = sqrt(No)*randn(1,M);
h13 = h13 - mean(h13);
h14 = sqrt(No)*randn(1,M);
h14 = h14 - mean(h14);
h15 = sqrt(No)*randn(1,M);
h15 = h15 - mean(h15);
h16 = sqrt(No)*randn(1,M);
h16 = h16 - mean(h16);

% Create degraded images and get geometric mean of Fourier Magnitudes
for i = 1:M
  hd = h + [h1(i) h2(i) h3(i) h4(i);h5(i) h6(i) h7(i) h8(i);h9(i) h10(i) h11(i) h12(i);h13(i) h14(i) h15(i) h16(i)];
  HM = HM + (abs(fft2(hd,xdim,xdim)));
  hm = hm + hd;
  % Run image through LSI Filter with additive noise set to 0
  dpict = distimag(pict,hd,0);
  DPICT = fft2(dpict);
  DEG = DEG.*abs(DPICT.^(1/M));    % Get geometric mean of degraded spectrum
end

% Plot restored image spectra
figure(1)
clf
imagesc(fftshift(abs(DEG)));
colormap(gray);
txt = ['Frequency Spectra of Restored Image'];
title(txt)
axis square
axis off

% Calculate Transfer function of average LSI filter
H = fft2(h,xdim,xdim);

% Restore Image
scale = exp((-H)/M);
index = find(scale == 1);
length(index)
RPICT = scale.*DEG;
rpict = ifft2(DEG);

% Plot restored image
figure(3)
clf
imagesc(abs(rpict));
colormap(gray);
txt = [num2str(name) ' restored using Homomorphic Restoration'];
title(txt)
axis square
axis off

% Calculate MSE of restored image
error = abs(rpict) - pict;
sqerr = sum(sum(error.^2));
RMSE = sqerr/(xdim^2)

% Calculate a degraded image to test
noise = sqrt(No)*randn(blur);
hd = h + noise
dpict = distimag(pict,hd,0);

% Plot abs(fft(h)) to see type of filter we have
H = fft2(hd,xdim,xdim);

figure(4)
clf
colormap(gray);
imagesc(fftshift(abs(H)));
txt = ['Magnitude of blurring filter'];
title(txt)
axis square
axis off

% Plot degraded image
figure(2)
clf
colormap(gray);
imagesc(dpict);
txt = [num2str(name) ' degraded with 4x4 averaging Filter with Noise'];
title(txt)
axis square
axis off

% Plot magnitude of original image
PICT = fft2(pict);
figure(5)
clf
imagesc(fftshift(abs(PICT)));
colormap(gray);
txt = ['Frequency Spectra of Original Image'];
title(txt)
axis square
axis off
% Calculate MSE of degraded image
error = dpict - pict;
sqerr = sum(sum(error.^2));
DMSE = sqerr/(xdim^2)
