%% 1D Convolution, Smoothing Reminder
u = -5:.1:5;
w = sin(u*pi)+0.5*randn(size(u));
plot(u,w)
%% Convolve with [1,1,1] kernel
wsm = conv(w,[1,1,1]);
whos wsm
%% Overplot
plot(u,w)
hold on; plot(u,wsm(2:102)/3,'r','linewidth',2); hold off
%% 2D Convolution (Reminder)
x = -5:0.1:5;
y = -5:0.1:5;
% gaussfit2d
ztest = params(1)*exp(-(((X-params(2))/params(4)).^2 + ((Y-params(3))/params(5)).^2));
k = ones(3);
zsm = conv2(ztest-Z,k)/9.;
colormap('default');
imagesc(x,y,zsm(2:102,2:102))
%% Other convolution kernels 
%% horizontal line
k = [1,1,1,1,1,1];
zsm = conv2(ztest-Z,k)/9.;
imagesc(x,y,zsm)
%% vertical line
k = [1,1,1,1,1,1]';
zsm = conv2(ztest-Z,k)/9.;
imagesc(x,y,zsm)
%% Sinc function
x = -5:0.1:5; y = -5:0.1:5;
[X, Y] = meshgrid(x,y);
Z = sinc(sqrt(X.^2 + Y.^2));
imagesc(x,y,Z);
%% Electric field vs. Power
imagesc(x,y,Z.^2);
%% Create some artificial stars
stars = randn(size(X))*0.1;
stars(50,50) = 1; stars(20,37) = 4; stars(87,74) = 2; stars(45,24) = 0.5;
colormap(gray(256));
imagesc(stars)
%% Convolve with point spread function
Z = sinc(sqrt(X.^2 + Y.^2)*5).^2;   % the *5 makes it smaller/sharper
imagesc(Z)
%% Convolved stars
imagesc(conv2(stars,Z))
%% Adjust axis to a square image
fuzzy = conv2(stars,Z); colormap(gray(256));
imagesc(stars); axis square
%% Overplot fuzzy stars
imagesc(fuzzy(51:150,51:150)); axis square
%% Blink compare
for i = 1:5; 
    imagesc(stars); axis square; 
    pause(0.5); 
    imagesc(fuzzy(51:150,51:150)); axis square; 
    pause(0.5); 
end
%% Read Image
[img map] = imread('fruit.gif');
%% Convolve with point spread function and clip to original size
fuzzy = conv2(single(img),Z)/sum(sum(Z));
fuzzy = fuzzy(51:515,51:750);
%% Display original image
image(img)    % Original image—observe the sharpness
%% Display image convolved with point spread function
image(fuzzy)   % fuzzy image
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(fuzzy) 
    pause(0.5); 
end
%% Convolve with larger point spread function
Z = sinc(sqrt(X.^2 + Y.^2)*2).^2;   % the *2 makes it larger/fuzzier
imagesc(Z)
%% Convolve with point spread function and clip to original size
fuzzy = conv2(double(img),Z)/sum(sum(Z));
%fuzzy = fuzzy(51:515,51:750);
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(fuzzy(51:515,51:750)) 
    pause(0.5); 
end
%% Deconvolution (regular deconvolution)
dconr = deconvreg(fuzzy,Z);
image(dconr)
%% Normalize by area of point spread function
image(dconr*sum(sum(Z)))
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(dconr(51:515,51:750)*sum(sum(Z))) 
    pause(0.5); 
end
%% Deconvolution (Lucy deconvolution)
dconl = deconvlucy(fuzzy,Z)*sum(sum(Z));
image(dconl)
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(dconl(51:515,51:750)) 
    pause(0.5); 
end
%% Image difference
imagesc(dconl(51:515,51:750)- single(img))
%% Equivalent Sinc and Gaussian functions
Zsinc = sinc(sqrt(X.^2 + Y.^2)*5).^2;
Zgaus = exp(-(X.^2 + Y.^2)*(5*1.913)^2);
%% Blink compare
for i = 1:5; 
    imagesc(Zsinc,[0,0.3]); axis square; 
    pause(0.5); 
    imagesc(Zgaus,[0,0.3]); axis square; 
    pause(0.5); 
end
%% Convolve with Gaussian (i.e. instrument + atmosphere blurring)
fuzzy = conv2(double(img),Zgaus)/sum(sum(Zgaus));
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(fuzzy(51:515,51:750)) 
    pause(0.5); 
end
%% Deconvolution again (Lucy deconvolution), but now with only instrumental
%  Sinc function
dconl = deconvlucy(fuzzy,Zsinc)*sum(sum(Zsinc));
image(dconl)
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(dconl(51:515,51:750)) 
    pause(0.5); 
end
%% Even larger Gaussian function (i.e. more atmospheric blurring)
Zgaus = exp(-(X.^2 + Y.^2)*(2*1.913)^2);
%  Convolve with Gaussian (i.e. instrument + atmosphere blurring)
fuzzy = conv2(double(img),Zgaus)/sum(sum(Zgaus));
%  Deconvolution again (Lucy deconvolution), but now with only instrumental
%  Sinc function
dconl = deconvlucy(fuzzy,Zsinc)*sum(sum(Zsinc));
image(dconl)
%% Blink compare
for i = 1:5; 
    image(img) 
    pause(0.5); 
    image(dconl(51:515,51:750)) 
    pause(0.5); 
end





