Revisiting PSF models: Unifying framework and high‐performance implementation
Yan Liu, Vasiliki Stergiopoulou, Jonathan Chuah, Eric Bezzam, Gert‐Jan Both, Michael Unser, Daniel Sage, Jonathan Dong
TL;DR
This paper introduces a unified framework and high-performance tool for computing point-spread functions in microscopy, enabling accurate and efficient simulations.
Contribution
A unified theoretical framework and PyTorch-based implementation that compares and optimizes PSF models for microscopy.
Findings
The Bessel strategy is optimal for axisymmetric beams, while the Fourier approach is more versatile.
The new PyTorch-based library enables efficient PSF computation on CPUs and GPUs.
The framework allows for accurate benchmarking and comparison of existing PSF models.
Abstract
Localisation 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 a 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 accurate 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 their 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…
Genes, proteins, chemicals, diseases, species, mutations and cell lines named across the full text — each resolved to its canonical identifier and authoritative record.
Click any figure to enlarge with its caption.
FIGURE 1
FIGURE 2
FIGURE 3
FIGURE 4
FIGURE 5
FIGURE 6
FIGURE 7- —Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung10.13039/501100001711
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Electron Microscopy Techniques and Applications · Digital Holography and Microscopy · Advanced Fluorescence Microscopy Techniques
INTRODUCTION
1
The point‐spread function (PSF), also referred to as the impulse response, is a fundamental concept that encapsulates key features of an optical microscope. Its monitoring and assessment have been longstanding routines in microscopy, supported by the development of specific tools1, 2, 3 and the creation of a consortium to standardise best practices.4 A detailed characterisation of the PSF is essential to the design of computational imaging processes such as single‐molecule localisation microscopy5 and super‐resolution microscopy, including 3D deconvolution microscopy,6, 7 structured illumination microscopy (SR‐SIM),8, 9 fluctuation microscopy,10, 11 stimulated emission‐depletion microscopy (STED),12 and MINFLUX.13 The PSF is at the heart of these localisation techniques and helps one to achieve a spatial resolution beyond the optical diffraction limit.
The PSF can be measured experimentally from a z‐stack of tiny fluorescent beads.6, 14 However, a theoretical PSF is often required. Here, the challenge lies in the accurate model of light propagation for imaging systems with high numerical aperture (NA). Detailed models have been proposed15, 16, 17, 18 to (a) go beyond simple paraxial approximations, (b) take into account the vectorial nature of the electric field, and (c) include various aberration factors such as the Gibson–Lanni aberrations due to refraction at planar interfaces.19 Several software packages have been developed to generate precise theoretical 3D PSFs that include various features and aberrations. Among these, the Huygens PSF software1 is well‐known in the microscopy community, and our Java‐based ImageJ/FIJI plugin, PSF Generator,20 has also been widely used in the past. Other open‐source alternatives include MATLAB‐based solutions3, 21, 22 as well as Python‐based tools.23, 24
Computational PSF models have often been applied in fluorescence microscopy. For example, the fitting of a theoretical PSF to a measured one can improve accuracy.20, 25 Detailed models enable the quantification of the 3D uncertainty,26 while the engineering of PSFs may take full advantage of parametric models.27, 28, 29 They are a key component of high‐quality single‐molecule localisation microscopy (SMLM) images,30, 31 including single‐particle tracking32, 33 and virtual SMLM microscopes.34, 35
Existing approaches to the computation of the PSF can be categorised into two classes of models: those based on 2D Fourier transforms3, 16 or 1D integrals of Bessel functions.15, 17, 18, 36 These two classes rely on different assumptions (with Bessel models being more restrictive) and bring different computational trade‐offs. Most works focus on one of the two approaches, and the relationship between them remains unclear. For example, the Gibson–Lanni aberrations19 or the apodisation factor15 are only applied in the Bessel case, although they could be generalised to both. To the best of our knowledge, no systematic benchmark of the two strategies has been published, either in terms of accuracy or computational speed.
In this work, we propose a unifying framework for PSF models in which we show the equivalence between the Fourier and the Bessel approaches as being two separate parameterisations of the same propagation integral. This enables us to generalise diverse correction factors and apply them to both models. Eventually, the choice is simplified to scalar versus vectorial models, optionally with additional correction factors. We distribute an open‐source PyTorch‐based library called psf‐generator. It inherits all the functionalities of PyTorch and allows for a seamless integration with modern learning‐based algorithms. We benchmark the two classes of models on CPU and GPU and show that the Bessel model is advantageous for axisymmetric PSFs while the Fourier one is more general and thus more suitable for applications such as PSF engineering. To support the open‐source community and facilitate the adoption of high‐NA PSF models, we also provide a graphical interface in the napari ecosystem37 and integrate our work in the chromatix optical simulation framework.38
BACKGROUND
2
Electromagnetic waves in optical systems are fundamentally described by the Helmholtz wave equations, yet solving these equations in full generality is computationally intractable. Various approximations and integral formulations have thus been developed. The cornerstone for precise PSF models in high‐NA systems is the Richards–Wolf integral,15 which can be viewed as a vectorial extension of the Debye integral.39 The conditions of validity for the Richards–Wolf integral have been thoroughly discussed by Wolf and Li.40 While alternative light propagation models exist, the Huygens–Fresnel approach has been shown to be equivalent, to some extent, for PSF calculations.41
This has been the basis of a line of work for PSF models based on Bessel functions, which we will later refer to as the spherical parameterisation of the Richards–Wolf integral. This approach includes simpler formulations like the Kirchhoff model and more sophisticated vectorial representations.17, 18 These models have progressively been refined by incorporating various correction factors to account for additional physical processes: the Gibson–Lanni model for spherical aberrations due to refractive index mismatch19 (later generalised in Refs. 36, 42), apodisation factors for energy conservation,15 and Fresnel transmission coefficients for accurate interface modelling.17 Models using this spherical parameterisation have been implemented in various software libraries in Java20 and Python,23 which makes it widely accessible to researchers, albeit with some limitations in computational efficiency and integration with deep‐learning frameworks.
Another line of work on PSF modelling is based on Fourier transforms, both in scalar43 and vectorial16 formulations. These models are based on a Cartesian parameterisation of the underlying Richards–Wolf integral and they represent a more general counterpart of the spherical parameterisation. Recently, these high‐NA Fourier models have been implemented in MATLAB3 and Tensorflow as part of a PSF fitting library.44 Adequate sampling of the Fourier transform is crucial for obtaining high‐resolution PSFs and avoiding aliasing. A common trick based on the chirp Z transform is usually implemented to achieve it.3, 16, 29
THEORY
3
The Richards–Wolf model
3.1
As depicted in Figure 1, one obtains the PSF by computing the propagation of light after going through a focusing element, typically a microscope objective or lens. The incident field einc(s), also called the pupil function, is represented by a disk with a maximal cutoff angle defined by the NA of the imaging system. The focusing element transforms the planar incident field into a spherical wave, e∞(s), evaluated on the Gaussian reference sphere. This corresponds to an ensemble of far fields that propagate with direction s and converge to the focal point O. Our goal is to compute the focused electric field E(ρ) near O. Due to the reciprocity of light propagation, this model can also be extended to the emission of a point source to the back focal plane of a microscope objective.
Geometry of the focusing of optical fields. An incident field einc is transformed by a focusing element into a converging spherical wave e∞. These fields are parameterised by a unit vector s, either in Cartesian coordinates (sx,sy,sz) or spherical coordinates (θ,ϕ). The focus field E is parameterised by ρ=(x,y,z), rewritten here to introduce cylindrical coordinates ρ and ϕ.
The model introduced by Ref. (37) is the starting point that allows us to derive all the precise PSF models described in previous works. The focal field is given by a sum of plane waves with direction s=(sx,sy,sz), as expressed by
where ρ=(x,y,z)=(ρcosφ,ρsinφ,z) is the position vector in the focal region of the lens, s=(sx,sy,sz)=(sinθcosϕ,sinθsinϕ,cosθ) is a unit vector that describes the direction of an incoming ray, f is the focal length of the lens, k=2πnλ is the wavenumber, λ is the wavelength, n is the refractive index of the propagation medium, and e∞(s) describes the field distribution on the Gaussian reference sphere. We integrate over the set Ω of solid angles defined on a region sx2+sy2≤smax2, where smax=NAni is the cut‐off determined by the NA and ni is the refractive index of the immersion medium. The angle θ is therefore defined within the immersion medium. We first describe the different classes of models based on this simplified concise equation and postpone to Section 3.4 the introduction of correction factors.
Scalar models
3.2
As a first step, it is common to rely on a scalar approximation to simplify calculations, especially in low‐NA scenarios. In this case, the far field is equal to the incident field, so that e∞(s)=einc(s). The focal field is then given by
This expression involves a two‐dimensional integral over the pupil disk. Two specific parameterisations yield the two classes of models described previously.
The Cartesian parameterisation utilises both sx and sy coordinates with dΩ=dsxdsy/sz, which results in
In this form, the focused field at a z transverse plane is given by the 2D inverse Fourier transform of (e∞(sx,sy)expikszz/sz), where sz=1−sx2−sy2. Thus, the Cartesian parameterisation of the Richards–Wolf integral leverages the speed and efficiency of the fast Fourier transform (FFT) algorithm.
Alternatively, the spherical approach parameterises the problem with two angles θ∈[0,θmax] (the maximum angle θmax is determined by the NA) and ϕ∈[0,2π], as depicted in Figure 1. With ρ=(x,y,z)=(ρcosφ,ρsinφ,z) and dΩ=sinθdθdϕ, the field in the focal region can be rewritten as
Equation (4) can be further simplified if one assumes that the pupil function is axisymmetric (rotational invariant), in the sense that e∞(θ,ϕ)=e∞(θ). In this case, the integral over ϕ in (4) can be computed explicitly using the Bessel function J0 2, which leads to:
Defocus is included in these models via the defocus phase factor expikszz=expikzcosθ where z is the defocus distance. This expression, also known as angular spectrum propagation,43 accurately models the propagation of an electric field in a homogeneous medium.
Vectorial models
3.3
The electric field is a vectorial quantity. Consequently, vectorial propagation models are necessary to accurately account for the propagation and crosstalk between the components of the vector field. Precise vectorial models are crucial for high‐NA systems, where the need to consider high angles arises.
In the vectorial model, the dependence on the incident field einc(s) of the far field e∞(s) now requires us to perform the basis change from a cylindrical to a spherical coordinate system, as in
where einc=[eincx,eincy,0]. The so‐called Fresnel transmission coefficients qs and qp have been introduced to account for the partial reflections at interfaces, which depend on the polarisation state and incidence angle. For each polarisation, the coefficients correspond to the product of all transmission coefficients for each interface from medium m to m+1, with
It is worth noting that in the low‐NA limit (θ=0) with qs=qp=1, (6) simplifies to e∞(θ,ϕ)=einc(θ,ϕ). This explains why the scalar model did not account for this geometric change of variables.
The Cartesian parameterisation of the vectorial model consists of the integral
which essentially boils down to the computation of the inverse Fourier transform of (e∞(sx,sy)expikszz/sz), similar to the scalar case.
With the help of coordinate transformations similar to the scalar case, we derive the spherical parameterisation of the field in the focal region as
Inserting (6) into (9) and using the axisymmetric assumption of the incident field, we obtain a simplified expression of the focal field as
where
with a∈{x,y}.
Correction factors
3.4
Precise PSF models commonly take into account several physical effects which may affect the true PSF. They are added as amplitude factors a(s) and phase factors W(s) in the original integral over solid angles in (1), leading to
Equation (12) enables us to express these correction factors in full generality for both the vectorial and scalar models and for both the Cartesian and spherical parameterisations. We present a detailed list of these correction factors in Section 3.4.1, 3.4.2, 3.4.3, 3.4.4 and their graphical description in 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 apodisation factor for energy conservation.
Gibson–Lanni aberrations
3.4.1
Microscopes typically have stratified layers of different refractive indices. The biological sample is usually aqueous, on top of which we place a coverslip made of glass, and the whole sample is then put in a water or oil immersion medium to increase the numerical aperture. The microscope objectives are designed to provide aberration‐free images in a specific setting with design values for refractive indices and thicknesses of each layer. Any mismatch introduces spherical aberrations due to refraction at the interface between the layers. These aberrations can be computed using the formula3
where ns, ni, ng are the refractive indices of the sample, immersion medium, and glass, respectively, ts, ti, tg are the thicknesses of the sample, immersion medium, and glass respectively, and their counterparts with stars are the design conditions.
In practice, it is challenging to assess the thickness of the immersion medium. Since this distance is manually tuned to obtain an optimal focus of a point emitter at depth ts on the camera, this focusing condition gives the relation
which can be inserted in 13. This particular expression has first been derived for the spherical scalar case in Ref. (13) and extended to the spherical vectorial case in Refs. (51, 52).
Arbitrary phase aberrations
3.4.2
Refined aberration models can be introduced to describe imperfections in the optical system, for the benefit of PSF engineering or wavefront shaping. Designed aberrations can be introduced purposefully via a phase mask or a spatial light modulator at the pupil plane. The aberrations are often parameterised by Zernike polynomials (a set of orthonormal polynomials defined on the pupil disk) or by a direct fixed phase mask to obtain desired PSFs. We write it in the most general case as
Equation (15) is composed of the sum of two terms. The first term is an inner product of the first K Zernike polynomials and their corresponding coefficients ck. The second term can be used to include special cases not covered by the Zernike polynomials (e.g., a vortex phase ramp that leads to a donut PSF as is typically used in STED).
Since the arbitrary phase aberrations described in Equation (15) may not necessarily be axisymmetric, they can only be applied to the most general, Cartesian parameterisation.
Apodisation factor
3.4.3
The apodisation factor is an amplitude‐correction factor that ensures energy conservation during the change of basis from cylindrical coordinates (incident field einc) to spherical coordinates (far field e∞), which takes prominence especially for high‐NA objectives. Since areas of cross‐sections are modified, the field is also rescaled accordingly. Such rescaling ensures that the differential areas dA1 on the plane and dA2 on the sphere, as shown in Figure 2, remain consistent under the change of coordinates. The correcting factor is
when going from cylindrical to spherical in the focusing configuration of Figure 1.
Gaussian envelope
3.4.4
The incident illumination can also depart from a perfect uniform plane wave. In particular, we often assume a Gaussian envelope and expressed it as
where the constant senv determines the size of the envelope.
IMPLEMENTATION
4
PyTorch library
4.1
We provide a high‐performance open‐source PyTorch library psf‐generator 4 to generate 2D and 3D PSFs. The library implements the four PSF models described in Section 3:
- ScalarCartesianPropagator (3);
- ScalarSphericalPropagator (5);
- VectorialCartesianPropagator (8);
- VectorialSphericalPropagator (10). Users can choose between these propagator types, as well as configure physical (e.g., numerical aperture, wavelength, and field of view), and numerical parameters (e.g., image dimensions and number of z‐planes). Our library also allows the users to freely apply any kind of correction factors tailored to their microscope on any chosen propagator as described in Section 3.4.2. The propagators use these parameters to first define the far field (pupil) and propagate it to obtain the focus field (PSF). Finally, the user can visualise, save, and export the generated PSF using our built‐in utility functions. Written in PyTorch, the library easily integrates into deep‐learning workflows, leveraging native features of PyTorch, such as automatic differentiation.
In our unifying framework, these models encompass all previously proposed PSF models based on the Richards–Wolf integral. As discussed in Section 3.2, the spherical propagators require the axisymmetric assumption. When it is valid, the Cartesian and spherical propagators perform equivalent integral computations and result in the same PSF. Both are proposed in our library as they differ in terms of computational efficiency and applicability.
Any tensor corresponding to a field has a shape of the form ( z , channel, x , y ) to benefit from PyTorch performance optimisation, both for data loading and computation parallelisation. Computation are performed plane‐by‐plane on the axial z axis, while x and y correspond to the transverse sizes of 2D images. Akin to greyscale versus RGB images, the channel dimension is equal to 1 for scalar models and 3 in the vectorial case. We use the ZernikePy library5 to generate Zernike polynomials.
The Cartesian parameterisation relies on multiple calls of 2D Fourier transforms with a desired field of view and pixel size. To account for the very small physical dimensions associated to pixels in localisation microscopy, we have implemented a custom variant of the 2D FFT in PyTorch.3, 16, 29 It enables arbitrary pixel sizes and is based on the chirp Z‐transform, at the cost of a convolution and three FFTs. The computational complexity of a single plane is still O(nlogn), with n the size of the transverse plane. The size of the Fourier transform is doubled as we only keep the valid convolution range.
The spherical parameterisation relies on the fast and accurate evaluation of multiple one‐dimensional integrals that involve the Bessel function of first order J0 over the interval [0,θmax], as defined in Equation (5) for the scalar case. To maximise speed, we vectorise the computation of the batch of 1D integrals at different defocus distances via torch.vmap. Then, for each 1D integral, we adopt the composite Simpson's rule45 to benefit from its high accuracy, good numerical stability, and minimal computational overhead. The computational complexity in this case is O(n), with n the number of integration steps in the interval [0,θmax]. As automatic differentiation of Bessel functions is not natively supported by PyTorch as of version 2.3, we also provide a differentiable version of the Bessel functions in our library.
Open‐source interoperability
4.2
Motivated by FAIR principles for data and code management,46 we not only open‐source our implementations but make them interoperable with other tools. Our goal is to further support the dissemination and adoption of accurate physical models for high‐NA light focusing. We present two complementary examples: an intuitive and rapid visualisation tool aimed at end users and an integration with a comprehensive optical simulation library intended for researchers.
napari is an open‐source, Python‐based image viewer designed for interactive visualisation and analysis of multidimensional biological and scientific image data. Built for researchers and developers working with large imaging data from optical microscopy and other imaging modalities, it provides a high‐performance interface for exploring, annotating, and processing images. napari is also extensible via plugins, benefiting from the scientific Python ecosystem.
We provide an interface to generate and visualise 3D PSFs using our PyTorch‐based framework, as depicted in Figure 3. Available as a plugin on the napari hub,6 users can select between the four different propagation models and configure a set of physical and numerical parameters through a user‐friendly graphical interface. This facilitates the introduction of high‐NA correction factors to all four propagators types. PSFs are computed in real‐time and rendered instantly, with support for both CPU and GPU backends. Our plugin allows immediate visual feedback and options to save the resulting volumes. It is designed to provide a code‐free solution for the generation, visualisation, and export of high‐NA PSFs for downstream applications.
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 visualisation 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.
Chromatix is an open‐source wave optics library capable of simulating a wide range of optical systems such as light field imaging, holography and quantitative phase imaging.38 Built on JAX, it supports GPU acceleration and full end‐to‐end differentiability, making it particularly suitable for large‐scale inverse problems and integration with deep learning applications. Compared to PyTorch, JAX offers more seamless support for automatic differentiation and just‐in‐time (JIT) compilation, which can be advantageous for optimising physics‐based simulations based on wave optics.
To ensure maximal performance and seamless integration with the rest of the library, the high‐NA focusing models have been reimplemented in JAX. We focused on the vectorial Cartesian parameterisation, which offers greater generality than the spherical formulation and supports arbitrary pupil fields and aberrations. This functionality is introduced via the high_na_ff_lens() function, which complements the existing ff_lens() implementation used to model the Fourier transform of low‐NA lenses. This addition enables Chromatix to accurately model high‐NA microscope systems and facilitates integration into existing simulation pipelines.
These two examples highlight the importance of actively promoting interoperability within the open‐source ecosystem. Nevertheless, all figures and benchmarks presented in this work were produced based on our comprehensive PyTorch‐based library, which serves as the foundation of our study.
RESULTS
5
PSF gallery
5.1
We provide in Figures 4 and 5a gallery of PSFs to showcase our PyTorch implementation. We used a wavelength of 632 nm, a circular polarisation for the vectorial models, and we display the amplitude (in arbitrary units) of the beams for better contrast. The 2D slices of the same 3D beam share the same dynamic range. In Figure 4, we present the focal electric field computed using both the scalar and vectorial models for objectives with different NAs. In the low‐NA case (top two rows), the resulting PSFs from both propagators are similar. The beams differ more from each other, as expected, in the high‐NA case (bottom two rows): the rings are blurred out as the energy is spread into different components of the focus field.
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, g) Slice of the x‐y plane. (b, f, d, h) Slice of the z–y plane. (i–l) Intensity profiles along a vertical line through the centre of the PSFs. The x‐axis shows the relative distance to the centre of the image in micrometres. Scale bars represent 3 μm and 0.6 μm in the low‐ and high‐ NA cases, respectively. For images of z–y planes (b, f, d, h), scale bars for the y and z axes are indicated by the horizontal and vertical bars, respectively.
PSFs with phase aberrations generated by the vectorial Cartesian propagator at high‐NA (1.3). (a, d, g) Introduced phase masks at the pupil plane. (b, e, h) Slice of the x–y plane at z=0. (c, f, i) Slice of the z–x plane at y=0. (j–l) PSF with vertical astigmatism at three z‐planes (in focus (k), 500 nm above (j), and 500 nm below focus (l)). All PSFs were generated with circular polarisation, except the Half‐Moon PSF where x‐axis linear polarisation is being used. The scale bars in all images represent 0.6 μm.
In Figure 5, we present additional PSFs computed with the vectorial propagator in the high‐NA setting. The impact of the Gibson–Lanni correction factor on a beam is shown in Figure 5a–c. The refractive indices and thicknesses of the sample, immersion medium, and glass coverslip are ns=1.3, ni=1.5, ng=1.5, ts= 1 μm, tg= 170 μm, while ti is computed using Equation (14). We observe that the spherical aberration it introduces degrades the quality of the focus. We also show the donut PSF, which has a vortex phase mask in the pupil plane (Figure 5d–f) and the half‐moon PSF, which has a π‐phase jump in its pupil plane (Figure 5g–i). Finally, we demonstrate an example with arbitrary phase aberrations using the Zernike polynomials in Figure 5j–l. Here, some amount of astigmatism is introduced, as is often done to encode defocus information,47 and we show the evolution of the beam shape along the z‐axis.
Computational performance
5.2
Speed benchmark
5.2.1
We benchmark the runtime to compute a single 201×201 2D PSF image against the size of the pupil. The image is captured at focus and all input parameters take default values in our library. We compare the runtime of all four propagators on CPU and GPU and show the results in Figure 6. The benchmarking was performed on a machine with an Intel i9‐10900X CPU and an NVIDIA GeForce RTX 3090 GPU.
Time to generate a 201×201 PSF in terms of the numerical size of the pupil and the propagator, on CPU (a) and GPU (b). Each data point in the plots is averaged over 10 runs.
We observe that the runtime of all propagators increases with the size of the pupil on CPU. The Cartesian propagators (red) are faster than the spherical ones (blue) on small sizes (<512 pixels) but slower on larger sizes. This illustrates the computational complexity of each method. Scalar propagators (dotted) are faster than their vectorial counterparts (solid), by roughly 1.5 times for the Cartesian and 3 times for the spherical cases. On GPU, Cartesian propagators are faster than on CPU at the same grid size and the curves (red) behave similarly. Spherical propagators, however, exhibit a flat curve (blue), which indicates that they benefit heavily from the GPU parallelisation. The speed improvement between scalar and vectorial propagators is small, especially for the Cartesian case. Hence, vectorial propagators should be preferred over the scalar ones as the accuracy gain does not come at a high computational cost. Moreover, Cartesian propagators are recommended if one works with images of small size for both scalar and vectorial cases. For larger sizes, however, spherical propagators should be preferred, especially when a GPU is available.
Accuracy benchmark
5.2.2
We benchmark the accuracy of the Cartesian and spherical scalar propagators with the analytic Airy disk for asymptotic limit (the Fourier transform of a perfect circular aperture) given by
where J1 is the Bessel function of the first order of the first kind. As the Airy disk is the Fourier transform of a perfect unit modulus circular aperture, an additional factor sz=cosθ is introduced in Equations (3) and (5). This corresponds to a paraxial approximation.
The spherical integral is computed using two numerical integration rules: Riemann rule (a first‐order method) and Simpson rule (a fourth‐order method). The error δ is the L2‐norm of the difference between the output electric field E of the propagator and FAD, as in
We observe in Figure 7 that the error decreases as the number of points in the integration domain increases in all cases. The spherical propagator has a linear convergence rate using Riemann rule (green) and fourth‐order convergence rate using Simpson rule (blue), which correlates well with their expected accuracy. The Cartesian propagator, based on our custom FFT, shows a convergence rate between first and second order.
Accuracy of the Cartesian (red) and the spherical propagators using Riemann rule (green) and Simpson rule (blue). The step size of integration is h.
CONCLUSION
6
In this work, we have introduced a unifying theory for accurate PSF models, revisiting previous approaches and generalising correction factors. This framework demonstrates the equivalence of Cartesian and spherical methods as different parameterisations of the same propagation integral, providing a simplified understanding of light propagation in high‐NA systems.
We have also developed a high‐performance implementation in PyTorch, which is open‐source and supported by comprehensive documentation. This library allows for efficient CPU/GPU computation of PSFs, with functionalities such as automatic differentiation that make it particularly suitable for optimisation pipelines.
In practice, the choice of propagators depends on the symmetry of the pupil function. Based on our benchmark, for axisymmetric pupil functions, spherical propagators are recommended thanks to their high accuracy and scalability. They are also particularly amenable to GPU parallelisation. For non‐axisymmetric cases, such as those involving Zernike aberrations or specialised phase masks, Cartesian propagators are required. In all these cases, the difference in computational time is relatively modest for pupil sizes up to a few hundred pixels. Hence, the vectorial Cartesian propagator could be a solid default choice in most applications.
We hope for our PyTorch library to contribute to the rapidly growing field of applying deep neural networks on physical imaging models.48 Learning‐based methods have demonstrated their effectiveness in various applications, such as deconvolution49, 50, 51 and 3D SMLM,52 with state‐of‐the‐art deep‐learning tools such as DeepSTORM53 and DECODE.54 For instance, some use cases of our framework include generating large reference datasets to train networks or adapting images based on physical rules to enable learning in self‐supervised inversions55 and generative adversarial networks.56
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Matthews, C. , & Cordelieres, F. P. (2010). Metrolo J: An Image J plugin to help monitor microscopes' health. In Image J User & Developer Conference Proceedings (pp. 1–6). https://imagej.net/events/conferences
- 2Theer, P. , Mongis, C. , & Knop, M. (2014). PS Fj: Know your fluorescence microscope. Nature Methods, 11(10), 981–982.25264772 10.1038/nmeth.3102 · doi ↗ · pubmed ↗
- 3Miora, R. H. D. , Rohwer, E. , Kielhorn, M. , Sheppard, C. , Bosman, G. , & Heintzmann, R. (2024). Calculating point spread functions: Methods, pitfalls, and solutions. Optics Express, 32(16), 27278–27302.39538569 10.1364/OE.523532 · doi ↗ · pubmed ↗
- 4Nelson, G. , Alexopoulus, I. , Azevedo, M. , Barachati, F. , Belavaev, Y. , Carvalho, M. T. , Cesbron, Y. , Dauphin, A. , Corbett, A. D. , & Dobbie, I. M. (2022). Monitoring the point spread function for quality control of confocal microscopes. 10.17504/protocols.io.bp 2l 61ww 1vqe/v 1 · doi ↗
- 5Lelek, M. , Gyparaki, M. T. , Beliu, G. , Schueder, F. , Griffié, J. , Manley, S. , Jungmann, R. , Sauer, M. , Lakadamyali, M. , & Zimmer, C. (2021). Single‐molecule localization microscopy. Nature reviews Methods Primers, 1(1), 39.10.1038/s 43586-021-00038-x PMC 916041435663461 · doi ↗ · pubmed ↗
- 6Sibarita, J.‐B. (2005). Deconvolution microscopy. Microscopy Techniques, 95, 201–243.10.1007/b 10221516080270 · doi ↗ · pubmed ↗
- 7Sage, D. , Donati, L. , Soulez, F. , Fortun, D. , Schmit, G. , Seitz, A. , Guiet, R. , Vonesch, C. , & Unser, M. (2017). Deconvolution Lab 2: An open‐source software for deconvolution microscopy. Methods, 115, 28–41.28057586 10.1016/j.ymeth.2016.12.015 · doi ↗ · pubmed ↗
- 8Heintzmann, R. , & Cremer, C. G. (1999). Laterally modulated excitation microscopy: Improvement of resolution by using a diffraction grating. In Optical Biopsies and Microscopic Techniques III (Vol. 3568, pp. 185–196). SPIE.