%% coh_image: Example for Coherent Imaging
close all; clear;clc;
circ=@(ri) (abs(ri)<=1.) ;

A=imread('LSsnap_cell','tif');  %read image file 
[M,N]=size(A);       %get image sample size 
A=flipud(A);           %reverse row order 
Ig=single(A);          %integer to floating 
Ig=Ig/max(max(Ig));    %normalize image 

%%
ug=sqrt(Ig);         %Prepare image field 
L=3.5e-3;              %image plane side length (m) 
du=L/M;                %sampling interval
u=-L/2:du:L/2-du; v=u; 

%%  
figure(1)              %check ideal image 
imagesc(u,v,Ig); 
colormap('gray'); 
xlabel('u (m)');
ylabel('v (m)'); 
axis square 
axis xy 
title('Object');

%% Define the imaging system parameters and 
% generate the coherent transfer function (CTF)
lambda=0.6*10^-6;    %wavelength 
xpd=1.0e-3;            % exit pupil diameter 
zxp=120e-3;              %exit pupil distance 
fc=xpd/(lambda*zxp); % cutoff frequency 
%
fu=-1/(2*du):1/L:1/(2*du)-(1/L); %freq coords 
fv=fu; 
[Fu,Fv]=meshgrid(fu,fv); 
H=circ(sqrt(Fu.^2+Fv.^2)/fc); 

%%
figure(2)              %check H 
surf(fu,fv,H.*.99) 
camlight left; lighting phong 
colormap('gray') 
shading interp 
ylabel('fu (cyc/m)'); 
xlabel('fv (cyc/m)'); 
title('Coherent transfer function');

%% Diffraction-Limited Imaging Simulation
H=fftshift(H); 
Gg=fft2(fftshift(ug)); 
Gi=Gg.*H; 
ui=ifftshift(ifft2(Gi)); 
Ii=(abs(ui)).^2; 
%
figure(3)              %image result 
imagesc(u,v,nthroot(Ii,2)); 
colormap('gray'); xlabel('u (m)'); ylabel('v (m)'); 
axis square 
axis xy 
title('Imaging by coherent optical system');
%
figure(4)                  %horizontal image slice 
vvalue=-0.8e-4;   %select row (y value) 
vindex=round(vvalue/du+(M/2+1)); %convert row index 
plot(u,Ii(vindex,:),'-',u,Ig(vindex,:),':'); 
xlabel('u (m)'); ylabel('Irradiance'); 
title('object: dot , image: line');





