```
In [1]:
```%pylab inline
from easy_plot import easy_plot
import scipy.ndimage as ndi
from scipy.signal import fftconvolve
%load_ext autoreload
%autoreload 2
from pyotf.otf import SheppardPSF, HanserPSF
from pyotf.utils import easy_fft, easy_ifft
import dphplotting as dplt
from dphutils import bin_ndarray

```
```

```
In [2]:
```plt.set_cmap("inferno");

```
```

For any type of fluorescence microscopy technique the total point spread function (PSF) can be expressed as the multiplication of the *excitation* PSF, $PSF_{exc}$, and the *detection* PSF, $PSF_{det}$. Both PSFs are *intensity* PSFs because fluorophores respond to the intensity of the field not the field itself and all detectors are square law detectors meaning they also respond to intensity and not field.

For a perfect confocal microscope, i.e. one with an infinitesimal pinhole (or various flavors of photon reassignment microscopy which act in much the same way with infinitely better efficiency), $PSF_{exc}$ and $PSF_{det}$ are simply widefield PSFs at the excitation and detection wavelengths.

We'll start with this.

```
In [3]:
```# We'll use a 1.27 NA water dipping objective imaging in water
psf_params = dict(
na=1.27,
ni=1.33,
wl=0.561,
size=128,
vec_corr="none"
)
# Set the Nyquist sampling rate
nyquist_sampling = psf_params["wl"] / psf_params["na"] / 4
# our oversampling factor, must be odd for easy integration (for peaked symmetrical funtions
# the kernel must be odd) increase this to have a more accurate simulation
oversample_factor = 1
# we need to be at half nyquist because, for the perfect confocal/airyscan/ISM, we'll get a factor of 2 improvement in OTF cutoff
psf_params["res"] = nyquist_sampling / oversample_factor / 2
psf_params["size"] *= oversample_factor

```
In [4]:
```# calculate 3D
psf_exc = SheppardPSF(**psf_params).PSFi
psf_params["wl"] = 0.585
psf_det = SheppardPSF(**psf_params).PSFi

```
In [5]:
```# total PSF is multiplication of excitation and emission
psf_con = psf_exc * psf_det

```
In [6]:
```labels = "Excitation", "Detection", "Confocal"
psfs = psf_exc, psf_det, psf_con
easy_plot(psfs, labels, oversample_factor, res=psf_params["res"])

```
```

When speaking about confocal pinholes we usual work in what's called *Airy Units* or *AU*. One AU is equivalent to the diamter of the first airy ring, which is defined as the first minimum in the airy function. That means that

Modeling the effects of the pinhole are actually quite simple: we just need to convolve the detection PSF in the lateral plane with the pinhole transmission function, which in the simplest (and most reasonable case) is simply a ciruclar top hat function. Why convolution? If you were to image a single point emitter perfectly aligned with the excitation PSF then you'd just multiply by the pinhole and integrate. But what if your emitter is off-axis? Then it's image is also off axis etc.

```
In [7]:
```# define airy units in terms of pixels
airy_unit = 1.22 * psf_params["wl"] / psf_params["na"] / psf_params["res"]
airy_unit

```
Out[7]:
```

```
In [8]:
```def disk_kernel(radius):
full_size = int(np.ceil(radius * 2))
if full_size % 2 == 0:
full_size += 1
coords = np.indices((full_size, full_size)) - (full_size - 1) // 2
r = np.sqrt((coords**2).sum(0))
kernel = r < radius
return kernel / kernel.sum()

```
In [9]:
```kernel = disk_kernel(1 * airy_unit / 2)
psf_det_au1 = fftconvolve(psf_det, kernel[None], "same", axes=(1,2))
psf_con_au1 = psf_det_au1 * psf_exc

```
In [10]:
```easy_plot((psf_exc, psf_det_au1, psf_con_au1), labels, oversample_factor, psf_params["res"])

```
```

```
In [11]:
```# now lets try 10 AU, 1 AU and 0.1 AU
psf_au = []
label_au = []
for au in (0.1, 1, 10):
kernel = disk_kernel(au * airy_unit / 2)
psf_det_au = fftconvolve(psf_det, kernel[None], "same", axes=(1,2))
psf_au.append(psf_det_au * psf_exc)
label_au.append(f"{au:} AU")

```
In [12]:
```easy_plot([psf_con] + psf_au, ["0 AU"] + label_au, oversample_factor, psf_params["res"])

```
```