INCPSF Coherent and Incoherent Transfer Functions and Point-Spread

Functions

SYNTAX: incpsf;

by Chuck DiMarzio Northeastern University December 2009

Requires: FFTAXISSHIFT,FFTAXIS

!! This file may be copied, used, or modified for educational and
!! research purposes provided that this header information is not
!! removed or altered, and provided that the book is cited in
!! publications, as DiMarzio, Charles A., Optics for Engineers,
!! CRC Press, Boca Raton, FL, 2011.
!! http://www.crcpress.com
!! Other distribution is prohibited without permission.

Contents

Setup

xax=fftaxisshift(0:1023); % set up axes with zero center
[x,y]=meshgrid(xax,xax);  % generate coordinates
rsq=x.^2+y.^2;  % Compute radius squared
fax=fftaxisshift(fftaxis(xax));  % generate frequency axis
[fx,fy]=meshgrid(fax,fax);  % generate frequency coordinates
fsq=fx.^2+fy.^2;  % Magnitude of frequency squared
ftop=max(max(fx));fmax=ftop/sqrt(length(fx));  % for mask size
xmax=length(xax)/2/sqrt(length(fx));

Compute the transfer functions

mask1=single(fsq<(fmax^2));  % Low Pass filter for coherent transfer
                             % function.
%mask1a=fftshift(ifft2(abs((fft2(fftshift(mask1)))).^2));
% The incoherent transfer function is the autocorrelation function of the
% coherent transfer function.

Compute the point--spread functions

psf1=fftshift(ifft2(ifftshift(mask1)));
%psf1a=fftshift(ifft2(ifftshift(mask1a)));
psf1a=abs(psf1).^2;
mask1a=fftshift(ifft2(fftshift(psf1a)));

Show a slice through the Transfer Functions

fig1=figure;plot(fax,real(mask1(:,1+floor(length(fx)/2))/max(max(mask1))),...
            fax,real(mask1a(:,1+floor(length(fx)/2)))/max(max(mask1a)));
grid on;
moose=axis;axis([fmax*4*[-1,1],moose([3,4])]);
xlabel('f_x, Spatial Freq., pixel^{-1}');
ylabel('Transfer Function(f_x,0)');
% Check that the imaginary part is zero
figure;plot(fax,imag(mask1(:,1+floor(length(fax)/2))/max(max(mask1a))),...
            fax,imag(mask1a(:,1+floor(length(fax)/2)))/max(max(mask1a)));
moose=axis;axis([fmax*4*[-1,1],moose([3,4])]);
grid on;
xlabel('f_x, Spatial Freq., pixel^{-1}');
ylabel('Transfer Function');
title('Sanity Check: Imaginary parts zero');

Show the transfer functions with imagesc

fig2=figure;imagesc(fax,fax,mask1);axis image;colormap(flipud(gray));
axis(fmax*4*[-1,1,-1,1]);
xlabel('f_x, Spatial Freq., pixel^{-1}');
ylabel('f_y, Spatial Freq., pixel^{-1}');

fig3=figure;imagesc(fax,fax,real(mask1a));axis image;colormap(flipud(gray));
axis(fmax*4*[-1,1,-1,1]);
xlabel('f_x, Spatial Freq., pixel^{-1}');
ylabel('f_y, Spatial Freq., pixel^{-1}');

Show the point--Spread functions with imagesc

fig4=figure;imagesc(xax,xax,20*log10(abs(psf1)));axis image;
axis(xmax*8*[-1,1,-1,1]);
moose=caxis;caxis(moose(2)+[-30,0]);colormap(flipud(gray));colorbar;
xlabel('x, Position, Pixels');
ylabel('y, Position, Pixels');

% Check angle to make sure psf is real
figure;imagesc(xax,xax,angle(psf1));colormap('hsv');axis image;colorbar;
axis(xmax*8*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
ylabel('y, Position, Pixels');
title('Sanity Check: Angle zero or 180');

% Now the incoherent one
fig5=figure;imagesc(xax,xax,10*log10(abs(psf1a)));axis image;
axis(xmax*8*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
ylabel('y, Position, Pixels');
moose=caxis;caxis(moose(2)+[-30,0]);colormap(flipud(gray));colorbar;

% Check angle to make sure psf is real
figure;imagesc(xax,xax,angle(psf1a));colormap('hsv');axis image;colorbar;
axis(xmax*8*[-1,1,-1,1]);
xlabel('x, Position, Pixels');
ylabel('y, Position, Pixels');
title('Sanity Check: Angle zero or 180');

Plot slices of the point--spread functions

fig6=figure;plot(xax,real(psf1(:,1+floor(length(fax)/2))/max(max(psf1))),...
            xax,real(psf1a(:,1+floor(length(fax)/2)))/max(max(psf1a)));
moose=axis;axis([xmax*8*[-1,1],moose([3,4])]);
grid on;
xlabel('x, Position, Pixels');
ylabel('PSF(x,0)');