objectIntensity = double(imread('cameraman.tif'));

%1) Coherent imaging

figure(1);

subplot(221); imshow(objectIntensity, []);

title('Input  object');

objectAmplitude = sqrt(objectIntensity);

waveLength = 0.5e-6; %0.5 um

k0 = 2*pi/waveLength;

pixelSize = 0.5e-6;  %0.5 um

NA = 0.1;

cutoffFrequence = NA*k0;

%coherent transfer function

objectAmplitudeFT = fftshift(fft2(objectAmplitude));

[m n] = size(objectAmplitude);

kx = -pi/pixelSize:2*pi/(pixelSize*(n-1)):pi/pixelSize;

ky = -pi/pixelSize:2*pi/(pixelSize*(m-1)):pi/pixelSize;

[kxn, kym] = meshgrid(kx,ky);

CTF = ((kxn.^2 + kym.^2)<cutoffFrequence^2); %uniform frequency response

subplot(222);

imshow(CTF,[]); title('OTF');

%filtering

outputFT = CTF.*objectAmplitudeFT;

subplot(223); imshow(log(abs(outputFT)),[]);

title('Filtered spectrum');

%output

outputAmplitude = ifft2(ifftshift(outputFT));

outputIntensity = abs(outputAmplitude).^2;

subplot(224);

imshow(outputIntensity,[]);

title('Output object');

suptitle('Coherent imaging');

%2) Incoherent imaging

figure(2);

subplot(221); imshow(objectIntensity, []);

title('Input  object');

%incoherent transfer function

cpsf = ifftshift(ifft2(CTF));

ipsf = (abs(cpsf)).^2;

ITF = abs(ifftshift(fft2(ipsf)));

ITF = ITF/max(max(ITF));

subplot(222);

imshow(ITF,[]); title('OTF');

%filtering

outputFT = ITF.*objectAmplitudeFT;

subplot(223); imshow(log(abs(outputFT)),[]);

title('Filtered spectrum');

outputAmplitude = ifft2(ifftshift(outputFT));

outputIntensity = abs(outputAmplitude).^2;

subplot(224);

imshow(outputIntensity,[]);

title('Output object');

suptitle('Incoherent imaging');

References:

https://ocw.mit.edu/courses/mechanical-engineering/2-71-optics-spring-2009/video-lectures/lecture-22-coherent-and-incoherent-imaging/

https://www.osapublishing.org/DirectPDFAccess/E895979C-D52F-0431-77653B998D6C1D5B_330658/FiO-2015-FW4G.1.pdf?da=1&id=330658&uri=FiO-2015-FW4G.1&seq=0&mobile=no

http://www.disp.duke.edu/~dbrady/ece469/notes/naii/naii.html

https://spie.org/samples/FG17.pdf