clear all;
close all;
%%%%%
% CTF parameters
%%%%%
Cs     = 2.;%(mm)

% 100 keV, 200 keV, and 300 keV is 3.70 pm, 2.51 pm and 1.96 pm, 
%lambda = 1.97;%(pm)
lambda = 0.00197;%(nm)
dZ1    = 2300;%(nm)
dZ2    = 500;%(nm)

w = 0.1;

f=linspace(-1.43,1.43,456);%-0.04:0.001:0.04;
[X,Y] = meshgrid(f,f);
R = sqrt(Y.^2+X.^2);
%1D ctf
figure(1)
CTF = ctf(R,Cs,lambda,dZ1,w);
imagesc(f,f,CTF);
colormap gray;
figure(2)
CTF = ctf(R,Cs,lambda,dZ2,w);
imagesc(f,f,CTF);
colormap gray;
%return

%%%%%%%%%%%%%%
%process einstein image
%%%%%%%%%%

Epsilon = 0.01
dev =10%128 is a good value
%read phantom
im = double(imread('einstein.png'));
%figure(1)
subplot(2,3,1); imagesc( im); title('Original Phantom');
axis image; axis off; colormap gray(256);

%%rand_image = dev * randn(size(CTF));
%rand_image = 0 * randn(size(CTF));
%%im =  im + rand_image;
%%subplot(2,3,2); imagesc(im); title('Noisy phantom');
%%axis image; axis off; colormap gray(256);

% convolve with ctf
IM = fft2(im);
CON = IM.*CTF;
con = ifft2(CON);
subplot(2,3,2); imagesc(real(con)); title('Phantom convolved with CTF');
axis image; axis off; colormap gray(256);

% add noise to CTF 
rand_image = dev * randn(size(CTF));
im =  con + rand_image;
IM = fft2(im);
CON = IM;
subplot(2,3,3); imagesc(real(im)); title('Noisy phantom convolved with CTF');
axis image; axis off; colormap gray(256);

%phase flip
RESTORE_1 = CON.*sign(CTF);
%figure(4)
%imagesc(real(ifft2(RESTORE_1)));
subplot(2,3,4); imagesc(real(ifft2(RESTORE_1))); title('Partialy recover using phase flip');
axis image; axis off; colormap gray(256);

%Wiener filtering_2
RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon));
subplot(2,3,5); imagesc(real(ifft2(RESTORE_2))); title(sprintf('Partialy recover using E=%f',Epsilon));
axis image; axis off; colormap gray(256);

%Wiener filtering_3
RESTORE_3 = CON.*(CTF./(CTF.*CTF + Epsilon*10));
subplot(2,3,6); imagesc(real(ifft2(RESTORE_3))); title(sprintf('Partialy recover using E=%f',Epsilon*10));
axis image; axis off; colormap gray(256);

% %Wiener filtering_2
% RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon*10));
% figure
% imagesc(real(ifft2(RESTORE_2)));
% title(sprintf('Partialy recover using E=%f',Epsilon*10))
% colormap gray
% %Wiener filtering_3
% RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon*100));
% figure
% imagesc(real(ifft2(RESTORE_2)));
% title(sprintf('Partialy recover using E=%f',Epsilon*100))
% colormap gray
% %Wiener filtering_4
% RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon*1000));
% figure
% imagesc(real(ifft2(RESTORE_2)));
% title(sprintf('Partialy recover using E=%f',Epsilon*1000.))
% colormap gray
% %Wiener filtering_5
% RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon*10000));
% figure
% imagesc(real(ifft2(RESTORE_2)));
% title(sprintf('Partialy recover using E=%f',Epsilon*10000.))
% colormap gray
% %Wiener filtering_6
% RESTORE_2 = CON.*(CTF./(CTF.*CTF + Epsilon*100000));
% figure
% imagesc(real(ifft2(RESTORE_2)));
% title(sprintf('Partialy recover using E=%f',Epsilon*100000.))
% colormap gray
return








