Revisiting PSF models: unifying framework and high-performance implementation

TL;DR

The paper addresses fragmented PSF models for high-NA microscopy by introducing a Richards-Wolf–based unifying framework that links Fourier (Cartesian) and Bessel (spherical) approaches via a shared propagation integral. It provides correction factors for Gibson-Lanni aberrations, arbitrary phase aberrations, apodization, and Gaussian envelopes, and implements these within a high-performance, differentiable PyTorch library (psf-generator) for CPU/GPU use. The authors demonstrate that the Bessel (spherical) approach is optimal for axisymmetric pupils, while the Fourier (Cartesian) approach is more versatile for non-axisymmetric configurations, enabling efficient integration into learning-based pipelines. Comprehensive benchmarks show the trade-offs in speed and accuracy across propagators and hardware, and the work includes open-source interoperability with napari and Chromatix, facilitating widespread adoption and integration into simulation and optimization workflows.

Abstract

Localization microscopy often relies on detailed models of point spread functions. For applications such as deconvolution or PSF engineering, accurate models for light propagation in imaging systems with high numerical aperture are required. Different models have been proposed based on 2D Fourier transforms or 1D Bessel integrals. The most precise ones combine a vectorial description of the electric field and precise aberration models. However, it may be unclear which model to choose, as there is no comprehensive comparison between the Fourier and Bessel approaches yet. Moreover, many existing libraries are written in Java (e.g. our previous PSF generator software) or MATLAB, which hinders the integration into deep learning algorithms. In this work, we start from the original Richards-Wolf integral and revisit both approaches in a systematic way. We present a unifying framework in which we prove the equivalence between the Fourier and Bessel strategies and detail a variety of correction factors applicable to both of them. Then, we provide a high-performance implementation of our theoretical framework in the form of an open-source library that is built on top of PyTorch, a popular library for deep learning. It enables us to benchmark the accuracy and computational speed of different models, thus allowing for an in-depth comparison of the existing models for the first time. We show that the Bessel strategy is optimal for axisymmetric beams while the Fourier approach can be applied to more general scenarios. Our work enables efficient PSF computation on CPU or GPU, which can then be included in simulation and optimization pipelines.

Figures (7)

• Figure 1: Geometry of the focusing of optical fields. An incident field $\mathbf{e}_\mathrm{inc}$ is transformed by a focusing element into a converging spherical wave $\mathbf{e}_\infty$. These fields are parameterized by a unit vector $\mathbf{s}$, either in Cartesian coordinates $(s_x, s_y, s_z)$ or spherical coordinates $(\theta, \phi)$. The focus field $\mathbf{E}$ is parameterized by $\boldsymbol{\rho} = (x, y, z)$, rewritten here to introduce cylindrical coordinates $\rho$ and $\phi$.
• Figure 2: Correction factors and their associated physical origins. Top part: phase correction factors, either introduced on the incident field or due to refraction. Bottom part: amplitude correction factors that model the incident beam envelope or the apodization factor for energy conservation.
• Figure 3: Graphical user interface of our PSF simulation plugin integrated with napari (v0.4.0, 2025-06-17). Left: napari viewer panel enabling interactive 3D navigation and rendering options. Middle: interactive visualization area, displaying here a slice of a high-NA 3D PSF. Right: plugin interface for selecting the propagator model, specifying physical and numerical parameters (e.g., NA, wavelength, Zernike aberrations, CPU or GPU backend), and exporting results.
• Figure 4: Unaberrated PSFs generated by the scalar ((a)-(d)) and vectorial ((e)-(h)) models, in the case of low-NA (0.5) (first two rows) and high-NA (1.3) (last two rows). (a), (e), (c), and (g): Slice of the $x$-$y$ plane. (b), (f), (d), and (h): Slice of the $z$-$y$ plane. (i)-(l): Intensity profiles along a vertical line through the center of the PSFs. The $x$-axis shows the relative distance to the center of the image in micrometers. Scale bars represent 3 and 0.6 in the low- and high- NA cases, respectively. For images of $z$-$y$ planes ((b), (f), (d), and (h)), scale bars for the $y$ and $z$ axes are indicated by the horizontal and vertical bars, respectively.
• Figure 5: PSFs with phase aberrations generated by the vectorial Cartesian propagator at high-NA (1.3). (a), (d), and (g): Introduced phase masks at the pupil plane. (b), (e), and (h): Slice of the $x$-$y$ plane at $z=0$. (c), (f), and (i): Slice of the $z$-$x$ plane at $y=0$. (j)-(l): PSF with vertical astigmatism at three z-planes (in focus (k), 500 above (j), and 500 below focus (l)). All PSFs were generated with circular polarization, except the Half-Moon PSF where x-axis linear polarization is being used. The scale bars in all images represent 0.6.