PSFSHIFT Fix problem with phase direction in Figure 1

Coherent and Incoherent Transfer Functions and Point-Spread Functions for a tilt in the pupil plane

by Chuck DiMarzio Northeastern University December 2009

!! 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

fcenter=1/128;  % Center frequency for bar chart in cycles per pixel
tilt=2*pi/4*(1/fcenter);  % Wavefront tilt 1/4 cycle per fcenter
xax=fftaxisshift(0:1023);
[x,y]=meshgrid(xax,xax);
rsq=x.^2+y.^2;
fax=fftaxisshift(fftaxis(xax));
[fx,fy]=meshgrid(fax,fax);
fsq=fx.^2+fy.^2;
ftop=max(max(fx));fmax=ftop/sqrt(length(fx));
xmax=length(xax)/2/sqrt(length(fx));

Compute the transfer functions and point--spread functions

theshift=floor(1/fcenter/4);  % We want to shift a quarter cycle.
psf1=(x==theshift)&(y==0);
mask1=fftshift(ifft2(fftshift(psf1))); % Go to pupil and generate the mask
                                       % for the displacement
mask1=mask1.*(fsq<fmax^2);  %  Make the aperture
psf1=fftshift(ifft2(ifftshift(mask1)));  % Recompute the actual psf1
psf1a=abs(psf1).^2;
mask1a=fftshift(ifft2(fftshift(psf1a))); % incoherent OTF

Show a slice of the phase through the Transfer Functions

a1=(abs(mask1)>1e-6*max(max(abs(mask1)))).*angle(mask1);
a1a=(abs(mask1a)>1e-6*max(max(abs(mask1a)))).*(angle(mask1a));
fig1=figure;plot(fax,unwrap(a1(1+floor(length(fx)/2),:)),...
            fax,unwrap(a1a(1+floor(length(fx)/2),:)));
grid on;
% why does the phase go in the opposite direction?
moose=axis;axis([fmax*6*[-1,1],moose([3,4])]);
xlabel('f_y, Spatial Freq., pixel^{-1}');
ylabel('PTF(f_x,0)');

Show a slice of the amplitudes through the Transfer Functions

fig2=figure;plot(fax,abs(mask1(1+floor(length(fx)/2),:)/max(max(mask1))),...
            fax,abs(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('MTF(f_x,0)');

Show the PTF with imagesc

fig3=figure;imagesc(fax,fax,a1);axis image;colormap(flipud(gray)); % CPTF
axis(fmax*4*[-1,1,-1,1]);
xlabel('f_x,  pixel^{-1}');
ylabel('f_y,  pixel^{-1}');
colorbar;

fig4=figure;imagesc(fax,fax,a1a);axis image;colormap(flipud(gray)); % IPTF
axis(fmax*4*[-1,1,-1,1]);
xlabel('f_x,  pixel^{-1}');
ylabel('f_y,  pixel^{-1}');
colorbar;

Plot slices of the point--spread functions

fig5=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)');

Now do the imaging: Coherent object

Show object on top and scaled image on bottom for comparison

cobject=cos(2*pi*fcenter*x);
cobjectf=fftshift(ifft2(fftshift(cobject))); % Go to pupil
cimagef=cobjectf.*mask1;
cimage=fftshift(ifft2(ifftshift(cimagef)));
fig6=figure;imagesc(xax,xax,[real(cobject(1:length(fax)/2,:));...
                real(cimage(length(fax)/2+2:end,:))/max(max(abs(cimage)))]);
colormap(flipud(gray));colorbar;
axis image;
xlabel('x, Position, Pixels');
ylabel('x, Position, Pixels');

Incoherent object

iobject=abs(cobject).^2;
iimage=abs(cimage).^2;
fig7=figure;imagesc(xax,xax,[real(iobject(1:length(fax)/2,:));...
                real(iimage(length(fax)/2+2:end,:))/max(max(abs(iimage)))]);
colormap(flipud(gray));colorbar;
axis image;
xlabel('x, Position, Pixels');
ylabel('x, Position, Pixels');