%%
Nx = 256;
sig = 6;
nlevel = 2e-8;
itol = 1e-12; % to avoid divide by zero.
MultiplicativeNoise = false; % true;
SheppLogan = true;
if(SheppLogan)
    im = phantom(Nx);
else
    pt = imread('../StandardTestImages/Cameraman256.png','png'); pt = double(imresize(pt,[Nx Nx]));
    
    im = pt;
end
cmap = 'bone'; % 'gray','bone';'hote' etc...

H = fspecial('gaussian',Nx,sig);
figure(1);clf; hold on;colormap(cmap);
subplot(2,6,1); imagesc(im);title('true image');
subplot(2,6,2); imagesc(H);title('blurring kernel');

fim = fft2(im); fH = fft2(H);
subplot(2,6,7); imagesc(fftshift(log(abs(fim)))); 
subplot(2,6,8); imagesc(fftshift(log(abs(fH)))); 


%%
fimH = fim.*fH;
imH = ifftshift(ifft2(fimH));
figure(1);colormap(cmap);
subplot(2,6,3); imagesc(abs(imH));title('blurred image');
subplot(2,6,9); imagesc(fftshift(log(abs(fimH)))); 

%% Add noise
if MultiplicativeNoise
imHn = imH.*(1+nlevel*randn(size(im)));
else
imHn = imH +nlevel*max(abs(imH(:)))*randn(size(im))/(Nx*Nx);
end

%%
fimH = fft2(imH);
fimHn = fft2(imHn);

figure(1);colormap(cmap);
subplot(2,6,3); imagesc((imHn));title(['blurred image + ',num2str(100*nlevel),'%noise']);
subplot(2,6,9); imagesc(fftshift(log(abs(fimHn)))); 


%%

fimR = (fimH./(fH+itol));
imR = ifft2(fimR);
imR = ifftshift(imR);
figure(1);colormap(cmap);
subplot(2,6,4); imagesc(abs(imR));title('deblurred from noiseless data')
subplot(2,6,10); imagesc(fftshift(log(abs(fimR))));


%%
fimRn = (fimHn./(fH+itol));
imRn = ifft2(fimRn);
imRn = ifftshift(imRn);
figure(1);;colormap(cmap);
subplot(2,6,5); imagesc(abs(imRn));title('deblurred from noisy data')
subplot(2,6,11); imagesc(fftshift(log(abs(fimRn))));


%%
alpha =  1e-4*nlevel*max(abs(fimH(:)))/(Nx*Nx);
fimRnalpha = fimHn.* (fH./(fH.^2+alpha));
imRnalpha = ifft2(fimRnalpha);
imRnalpha = ifftshift(imRnalpha);
figure(1);
subplot(2,6,6); imagesc(abs(imRnalpha));title(['deblurred from noisy data \alpha=',num2str(alpha)])
subplot(2,6,12); imagesc(fftshift(log(abs(fimRnalpha))));