# Point spread function (PSF) calculation (diffraction limited)

The point spread function (PSF) is a basic characterization of spatial image resolution. The PSF of a diffraction limited optics defines an upper limit on resolution. This script illustrates the PSF calculations embedded in oiPlot/OTF.

Copyright ImagEval Consultants, LLC, 2012.

## Not needed, but often convenient.

```ieInit
```

## Create an optical image and optics with a large f number

```oi     = oiCreate;
optics = oiGet(oi,'optics');

% This large is a little blurry
fNumber = 12;
optics = opticsSet(optics,'fnumber',fNumber);
oi     = oiSet(oi,'optics',optics);
```

## This is the oiPlot/OTF code that converts the OTF to PSF conversion

```% Specify units
units = 'um';
nSamp = 100;   % Number of frequency steps from 0 to incoherent cutoff
thisWave = 600;

% Get the diffraction limited frequency support
val = opticsGet(optics,'dlF Support',thisWave,units,nSamp);
clear fSupport
[fSupport(:,:,1),fSupport(:,:,2)] = meshgrid(val{1},val{2});
```

## Oversampling

```%  We over sample the frequency to get a smoother PSF image. You can
%  specify the factor for oversampling if you like in the calling
%  arguments to oiPlot.
s = 4;
fSupport = fSupport*s;

% Frequency units are cycles/micron The spatial frequency support runs from
% -Nyquist:Nyquist. With this coding, the Nyquist frequency is actually the
% peak frequency value. There are two samples per Nyquist, so the sample
% spacing is 1/(2*peakF)
deltaSpace = 1/(2*max(fSupport(:)));
```

## Calculate the diffraction limited MTF

```otf = dlMTF(oi,fSupport,thisWave,units);
psf = fftshift(ifft2(otf));

samp = (-nSamp:(nSamp-1));
[X,Y] = meshgrid(samp,samp);
clear sSupport
sSupport(:,:,1) = X*deltaSpace;
sSupport(:,:,2) = Y*deltaSpace;
```

## Determine first zero crossing

```fNumber = opticsGet(optics,'f number');
radius = (2.44*fNumber*thisWave*10^-9)/2 * ieUnitScaleFactor(units);

x = sSupport(:,:,1); y = sSupport(:,:,2);
```

## Plot the PSF. Zoom and rotate the image to see the details.

```vcNewGraphWin;
mesh(x,y,psf);

% For the diffraction limited case, we draw the Airy disk 