How to build a function in Python (Jupyter) for calculation of Point Spread Function (image processing?)

Question

In Image processing, the image is fourier transformed, then center transformed in fourier, and then under-sampled, and back to inverse center transform and back to image by inverse fourier transform. So, to calculate Point Spread Function, is the below function alright?

def apply_Fu(sampling_pattern, x):
    #compute subsampled FFT    #Sampling pattern is matrix where 1 means we take that sample for undersampling

    return sampling_pattern*np.fft.fftshift(np.fft.fft2(x))

def apply_Fu_adjoint(sampling_pattern, y):
    #Compute adjoing of subsampled k space
    return np.fft.ifft2(np.fft.ifftshift(sampling_pattern*y))

def spr(gridsize, sampling_pattern):
    # gridsize should be 2-element tuple, e.g. gridsize = (10, 10)
    maxima = np.zeros(gridsizedtype = np.complex_)
    for x in range(gridsize[0]):
        for y in range(gridsize[1]):
            # in this iteration, the index "i" corresponds to the gridpoint (x, y)

            # construct basis vector
            e_i = np.zeros(gridsize, dtype = np.complex_)
            e_i[x, y] = 1

            # compute psf_i = Fu* Fu e_i
            # psf_i[xx, yy] is PSF(i,j) if index "j" corresponds to gridpoint (xx, yy)
            psf_i = apply_Fu_adjoint(sampling_pattern, apply_Fu(sampling_pattern, e_i))

            # normalize; psf_i[x, y] is PSF(i,i)
            psf_i = psf_i / psf_i[x, y]

            # trick to exclude point "i" itself from maximum: set it to -infinity
            psf_i[x, y] = -np.inf

            # "inner" maximum, over "j"
            maxima[x, y] = np.max(psf_i)
    spr = np.max(maxima)
    return spr

spr = spr(img.shape, random)

np.abs(spr)```

Answer 1

You can use this https://github.com/cgohlke/psf, a Point Spread Function calculations for fluorescence microscopy

Psf is a Python library to calculate Point Spread Functions (PSF) for fluorescence microscopy.

>>> import psf
>>> args = dict(shape=(32, 32), dims=(4, 4), ex_wavelen=488, em_wavelen=520,
...             num_aperture=1.2, refr_index=1.333,
...             pinhole_radius=0.55, pinhole_shape='round')
>>> obsvol = psf.PSF(psf.GAUSSIAN | psf.CONFOCAL, **args)
>>> print(f'{obsvol.sigma.ou[0]:.5f}, {obsvol.sigma.ou[1]:.5f}')
2.58832, 1.37059
>>> obsvol = psf.PSF(psf.ISOTROPIC | psf.CONFOCAL, **args)
>>> print(obsvol, end='')  # doctest:+ELLIPSIS
PSF
 Confocal, Isotropic
 shape: (32, 32) pixel
 dimensions: (4.00, 4.00) um, (55.64, 61.80) ou, (8.06, 8.06) au
 excitation wavelength: 488.0 nm
 emission wavelength: 520.0 nm
 numeric aperture: 1.20
 refractive index: 1.33
 half cone angle: 64.19 deg
 magnification: 1.00
 underfilling: 1.00
 pinhole radius: 0.550 um, 8.498 ou, 1.1086 au, 4.40 px
 computing time: ... ms
>>> obsvol[0, :3]
array([1.     , 0.51071, 0.04397])
>>> # save the image plane to file
>>> obsvol.slice(0).tofile('_test_slice.bin')
>>> # save a full 3D PSF volume to file
>>> obsvol.volume().tofile('_test_volume.bin')