In [1]:
%pylab inline

%load_ext autoreload
%autoreload 2

from pyotf.otf import SheppardPSF, HanserPSF
from pyotf.utils import easy_fft, easy_ifft

from easy_plot import easy_plot
plt.set_cmap("inferno")


Populating the interactive namespace from numpy and matplotlib
<Figure size 432x288 with 0 Axes>

In [2]:
# We'll use a 1.27 NA water dipping objective imaging in water
psf_params = dict(
    na=1.27,
    ni=1.33,
    wl=0.585,
    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.45 / oversample_factor
psf_params["size"] *= oversample_factor

In [3]:
psf = SheppardPSF(**psf_params).PSFi

In [4]:
# centered coordinate system.
x = (np.arange(psf_params["size"]) - (psf_params["size"] + 1) // 2) * psf_params["res"]

In [5]:
zz, xx = meshgrid(x, x, indexing="ij")

In [6]:
freq = 2 * pi * psf_params["ni"] / psf_params["wl"]
alpha = np.arcsin(psf_params["na"] / psf_params["ni"])
print(freq, np.rad2deg(alpha))
s = np.zeros_like(xx * 1j)
for theta in (-alpha, 0, alpha):
    s += exp(1j*((xx * sin(theta) + zz * cos(theta)) * freq))


14.284848647092053 72.72437394563345

In [7]:
sim_psf = psf * (abs(s)**2)[:, None]

In [8]:
easy_plot((psf, sim_psf), ("Epi", "SIM",), oversample_factor=1, gam=0.5)



In [9]:
zz, yy, xx = x[:, None, None], x[None, :, None], x[None, None, :]

In [10]:
freq = 2 * pi * psf_params["ni"] / psf_params["wl"]
alpha = np.arcsin(psf_params["na"] / psf_params["ni"])
print(freq, np.rad2deg(alpha))
sim_psf = np.zeros_like(psf)
for orientation in (0, 2 * pi / 3, 4 * pi / 3):
    s = np.zeros_like(psf, dtype=complex)
    for strength, theta in zip((1, 0.5, 1), (-alpha, 0, alpha)):
        s += strength * exp(1j*(((xx * cos(orientation) + yy * sin(orientation)) * sin(theta) + zz * cos(theta)) * freq))
    sim_psf += psf * (abs(s)**2)


14.284848647092053 72.72437394563345

In [11]:
otf = abs(easy_fft(psf))**2

In [12]:
w = otf.max() / 100
wiener_otf = otf / (otf + w)
# wiener_otf = otf > 1e-16
wiener_psf = easy_ifft(wiener_otf).real

In [13]:
freq = 2 * pi * psf_params["ni"] / psf_params["wl"]
alpha = np.arcsin(psf_params["na"] / psf_params["ni"])
print(freq, np.rad2deg(alpha))
wiener_sim_psf = np.zeros_like(wiener_psf)
for orientation in (0, 2 * pi / 3, 4 * pi / 3):
    # s = np.zeros_like(xx * 1j)
    s = exp(1j * zz * freq) * np.ones(psf.shape[:2], dtype=complex)
    for theta in (-alpha, alpha):
        s += strength * exp(1j*(((xx * cos(orientation) + yy * sin(orientation)) * sin(theta) + zz * cos(theta)) * freq))
    wiener_sim_psf += wiener_psf * (abs(s)**2)
wiener_sim_psf -= 8 * wiener_psf


14.284848647092053 72.72437394563345

In [14]:
easy_plot((psf, sim_psf, wiener_psf, wiener_sim_psf), ("Epi", "SIM", "Wiener Epi", "Wiener SIM",), oversample_factor=1)



In [15]:
freq = 2 * pi * psf_params["ni"] / psf_params["wl"]
s = exp(1j * zz * freq) * np.ones(psf.shape[1:], dtype=complex)
for orientation in (0, pi / 2):
    rr = xx * cos(orientation) + yy * sin(orientation)
    for theta in (-alpha, alpha):
        s += exp(1j*((rr * sin(theta) + zz * cos(theta)) * freq))
lattice_sim_psf = wiener_psf* (abs(s)**2)

In [16]:
easy_plot((psf, wiener_sim_psf, lattice_sim_psf,), ("Epi", "SIM", "Lattice SIM"), oversample_factor=1, res=psf_params["res"], gam=0.5)