``````

In [1]:

%pylab inline

from easy_plot import easy_plot

from pyotf.otf import SheppardPSF, HanserPSF

``````
``````

Populating the interactive namespace from numpy and matplotlib

``````
``````

In [2]:

plt.set_cmap("inferno");

``````
``````

<Figure size 432x288 with 0 Axes>

``````

## Difference between models

Because the `HanserPSF` is a mixed frequency space and real space model, with the lateral dimensions being modeled in frequency space and the axial dimension in real space, the results are slightly different from the `SheppardPSF`. In particular, the `HanserPSF` model obeys conservation of energy, meaning that the sum of all values for each \$z\$ plane is constant. While this is correct in theory it isn't correct in practice because of the finite extent of the simulation. What happens is when the lateral extent is too large the energy get's reflected back towards the center.

Let's look at these directly

``````

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)
oversample_factor = 1

# Over sample to
psf_params["res"] = nyquist_sampling * 0.99 / oversample_factor
psf_params["size"] *= oversample_factor

``````
``````

In [4]:

# calculate 3D
psf_s = SheppardPSF(**psf_params)
psf_h = HanserPSF(**psf_params)

``````
``````

In [5]:

for p, l in zip((psf_s.PSFi, psf_h.PSFi), ("Sheppard", "Hanser")):
p = p/p.max()
plot(p.sum((1,2)), label=l)
legend()

``````
``````

Out[5]:

<matplotlib.legend.Legend at 0x1a1f9207f0>

``````

Clearly the `SheppardPSF` model has decreasing energy as the PSF gets spread out. What happens if we normalize it?

``````

In [6]:

# plot
easy_plot((psf_s.PSFi, psf_s.PSFi/psf_s.PSFi.sum((1, 2), keepdims=True), psf_h.PSFi), ("Sheppard", "Sheppard Norm", "Hanser"), res=psf_params["res"])

``````
``````

``````
``````

In [7]:

psf_params["zsize"] = 128

for s in (32, 64, 128, 256):
psf_params["size"] = s
p = SheppardPSF(**psf_params).PSFi
p = p.mean((1,2))
p = p/p.max()
plot(p, label=f"{s:} pixels")
legend()

``````
``````

Out[7]:

<matplotlib.legend.Legend at 0x1a3533a2b0>

``````
``````

In [8]:

psf_params["zsize"] = 128 * 3
psf_params["size"] = 128
psf_norm = SheppardPSF(**psf_params).PSFi
psf_norm /= psf_norm.sum((1,2), keepdims=True)
psf_norm = psf_norm[128:-128]

``````
``````

In [9]:

# plot
easy_plot((psf_s.PSFi, psf_s.PSFi/psf_s.PSFi.sum((1, 2), keepdims=True), psf_norm, psf_h.PSFi), ("Sheppard", "Sheppard Norm", "Sheppard Norm2", "Hanser"), res=psf_params["res"])

``````
``````

``````