% Create a script that:
% ? reads in ElliottErwitt_Provence.jpg, converts it to grayscale (if necessary) and
%     converts it to double (right now with im2double -- this will change the range of values as well to [0 1])
% ? blurs the image (motion blur with pixel shifts of 21 and rotation of degrees 11; use fspecial to create the
%     psf of our ‘capturing system’, then use imfilter)
% ? create a noise layer separately, and add it to the blurred image (imnoise,
%     Gaussian-type with zero mean and 0.001 variance)
% ? calls your Wiener filter function in two different cases:
%     ? case1: with noise-to-signal power ratio (which is the variance of noise divided by the variance of the
%         input image; to compute the variance of the input please use var),
%     ? case2: with the ACF of your noise-layer (what you created separately) and with the ACF of the
%         input image (in this case use edgetaper on the noisy-and-blurred image to make smoother the
%         edge/frame-regions of your image; later we will reproduce this edgetaper function),
% ? pops up a figure with 3 subplots: original img, noisy img, reconstructed img. 

im = imread('ElliottErwitt_Provence.jpg');
if (size(im,3) ~= 1)
    im = rgb2gray(im);
end
im = im2double(im);

psf = fspecial('motion',21,11);
blurredIm = imfilter(im,psf,'replicate');

noise = imnoise(zeros(size(im)), 'gaussian',0,0.001);

blurryNoisy = blurredIm + noise;

figure(1)
subplot(1,3,1)
imshow(mat2gray(im*255))
subplot(1,3,2)
imshow(uint8((blurryNoisy)*255))
subplot(1,3,3)
imshow(uint8(wiener(blurryNoisy, psf, 0.001./var(im))*255))