% Optical Imaging and Spectroscopy
%
% David J. Brady
% Duke University
% www.opticalimaging.org
%
%
%
% goofing off with lucy richardson
%
figure(2);
set(gcf,'color','white');
colormap 'gray';
%mwir=double(imread('MWIR087c.png'));
%pip=size(mwir);
%mwir=interp2(mwir,linspace(1,pip(2),256),linspace(1,pip(1),256)');
%imagesc(mwir);
%I = double(imread('cameraman.tif'));
%range=linspace(20,20,256);
%[x,y]=meshgrid(range,range);
%I=(1+cos(((x+10).^2+(y+10).^2)/7))/2;
I=round(0.501*rand(256));
%I=mwir;
I=I-min(min(I));
I=I/max(max(I));
I=[zeros(128,256);I;zeros(128,256)];
I=[zeros(512,128),I,zeros(512,128)];
subplot(2,2,1);imagesc(I(129:385,129:385));title('Original Image');
PSF=round(.51*rand(32));
PSF=PSF/sum(sum(PSF));
%PSF = fspecial('gaussian',5,5);
%PSF=[1 0;0 0];
Blurred = imfilter(I,PSF,'symmetric','conv');
V =.00001;
BlurredNoisy = imnoise(Blurred,'gaussian',0,V);
ds=kron(eye(128),[1 1 1 1]);
BlurredNoisyDS=ds*BlurredNoisy*ds';
subplot(2,2,2);imagesc(BlurredNoisyDS);title('Blurred Downsampled Image');
luc1 = deconvlucy(BlurredNoisyDS,PSF,25,0,ones(128),0,4);
subplot(2,2,3); imagesc(luc1(129:385,129:385));
title('Restored Image');
luc1=luc1-min(min(luc1));
luc1=luc1/max(max(luc1));
subplot(2,2,4);imagesc(PSF);title('PSF');
%deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)
figure(1);
set(gcf,'color','white');
colormap 'gray';
subplot(2,2,1);imagesc(I);title('Original Image');
%PSF=round(.504*rand(32));
 %PSF=PSF/sum(sum(PSF));
%PSF = fspecial('gaussian',5,5);
PSF=[1 0;0 0];
Blurred = imfilter(I,PSF,'symmetric','conv');
V =.00001;
BlurredNoisy = imnoise(Blurred,'gaussian',0,V);
ds=kron(eye(128),[1 1 1 1]);
BlurredNoisyDS=ds*BlurredNoisy*ds';
subplot(2,2,2);imagesc(BlurredNoisyDS);title('Blurred, Noisy downsampled Image');
luc1 = deconvlucy(BlurredNoisyDS,PSF,30,0,ones(128),0,4);
subplot(2,2,3); imagesc(luc1);
title('Restored Image');
luc1=luc1-min(min(luc1));
luc1=luc1/max(max(luc1));
subplot(2,2,4);imagesc(PSF);title('PSF');
% %deconvlucy(I,PSF,NUMIT,DAMPAR,WEIGHT,READOUT,SUBSMPL)