Calculating Point Spread Functions: Methods, Pitfalls and Solutions

The knowledge of the exact structure of the optical system PSF enables a high-quality image reconstruction in fluorescence microscopy. Accurate PSF models account for the vector nature of light and the phase and amplitude modifications. Most existing real-space-based PSF models fall into a sampling pitfall near the centre position, yielding to the violation of the energy conservation. In this work, we present novel, to the best of our knowledge, Fourier-based techniques for computing vector PSF and compare them to the state-of-the-art. Our methods are shown to satisfy the physical condition of the imaging process. They are reproducible, computationally efficient, and easy to implement and easy to modify to represent various imaging modalities.


Introduction
Many studies have been conducted to represent, model and approximate the impulse response of a fluorescent microscope, the point spread function (PSF) [1][2][3][4][5][6][7][8].An accurate PSF model is an important requirement for a successful image reconstruction [9].A PSF model can be based on scalar diffraction theory or account for the vectorial nature of the electromagnetic field, such as the polarization state of the field.The former generally outperforms the latter in terms of computation time, but its validity is limited to the case of low numerical aperture (NA) optical systems, i.e. a small maximum half opening angle,  max , which corresponds to sin( max ) ≲ 0.4 [10].A very popular scalar PSF model was developed by Gibson and Lanni (GL) [1,11].Li et al. implemented this scalar PSF model computationally by approximating Kirchhoff's integral using a Bessel series [2].The method is fast as it assumes radial symmetry in the PSF, however this assumption limits the generality of the type of optical aberrations that can be included in the model.
In high NA systems, sin( max ) > 0.7 [10], the vector nature of light is required to produce a realistic PSF [3,4,12].For this it is common to use the Richards and Wolf (RW) model.The RW model mathematically formalizes the most accurate approximation of the exact representation of the vector diffraction pattern in optical systems [3,4].The three integrals in its formulation are however expensive to compute.The model was modified by Török and Varga (TV) to include planar interfaces and refractive index mismatches [13].The TV model was later generalized to account for a stratified medium [8].Additionally, Haeberlé combined the TV formulation with the GL model to compute the PSF of a conventional microscope [14].Ghosh and Preza adapted the method developed by Haeberlé to account for aberrations due to spatially variant refractive indices in the sample and validated the model experimentally [15].The consideration assumes patches (tiles) of the sample for which the assumption of shift-invariant imaging holds.The RW model was also computed by Kirshner et al. [5].A toolbox that allows free access to this model is freely available online under the name PSFGenerator [16,17].This toolbox has been widely used for deconvolution [9,18].However, based on our experience, the software does not support large window sizes and the computation is not time-and memory efficient.Compared to the other PSF models, we also observed that the violation of the missing cone near the focus is more pronounced in PSFs computed using this toolbox.Aguet et al. also presented a method for PSF computation [6].They use a numerical integration based on Simpson's rule to compute the scalar and vector PSF.A recent work in PSF modeling consists of computing the PSF by combining several weighted Gaussian kernels shifted to different positions [7].The efficiency and accuracy of this model critically depends on the number of single Gaussian functions used.Lastly, a free and open-source Python software package for PSF calculation, PyFocus, was also released [19].This software calculates the vector field beyond the paraxial approximation using a custom-made 2D trapezoidal integration algorithm.
Despite the numerous advantages of many of the PSF models described above, each model still has its limitations.Many models exploit radial symmetry to achieve fast computation.Yet, such models cannot directly support asymmetric aberration that may be present in the system.Non-radially symmetric aberrations often include coma and astigmatism.In this work, we present four Fourier-based, comprehensive, easy to implement and more realistic methods for computing PSFs, two of which are completely new to the literature to the best of our knowledge.The use of Fast Fourier Transform (FFT) on a regular grid allows us to efficiently include non-radially symmetric effects.Nevertheless, the FFTs have their drawbacks (e.g.border wrap-around, a form of aliasing) which must be carefully accounted for.We therefore present the pitfalls that one may encounter in the calculation and ways around them.Our models also account for losses by reflection on non-coated optical surfaces (e.g. the coverslip) via Fresnel coefficients.Additionally, the Fourier plane can be easily accessed and sample-induced aberrations that may arise during imaging can also be introduced into the PSF calculation.The models are demonstrated to be fast compared to state-of-the-art models supporting similar features.They can easily be modified to represent different imaging modalities or to include a tilt of optical surfaces in the imaging system.Common situations such as a tilt of the coverslip can thus be modelled.A further aim of this manuscript is to release a toolbox to the scientific community, which others can benefit from for calculating PSFs using uniform or modified (aberrant) apertures [20].The software code is written with MATLAB programming language and is released in an open source format.It can easily be converted to any other programming language such as Python.
Before we introduce the four Fourier-based models in Section 3, we describe the necessary tools that are needed for their computation in Section 2. We then conduct a quantitative comparison of the models with the state-of-the-art in Section 4. We conclude our results in Section 5.

Abbe sine condition
An optical imaging system usually fulfills the Abbe sine condition to image a plane in the object to the detector plane (e.g. a CCD or CMOS camera) [10,21].The Abbe sine condition relates any angle of a beam emitted by the sample ( em ) to its corresponding angle reaching the detector ( det ).It is given by sin( em )/sin( det ) = , with  being the magnification of the optical system in terms of geometric optics [10,22].As shown in the thin lens approximation in Fig. 1, a simple 4  -imaging microscope imaging system with a magnification of  = 2 would not fulfill the above-mentioned Abbe sine condition.
The concept of Gaussian reference sphere, defined as the collection of equivalent refractive loci for an aplanatic imaging system, is introduced to maintain the simplicity of such drawings (see Fig. 2).This sphere is centred at the intersection point of the object plane or the image plane Geometrical representation of a 4f-imaging system using the thin lens approximation.Note that here tan( em )/tan( det ) =  2 /  1 , which violates Abbe's sine condition.BFP stands for back focal plane.Here lens 1 corresponds to the microscope objective lens while lens 2 represents the tube lens.
with the optical axis.Its radius is equal to the focal length  1 or  2 of the respective imaging lenses [22].The rays propagate to the Gaussian reference sphere and get projected to a plane normal to the optical axis without acquiring any extra phase.This projection is indicated by the red dashed lines in Fig. 2. The same principle is applied in reverse to the tube lens (lens 2 in Fig. 2), where the effect can often be neglected due to the usually large imaging magnification and thus small angles near the image plane (see right side of Fig. 2).
em  det lens 1 lens 2 ℎ Gaussian reference sphere Gaussian reference sphere ''projection'' Fig. 2. Representation of the Gaussian reference sphere, forcing sin( em )/sin( det ) =  2 /  1 to be constant in agreement with Abbe's sine condition."Teleportation" means that the beams are continued at the connected surface without acquiring any phase for the space in between the plane normal to the optical axis and the Gaussian reference sphere.

Energy conservation
The energy of light propagating through the aplanatic focusing system must be conserved.
Integrating over the radial position in the back focal plane (see Fig. 1 and 2) must yield to the same quantity of energy as integrating over the corresponding angle  in the equivalent refractive locus.Given the projection of the field in the Gaussian reference sphere onto the normal planes, the aplanatic factor (AF) needs to be accounted for to conserve the energy.
To understand and derive the AF, let us assume an isotropic emitter placed at the centre of the Gaussian reference sphere,  0 (see Fig. 3).Since the emitter is emitting uniformly in all directions, the strength of the amplitudes on the Gaussian reference sphere is uniform.Let us also consider a small parallel ray which is redirected from the Gaussian reference sphere at a given incidence angle .If we attribute a given irradiance  0 as a power  0 per unit area  0 to such a beamlet from  0 , the same power will have to be distributed to a smaller area  1 =  0 cos() after the projection of the beam from the Gaussian reference sphere to the plane parallel to the pupil plane in the BFP (see Fig. 3).This means that the irradiance measured perpendicular to the local direction of propagation changes to  0 /cos(), and each of the field vector component for this beamlet thus needs to change by 1/ √︁ cos() to be consistent with this intensity change.
Fig. 3. Schematic diagram illustrating the "aplanatism effect".The two mutually parallel planes Ξ 0 and Ξ 1 are perpendicular to the optical axis.The area  0 is a projection of the elementary area  1 of the plane Ξ 1 onto the Gaussian reference sphere.
If a large enough magnification of the imaging system is assumed, the vectorial and aplanatic factor effects caused by the tube lens can be neglected.
To understand the factor that needs to be applied to conserve the energy when focussing a uniformly illuminated 2D pupil (BFP) with a high-NA objective (excitation PSF), let us consider the same Fig. 3 but with the light travelling from the right side to the left.The point  1 denotes a point in the pupil plane from which a beamlet emerges.The parallel beamlet will focus onto  0 after the objective lens.We can consider the Helmholtz reciprocity theorem: a lossless (non-magnetic) monochromatic optical system in which a field (  0 ,   0 ,   0 ) (here an isotropic emitter) at object plane position  0 gives rise to a field (  1 ,   1 ,   1 ) at a point in the image plane via a virtual aperture around  1 , warrants that placing the isotropic emitter (  0 ,   0 ,   0 ) as a source at that same point in the image plane, generates the field (  1 ,   1 ,   1 ) at the object plane position  0 [12] with the same virtual aperture at  1 .A uniform emitter near the focus of the tube lens, leads, due to the low NA of the tube lens to a uniform illumination of the pupil plane (BFP), corresponding to our assumption for calculating the excitation PSF.Therefore the excitation PSF should have the same aplanatic correction as the emission PSF as long as the magnification is large so that we can neglect the aplanatic factor of the tube lens.This is confirmed also by considering the √ cos  dependence on the Gaussian reference sphere and thus on the McCutchen pupil (detailed in the next section).To arrive at the 2D Fourier-transform of the in-focus excitation field, we need to project the McCutchen pupil along   and thus apply a projection factor of 1/ which leads to an overall amplitude of √ cos /cos  = 1/ √ cos  confirming the above reciprocity argument.

The McCutchen pupil
To calculate the PSF using a Fourier optics, let us consider the diagram in Fig. 4 depicting a Huygens wavelet contributing to the focusing of a monochromatic coherent plane wave by a high-NA microscopy objective.The beam enters the objective system from the right side and is spatially limited by the entrance pupil of the optical system.An "aperture stop", is in practice either intentionally introduced to warrant the linear shift invariant performance of the system and to avoid aberrations from unwanted beams, or effectively provided by the inner geometry of the objective.The limitation of the beams is therefore approximated to be at the limit of that aperture stop.At this pupil plane, every point   on the wavefront is considered as a source of a Huygens wavelet, denoted by W [23].Each such Huygens wavelet is approximated by a plane wave in object space, as seen by the phase fronts within the Gaussian reference sphere.For the process of emission, which we now consider, we can interpret Fig. 4 by decomposing the emitted wave into the same plane waves in object space, corresponding to pupil-plane wavelets.At first, we limit ourselves to a scalar electric field, where the field is described by a single amplitude value as a function of spatial position.The vector nature of the electric field will be considered further down below.The aperture plane can be seen as a coherent superposition of spherical wavelets, each of which gives rise to a plane wave in object space directed towards the nominal focus point S (see the wavelet labelled "W" in Fig. 4).According to the Huygens-Fresnel principle, the spherically converging wave is obtained by superimposing all these plane waves [23].These superimposed wavelets have to acquire exactly the same optical path length and constructively interfere at S. In other words, the phase at the nominal focus is identical for all such wavelets and can thus be set to zero in our simulation.For convenience, we choose S as the centre of our real-space coordinate system.The wavelet W giving rise to a plane wave in focus (Fig. 4) can now conveniently be described in Fourier space as a single point P W (Fig. 4), i.e. a single 3-dimensional vector ( ì   ) in Fourier space.All such vectors necessarily reside on a sphere of a radius  0 = 2/ em ,  em = /,  em and  being the vacuum wavelength corresponding to the emitted wave respectively, and  the refractive index of the embedding medium.
A pupil position in real space corresponds to the lateral components ì  , of the wave-vector ì .This linear correspondence ì  , =  ì  , / 0 with the focal length  of the objective is forced by the Abbe sine condition between the pupil plane coordinate and the -vector position of the wave near the focal plane.
The pupil plane aperture stop thus gives rise to a three-dimensional cap residing on the -sphere in Fourier space (see a section though the cap as the bold curve in Fig. 5 in the Fourier space representation).As the aperture is limited by the NA of the objective, the 3D frequency spectrum in Fourier space is represented in a segment of the -sphere sphere.This segment is called "generalized aperture" or "McCutchen pupil" [24].
According to McCutchen, the (complex-valued) amplitude distribution of the electromagnetic field in real space near the focus S corresponds to the three-dimensional Fourier transformation of the amplitude distribution on the McCutchen pupil [24].To compute the intensity PSF, we can thus project the amplitude on the McCutchen pupil and perform an inverse three-dimensional Fourier transformation.on this basic understanding.The four different methods differ in how such a distribution is obtained and how the field is propagated in the homogeneous medium.

Sampling condition of the computation
A PSF calculation typically samples the continuous mathematical function at discrete locations (delta-shaped points).The imaging of a point emitter, which is practically detected on a pixelated imaging device such as a CCD or CMOS camera.Those devices integrate, in each of their rectilinearly spaced pixels, over the continuous signal weighted by a (box-shaped) pixel sensitivity function.The local integration of the PSF in every pixel by the camera can be rewritten as first convolving the PSF with the pixel sensitivity function and then sampling it at regularly spaced positions.Due to the convolution theorem, the effect of detector integration can be represented by a multiplication of the Fourier transform of the PSF, the optical transfer function (OTF), with the Fourier transform of the pixel sensitivity function.If we assume box-pixels with uniform sensitivity, the OTF gets modified by a multiplication with a sinc   / samp sinc   / samp .This means that at the current sampling frequency  samp = 2/ samp ,  samp being the pixel pitch, the overall transfer would cross zero.To avoid aliasing, the sampling of the PSF  samp must satisfy the Nyquist Shannon theorem given by:  samp < 2/(2 limit ), with the limit frequency  limit referring to the transferred frequency limit, to avoid potential aliasing of information within the frequency band [25].
In a wide-field system, the maximal in-plane spatial frequency is derived from Abbe diffraction limit and is given by   ,max = 4NA/ em [25].Similarly, the maximal axial frequency in real space for a wide-field microscope is given by  ,max = (2(1 − cos()))/ em .The highest frequency of the intensity result has to be sampled with at least two positions per shortest period that can be transmitted by the system [25].This requires the pupil to fit into half the digital Fourier-space representation such that its autocorrelation (i.e. the incoherent OTF) fits in the digital Fourier space.The maximal pupil radius in Fourier space along  or  should be lower than half the maximally represented frequency along   or   in our Fourier-space representation.

Jinc aperture
In addition to the potential error with the aforementioned sampling, a digitization of the usually round pupil in Fourier space as a hard aperture onto a rectilinear grid may also induce severe artefacts that we need to avoid when using the FFT operator.To illustrate this problem, let us consider a field distribution of a high-NA PSF with numerical aperture NA = 1.4,refractive index  = 1.518 and emission wavelength  em = 580 nm.We calculate the pupil radius, which also corresponds to the Nyquist frequency, using the theory stated in Section 2.4.We denote  max this pupil radius.We generate a hard aperture with radius equal to  max /8 and calculate the corresponding field distribution in real space by generating the Fourier transform of the hard aperture.The window size for this first experiment is 128 × 128 pixels.
For symmetry reasons, we would expect a perfect circularly symmetric PSF.A significant deviation from circular symmetry is however clearly visible in Fig. 6a.By repeating the same calculation for 1024 × 1024 pixels, the discrepancy is significantly reduced even though it is still not totally spherically symmetric (see Fig. 6c).However, calculating on such large grids causes about significant computational time overhead (0.07 s vs 0.025 s i.e. about 3 times slower).This computation might also be unnecessary as we may not need so many pixels of the PSF far away from its centre.
Interpolation in Fourier-space using a sinc function to obtain a better representation of the pupil may be one route to re-establish spherical symmetry.However, this is a tricky business [26].Here we present a different approach.We calculate the two-dimensional (2D) Fourier-transform of the uniform pupil disk, for which the analytical solution in real space is well known: jinc() =  1 ()/,  1 being the Bessel function of the first kind and  the radial coordinate.We therefore obtain an "ideally" representation of a disk in Fourier-space by Fourier-transforming a two-dimensional jinc function.This "interpolated" disk is then appropriately modified with -space dependent phase and magnitude alterations.
The computation time of the 1024 × 1024 pixels ideal representation of disk in Fourier-space using the jinc trick is 0.3 s on average (Windows 10, 64-bit, Intel ® Core TM i5-3570S CPU @ 3.10GHz, 8,0GB RAM, Intel HD Graphics).
As the jinc-function exhibits first order discontinuities in real space at the border, causing unwanted high-frequencies [27] in its Fourier-transform.To reduce this effect, we modify the jinc-function at the outer rim by appropriately smoothing the 15 % of its edges towards zero ("DampEdge" function in the PSFToolbox [20]).As seen in Fig. 6b, the real-space representation of the field distribution is perfectly symmetric and spherical by design even for images with relatively few pixels.The circular edge of this field in Fig. 6b has been dimmed down using the DampEdge function.We refer to this method of generating an interpolated disc in Fourier space as the jinc-FT trick.

Fourier wrap-around
The grid on which the field propagation is simulated is finite.There exists an axial depth position at which the disk of defocus no longer stays well within the available lateral space provided by the real-space grid.Upon propagation, waves leave on one side of the lateral sampling grid and enter on the opposite side due to the periodic boundary conditions of the Fourier-transform.This causes severe standing-wave effects called Fourier wrap-around (see Fig. 7a and 7d).
-10 0 10 x/µm Three possible strategies can help to avoid this wrap-around effect by: A. Zero-padding the in-focus plane to extend the axial depth from whereon the standing wave patterns occur (see Fig. 7b and 7e), B. Establishing an ideal absorptive boundary condition to the outside boundary and continuing the propagation by re-projecting the filtered field onto the pupil plane [28], C. Using a chirp Z-transform (CZT) to calculate only a part of the propagated field while sapling the pupil at the finest spacing (see Fig. 7c and 7f).
Option A is computationally expensive since twice the initial image window size slows the calculation down by a factor ∼5.However, padding with zero to twice the original size may still yield unacceptable artifacts in 3D PSF calculations.Option B has the disadvantage of sacrificing a good PSF for a portion of pixels near the edge of the lateral grid of the calculation.In addition, two FFTs are required for this approach instead of only one for each propagation step.With option C, wrap-around artefacts can partially be avoided at the expense of more than twice increase of computation time.For our work, we selected option C.This can be observed in Fig. 7c and 7f where the standing waves are seen to be considerably reduced.The use of CZT for Fourier optics and PSF modelling is not new in the literature and has been shown to be more efficient than FFTs with zero-padding without loss of accuracy [29][30][31].

Fourier-based methods for PSF calculation
We present four Fourier-based methods for computing the PSF of a conventional fluorescence microscope in this section.The four methods are based on the finding of McCutchen [24] and differ how the field is propagated in free space.

The electric field on the 𝑘−sphere
To calculate the electric field on the −sphere, we generalize the theory described in Section 2.3 from scalar to vector formulation, accounting for the polarization state of the input field.We associate each plane wave arriving at the focus with the refractive effect ("bending") applied to the corresponding "ray" at the equivalent refractive locus on the Gaussian reference sphere.We assume a perfect anti-reflection coated objective lens (system) transmitting all the energy of rays.We assume a system satisfying the Abbe sine condition, which requires the beams to change direction at the Gaussian reference sphere.Knowing that the electrical field of the plane wave is a transversal wave and neglecting any axial components of the field, we thus need to project the electric field at the 2D pupil plane (  ,   ) onto a plane perpendicular to the new propagation direction, to obtain the electric field (  ,   ,   ) of each position on the 3D McCutchen pupil.The diagram in Fig. 8

Let ì
= (  ,   , 0) denote the incident polarized wave from infinity at the left side of Fig. 8.It is useful to consider a locally varying coordinate system along azimuthal ( ì   ) and radial ( ì   ) directions respectively.Let ì   denote the field amplitude transmitted along the wave vector ì  towards a point (, , ) near the focus where the field is evaluated (see Fig. 8).The unit vector corresponding to the radial component ì   is refracted by  and becomes ì   while the azimuthal component oriented along ì   remains unchanged.
The field amplitude distribution on the −sphere oriented towards a point (, , ) in the image plane is therefore given by: The incident wave field ì   is with a given polarization state.In the presence of refractive index mismatch in the system, a phase aberration, Φ, is introduced into the system and the  and −polarized light are transmitted at a rate determined by the transmission coefficient of each polarization state.The amplitude field on the McCutchen pupil is generalized as follows: 0 = 2/ em being the wavenumber and   and   are the transmission coefficients for  and  polarization respectively.With the assumption of an anti-reflection coated objective lens, we can further assume that the reflection coefficients are smaller than one.Omitting the phase terms in the complex transmission coefficients and add it as part of the phase term in Φ, the expressions of   and   for  number of layers in a stratified medium and  − 1 of interfaces are therefore formulated as follows: with . Any additional phase or amplitude modulation can easily be introduced in the formulation of the electric field in Eq. ( 3).In the next sections, we present the different methods from this work to propagate the fields in free space.

Slice propagation method with FFT (SP-FFT)
This method uses the well known angular spectrum method to propagate the field in free space slice by slice without zero-padding of the input field.Based on the Fourier slice theorem, a −slice in real space corresponds to an integral of the amplitude over the axial spatial frequency   in Fourier space.A defocus by a distance  in real space corresponds to a phase shift equivalent to    in Fourier space.Our method here accounts for the jinc-FT trick described in Section 2.5.1 to avoid the pitfalls of the digital Fourier transform and the aplanatic correction for energy conservation.The steps for its computation is summarized in Algorithm 1: In the presence of stratified medium of different refractive indices, the effective numerical aperture is equal to NA eff = min{NA,   , ∀  ∈ [1, ]},  being the number of layers,   the refractive index of the  th layer and NA is the numerical aperture of the objective lens as given by the manufacturer.The pixel pitch in , ,  is stored in the variables [   ,   ,   ] and are in the same units as the emission wavelength  em .We denote [  ,   ,   ] the size of the grid in , ,  respectively in pixels.

Slice propagation method with CZT (SP-CZT)
The SP-FFT still has its limitations due to Fourier wrap-around at higher depth as described in Section 2.5.2.The CZT is an alternative operator to reduce or totally avoid the wrap-around effect since it allows to specify the (zoomed) region after the transform implicitly applying zero-padding to the field to transform.For a 1D signal   ,  ∈ [0,  − 1] ∩ N with  being the number of points of the signal and N the set of natural numbers, the Z transform X ,  ∈ C is given as follows: where   =  − ,  ∈ N is a spiral path (−path) in the complex plane with  being the starting point and  = exp(−Δ) the ratio of two consecutive points with a given angular increment phase Δ.For  = 1 and Δ = 2/,   is computed over an unit circle and the CZT operation becomes a discrete FFT.To zoom the signal   in by a scalar factor ,  = exp(−/) and  = exp(−2/()) [33].Therefore, Eq. 5 can be expressed in terms of convolution as follows [29]: The inverse CZT of a signal X  in a frequency-domain representation is defined as the complex conjugate of the CZT of the complex conjugate X *   of X  within some scaling factor for a CZT operating on a unit circle [34].
In this approach, we calculate a zoom factor  to perfectly fit the pupil near the edge of the array to transform, corresponding to Nyquist sampling of the amplitude in real space.By doing so, we use all the available lateral space in Frequency space.We calculate the field in Frequency space according to Step 1 to 5 in Algorithm 1.If necessary, we can calculate the field on the pupil plane on a finer grid to obtain an even finer sampling in Fourier space.We use the inverse CZT to zoom-in (i.e.zoom out in the pupil plane) using the zoom factor  to obtain the corresponding real space field with the desired sampling.To minimize wrap-around, the appropriate zoom factor  is calculated such that the lateral window size of the calculated PSF is slightly bigger or equal to the lateral dimension of the PSF at the highest depth position Δ from the focus.The maximum angular aperture of the optical emission determines the maximum angle at which an oblique emission will propagate.We denote D the lateral radius of the propagated beam at Δ (see Fig. 9).
The factor  is therefore calculated as  = ( +    /2)/(   /2), where tan  max = /Δ in real space and tan  max =    /  in Fourier space with  max being the maximal angular aperture and    = [  ,   ] the number of pixels in the in-focus -plane (see Fig. 9).If the number of pixels in  and  directions are not equal, the zoom-in factor  will also be different in both directions or can be set to be the minimal value.Given that the in-focus plane is calculated over    pixels, the PSF may require  ′   = 2(   /2 + ) lateral grid.We modify Algorithm 1 to introduce the wave propagation using CZT method and the result is summarized in Algorithm 2. This method aims at representing the useful part of the McCutchen pupil directly in 3D-Fourierspace and projecting the two-dimensional pupil modifications (the aplanatic factor and vectorial projection effects) onto this three-dimensional shell.The difficulty is that the shell, at each integer [  ,   ] position, has a non-integer   position leading to an interpolated representation of   in Fourier space.As a credible representation of such a non-integer   would require essentially the full   -range.To keep the computation efficient, an appropriate compromise was therefore made.We calculate a table of interpolation kernels ("interpolators"), each of them applicable for a specific (small range of) sub-pixel shifts.We aim to represent only the central part of the corresponding real-space representation as faithfully as possible and label a border region as "don't care" region (see Fig. 10b).This "don't care" region is limited by a chosen factor  reg (here it is chosen to include 8 pixels from both edges).The part of real space within this border region is iteratively updated, while the central part is forced to the expected values in each iteration in this iterative Fourier transformation algorithm (IFTA) to obtain the interpolation kernels.
In addition, a pre-defined cut-off frequency   cut-off is chosen.This cut-off frequency limits the number of interpolation coefficients given by   = 2  cut-off + 1, which can be used to fill the voxels along   in Fourier-space adjacent to the one nearest to the non-integer   (  ,   ) position of the McCutchen pupil.The required interpolation coefficients are generated with the help of the IFTA [35].We denote  subpix the number of sub-pixel positions along   .As initialization, ideal exp(2  ) waves are generated in real space corresponding to the respective sub-pixel frequencies in Fourier space.The ideal waves are then Fourier-transformed and only   interpolator values are kept while all others are set to zero.The result is transformed back to real space, where the central area is replaced by the original perfect waves, but the "don't care" region is not touched.This is repeated over  iter iterations (typically 500 times) until convergence.The so-generated interpolation table of size   ×  subpix is stored for later use (see Fig. 10c).In the example presented here,   cut-off = 8 leading to   = 17 interpolation coefficients to be determined.An interpolation table of  subpix = 60 sub-pixel positions of the 17 complex valued coefficients is pre-computed via IFTA.The inner part in the real space regime which represents the "do care region" is about 66 % of the given -range.A typical example for the offset of 0.25 pixels is shown in real and Fourier space in Fig. 10b and 10d respectively, overlayed with the ideal subpixel wave (solid line which corresponds to the legend 'Original signal').The border of the "don't care region" is indicated by the dashed vertical lines.A real space representation of the full interpolation table is shown in Fig. 10a with the "don't care region" also indicated by the vertical red dashed lines.
The size of the border factor  reg (in pixels) and the cut-off frequency   cut-off defining the number of interpolation coefficients   should be roughly the same.If the "don't care region" is far bigger than the "do care region", fewer non-zero interpolation coefficients are needed in Fourier space, speeding up the algorithm in the application phase.However, the integrated intensity of the PSF is less uniform, increasing the violation of the missing cone property, which will be further discussed in Section 4. A small "don't care region" on the other hand can lead to inaccuracies inside the "do care region" if the number of interpolation coefficients is kept small.The principle for the PSF generation using this method is summarized in Algorithm 3.
This PSF calculation method (step 5) can be performed fast and memory efficient as a single access operation in MATLAB by exploiting its indexed addressing capabilities.In this way, the complex-valued 2D pupil can be rapidly written into the appropriate 3D Fourier space region with the optimized interpolation coefficients as described above and the "don't care" region can be later removed by cropping.The required   -range can be kept to a minimum.

Algorithm 3 Fourier shell interpolation (F-Shell)
Input: NA eff ,  em , , [   ,   ,   ], [  ,   ,   ],  reg ,   cut-off ,  subpix ,  iter Output: ℎ: PSF intensity, ì ℎ   : complex amplitude field 1: Generate the McCutchen pupil projections ì  ′  (  ,   ) as described in Algorithm 1, Step 1 to 5, using the jinc-FT trick as an aperture delimiter 2: Calculate   (  ,   ) for every pixels within the pupil and round it to the nearest 1/ subpix subpixel   position 3: Generate a 3D index of the full −sphere 4: Store the interpolation coefficients for future use 5: For each PSF to calculate: Write the field ì  ′  calculated in Step 1 with the appropriate interpolation kernel for the sub-pixel at the 3D index position of the −sphere 6: Perform a 3D FFT of the result from Step 5 to obtain the three-dimensional field distributions ì ℎ   (with expected errors in the "don't care region") 7: Calculate the PSF intensity: This method is derived using the knowledge that the three-dimensional Fourier transform of a complete spherical shell is a sinc( 0 | |) function,  the spatial radial coordinate.To compute the −sphere representation, in analogy to the jinc-trick, we apply the three-dimensional FFT to a sinc( 0 | |) calculated in real space.The sinc-shell method is described in Algorithm 4.

Algorithm 4 Sinc-R method
Output: [ℎ : PSF intensity; ì ℎ   : complex amplitude field] 1: Set the calculation grid to be 1,25 times bigger in both direction  and  than the desired grid 2: Generate a sinc( 0 | |),  being the radial coordinate, amplitude distribution in three dimensions in real space over the region described in Step 1 3: Perform an appropriate DampEdge to the edge region of the image (e.g. 5 % on each side of the image border) 4: Compute the 3D FFT of the distribution in Step 3 5: Keep only the (positive frequency)   -range, which contains valid (non-zero amplitude ) ì  vectors.This cropping in Fourier space changes the -sampling and causes a phase ramp in real space.Note that the latter does not affect the intensity values.For a given sampling of the final intensity PSF, these steps can be adjusted accordingly.This method has the attractive property that it does not suffer from the Fourier wrap-around effect since the field was directly generated in real-space.A disadvantage is each step has to be performed observing Nyquist sampling along   for the full field including its -propagation.This method is also not readily applicable to a single slice (in or out-of-focus).

Comparison methods
We compare the four models described in this project with state-of-the-art (see below) PSF models.The Richards and Wolf (RW) model was chosen as a gold standard for the comparison because it is a well accepted model of the field in an aplanatic system and its accuracy relies on the computation of the three integrals in its analytical expression.We denote the gold standard by RW-GS.To avoid interpolation problems in the calculated PSF especially near the focus position, we sample in our RW-GS 10× higher than the standard sampling (see below) in  and .This oversampled RW-GS is then subsampled by considering only every 10 th pixel along  and  to correspond to what we choose as "standard sampling".The standard sampling corresponds to a voxel size of 83 × 83 × 100 nm 3 and a calculation grid of 127 × 127 × 65 pixels.We only consider 90% of the calculation grid, the central 115 × 115 pixels in each slice in our quantitative comparison to discount any artefacts that may be present at the edge of the border of the calculation grid.We assume a water objective lens of numerical aperture equal to 1.2 and an emission wavelength of 510 nm in the simulation.
We choose four state-of-the-art commonly used PSF models to compare our methods to.The first model is the scalar PSF based on the work of Gibson and Lanni [1], and further developed by Li et al. [2].This technique calculates the PSF fast by using a combination of Bessel series.We denote this model GL in our comparison.The second and third models are a scalar and vector PSF as described in [6], computed using a numerical integration based on Simpson's rule.We denote sPSF the scalar PSF and vPSF the vector PSF.The GL, sPSF and vPSF firstly compute a -map,  and  being the radial and axial coordinate, of the PSF.The radial symmetry is used to fill the whole volume linearly interpolating in the −map.The last state-of-the-art PSF to compare to is the vector PSF calculated with the Richards and Wolf 3D optical model in the PSFGenerator toolbox in [5].We denote this PSF by PSFGen.The four models described in this manuscript are distinguished from those state-of-the-art using an asterisk: SP-FFT*, SP-CZT*, F-Shell* and Sinc-R*.GL, sPSF, vPSF and PSFGen are not centred at the same centre position as our PSFs models for even sizes if the calculation window.We therefore choose odd-sized calculation window grids throughout this document.SP-CZT* is however calculated on even grids and is further cropped to obtain the same window size as the other PSF models since the CZT function does not center the PSF at the same position as the other models for odd-sized grids.
In Section 4.2, we investigate and quantify the accuracy of the intensity values of the PSFs.We use as error metric the relative square error (RSE) of a particular calculation with the gold standard RW-GS to quantify the performance-error of each simulation in comparison to the RW-GS with PSF sim being the simulated PSF to compare with the RW-GS PSF.We further evaluate to what extent each PSF model satisfies the conservation of energy by quantifying the violation of the missing cone in Section 4.3.
As the models described in this work are Fourier-based, the effect of any remaining FFT wrap-around effects needs to be quantified to conclude on their accuracy.To achieve this for each PSF model, we calculate the error  between a reference PSF, ℎ ref computed in a window demonstrated to contain least standing waves due to Fourier wrap-around, and the PSF calculated with the same method in a different calculation grid .Details on this approach are given in Section 4.4.The error  is calculated using Eq. 8.It quantifies the amount of standing waves due to Fourier wrap-around in a given simulated PSF.

𝜖 = ∑︁
,, Finally, we present the computation cost of each method in Section 4.5.[5], SP-FFT* to the method described in Algorithm 1, SP-CZT* to the method employing Algorithm 2, F-Shell* to the Fourier shell interpolation method in Algorithm 3, GL to the scalar PSF based on the Gibson and Lanni model [2], sPSF and vPSF to the scalar and vector PSF from [6].

Intensity profiles of the PSFs
To observe finer differences between our chosen gold standard RW-GS and the intensity profiles of the PSFs of various methods, we display in Fig. 12 the radial mean intensity in log scale at the focus position.We observe in Fig. 12c a clearly visible discrepancy between the scalar PSF models (GL and sPSF) and vector models (RW-GS and vPSF).PSFs calculated by vPSF and PSFGen fit the profile of the RW-GS at high precision.Sinc-R* PSFs departs from the profile of the RW-GS beyond a distance of about half micron from the centre.The SP-FFT*, SP-CZT*, and F-Shell* methods are slightly inaccurate near the edge of the radial position.
We compute the relative local error map between the simulated PSFs, PSF sim , relative to the RW-GS using the following equation Eq. ( 9) to visualize the error in the volume PSF and to indicate how precise the local prediction is.
This map represents the absolute of the difference between the simulated PSFs and the RW-GS relative to the intensity values of the RW-GS at each pixel of the volume PSFs.The map for each of the PSFs is displayed in Fig. 13.
It is noticed from the error map that the error in the Fourier-based models are more dominant towards the edge of the grid while the error is less significant closer to the optical axis .This Radial Position/µm - Radial Mean Intensity in log scale   error in the Fourier-based models is due to the FFT and CZT operation.A discontinuity is observed in the error of the RW with normal sampling along the two diagonals of the grid.This is due to the mapping of the −map to form the whole volume PSF.The error near the axial axis is also higher than the error at different position of the grid for the RW due to its sampling.A discontinuity is observed in the −cross section of the scalar PSFs GL and sPSF where the angular aperture  is maximal.This discontinuity is not observed in the vector vPSF, which represents a vector formulation of the sPSF.The error in vPSF is more concentrated near the region closer to the axial axis, similar to the case of the RW.
To conduct further investigation of the error in each model, we compute the RSE to the gold standard per slice () along the axial axis.The results are displayed in Fig. 14.We also compute the RSE over the full 3D calculation volume (see RSE 3D in the labels in each panel of Fig. 14).
-5 0 5 Axial Position/µm It is shown that the SP-CZT* PSF has the least RSE relative to the RW-GS with a 3D RSE value of 0.001910 −3 .This is followed by the RW, Sinc-R*, F-Shell* then SP-FFT*.The 2D RSE in those Fourier-based models become higher at higher depth due to possible Fourier wrap-around present in the PSFs.The vector PSFs, vPSF and PSFGEn, have higher RSE values than the Fourier-based models with the RW-GS.The error in those models are more significant near the focus due to the violation of the missing cone.This concept of missing cone is discussed in details in Section 4.3.The highest RSE are in the scalar PSFs and the errors stand at 4.864410 −3 and 3.697110 −3 for the GL and sPSF respectively.These reflect the limitation of the validity of the scalar approximation for a system with high numerical aperture.

Violation of the missing cone
In wide-field microscopy, the frequency spectrum along the   -axis is missing apart from the   = 0 position due to energy conservation in the system.Mostly independent of the focus position, all power emitted into angles of acceptance, as defined by the pupil, reaches the detector (if sufficiently large).This effect corresponds to a missing cone in the optical transfer function (OTF) that prevents the transmission of information about the object within this cone region (yellow cone in Fig. 15a).Out-of-focus light is distributed to different regions but, the power is still conserved as long as it resides somewhere on the detector.The integrated intensity at each  position should therefore remain constant.This should be the case for any wide-field PSF as long as the PSF remains confined well within the calculation grid.Monitoring the integrated intensity for each −plane at each −position to check for a violation of the missing cone, can, especially near the focus, be conclusive to reveal interpolation errors caused by mapping radial results to the rectilinear grid.In the calculation of the vector field in a finite grid, interpolation errors caused by insufficient data points for the computation can lead to the violation of the missing cone in the corresponding transfer function of the system.A more precise observation of the violation of the missing cone can be made by zooming on the −range around the focus and displaying the integrated intensity of each slice along the axial position  (see Fig. 16).To quantify this effect, the standard deviation of the integrated intensities within 20 % around the focus of the given −range (Std) is calculated and is plotted in Fig. 16.This 20 % region is delimited by the two dashed vertical lines in Fig. 15b, 15c, and 15d.The Std measures the non-uniformity of the laterally integrated intensity over the axial position range.
The Std , which corresponds to the SP-FFT* and the F-Shell* methods are very close to the Std of the highly sampled RW-GS while the SP-CZT* demonstrates the lowest value of Std.The methods described in this work, SP-FFT*, SP-CZT*, F-Shell* and Sinc-R* all have relatively low Std compared to the state-of-the-art PSFs: GL, sPSF, vPSF, and PSFGen.PSFGen is in the first position in violating the missing cone condition followed by RW computed with a standard sampling and GL.In 3D deconvolution routines applied to widefield data, preserving the missing cone property is paramount.

FFT wrap-around effect
To quantify the wrap-around effect, we calculate densely sampled PSFs using each technique on a large calculation window, denoted by  ref using the formula in Eq. ( 10).This PSF is denoted by ℎ ref .
with   |max and   |max are the maximal radial and axial frequency components,   |lim being the resolution limit,   the pixel pitch and 1.3 is a heuristic factor.The first factor 2 in the expression of  ref|  is considered to double the half window and the second factor 2 is to sample the frequency space two times higher.The variables  max ,   and   |lim have the same units and  ref|  is expressed in pixels.The same amount of padding is also applied along −axis.Substantial wrap-around effects should be avoided within this large computation grid  ref .
In order to quantify the wrap-around effect of each PSF model, we then compute the PSF ℎ  with the (typically smaller) window  and compare it to ℎ ref calculated with the window  ref .
The wrap-around effect in ℎ  in relation to ℎ ref is calculated by subtracting ℎ  and ℎ ref over the smaller of the two windows (see Eq. ( 8)) to obtain the amount of standing waves expressed in intensity values in the PSFs.No wrap-around effect is expected for  >  ref .However, since the reference window may already show minor wrap-around effects, the error  is expected to converge to a small non-zero constant for  >  ref .
The  ref is calculated to be 337 × 337 in  and .We compute different window size: 63 × 63, 127 × 127, 257 × 257, 511 × 511, and 673 × 673.The −range is 65 pixels (i.e.65 −slices along ).The rest of the parameters are the same as previously described in the introduction of Section 4.1.GL, sPSF, vPSF, PSFGen, which are state-of-the art methods and RW do not employ any FT operators and are therefore not prone to wrap-around problems.We observe from Fig. 17  red dashed lines passing at the RSE of RW.The RSE of the SP-CZT* falls on the line passing horizontally on RW, vPSF and PSFGen for all the different windows.The Sinc-R* presents an error  slightly above the red lines while SP-FFT* and F-Shell* fall at about 0.05 in RSE above the red lines.This means the wrap-around effect in the SP-CZT* has been completely suppressed while there is still some remaining wrap-around problem in the SP-FFT* and F-Shell* methods.The error in the value of  is less than the errors in SP-FFT* and F-Shell* for the Sinc-R*.This is still an error in the Sinc-R* method that should be improved in the future.All the  errors converge to a constant as expected for window bigger than the reference window 327 × 327.

Computational time
The estimation of the speed of the various algorithms in this section was performed using Windows 10, 64-bit, Intel ® Core TM i5-3570S CPU @ 3.10GHz, 8,0GB RAM, Intel HD Graphics.The profiles of the computation time per voxel for each technique at four different window sizes are displayed in Fig. 18.In the RW model, the three integrals of the electric field [4] are evaluated numerically by exploiting a cylindrical coordinate system for the integration.This numerical integration is computationally advantageous since only the field in the centred radial axis versus the axial position  is calculated.RW, GL, sPSF and vPSF exploit the radial symmetry in their computation.This has the disadvantage that non-spherical aberrations cannot easily be represented.The accuracy of the RW model is a function of the number of data points used in the numerical integration.3D SP-FFT* and SP-CZT* compute PSFs slice by slice for each  position.Using a code profiler under Matlab, we observed that about 70% of the computational time of the SP-FFT* technique is dominated by the FFT operation.The algorithmic complexity of each 2D slice in the SP-FFT* is given by  (    log(    )).The computational time per voxel for different window sizes therefore remains constant, while the total computation time of the volume PSF scales linearly with the total number of pixels.
Reflecting on Eq. ( 6), the CZT operation requires 3 FFTs to compute.In the SP-CZT* method, a new and bigger window is calculated such that the field at higher depth would not suffer from wrap-around effects.The run-time depends on both, lateral size and −depth.This explains the relative higher computational time per voxel in SP-CZT* compared to SP-FFT*.
The Fourier shell method (F-Shell*) computationally depends on the interpolation of the coefficients of   at different sub-pixels.By analysing the algorithm with the code profiler under MATLAB, this interpolation process takes up about 70% of the total computation cost.If the −range becomes significantly larger, a large number of coefficients also needs to be interpolated.The computation time therefore mainly depends on the number of iterations and the kernel size.Those coefficients are stored on the disk of the computer and can be accessed easily without any heavy computation, reducing the computation time of a second run by up to 20×.
In the computation of the Sinc-R* PSF, the Fourier wrap around in the calculation of wave propagation is avoided by starting from the real-space (sinc()) representation of a Fourier-shell and limiting and modifying it further in Fourier space.Artefacts that may be present in the calculation from the edge of the grid during this 3D FFT.This can be avoided by damping the edge or extending the window by 25%.The difference in the computation with or without this 25% extension is not significant.It is advisable to consider it in the computation as more information can be preserved by it.It should be noted that modifications in Fourier-space lead to convolutions in real space, which can also cause wrap-around artefacts, but these seem to be relatively minor.In the Sinc-R method, the required −range needs to initially be highly sampled leading to an heavy computation, which is not required for the other propagation methods.
The state-of-the-art GL method discussed here is implemented using a Bessel series approximation of the GL method [2] over a single  slice.A radial asymmetry in the PSF is assumed and a piece-wise linear interpolation is used to compute the whole 3D volume PSF on a rectilinear grid.The model's computational time and accuracy are inversely dependent.Both depend on the number of basis in the Bessel series and the sampling number.The computation time of the PSF corresponds to the default number basis vectors, which is equal to 100 Bessel functions and number of coefficient parameters equal to 1000.These make the GL PSF calculation the fastest, albeit not very precise scalar method.
sPSF, vPSF and RW compute the numerical integration of the electric field in the image plane using Simpson's rule [6].The methods are relatively fast and the computational time per voxel seems to be independent of the −range.Lastly, the software PSFGen is the most expensive in time and in memory among the PSF models, which we compare to in this work.

Conclusion
In this work, we presented a general approach for calculating the 3D PSF of a system satisfying the Abbe sine condition using Fourier based techniques.The PSFs models, with other state-of-the-art PSF calculation methods, are compared to a defined gold standard, which consists of an highly sampled Richard and Wolf model.We account for the digital Fourier transform pitfalls in the calculation and present different strategies to avoid them.The SP-CZT* model is proven to be the most accurate among our models.This accuracy comes at a cost in terms of its computation particularly at large field depths.The SP-FFT* and F-Shell* have similar accuracy and behaviour.In their computation, standing waves caused by wrap-around (which can also be interpreted as undersampling of the phase near the edge of the pupil) in the PSF are not fully avoided when the calculation grid is insufficient for the given depth of field but reduced.These two methods are however accurate and fast enough when a large depth of field is not required.The Sinc-R* model presents an attractive feature for easy implementation of a volume PSF largely avoiding the wrap-around problem.Overall the methods described in this work satisfy the physical condition of the imaging better than the state-of-the-art PSFs.They are already computationally efficient compared to the state-of-the-art given the fact that there is no use of radial symmetry to speed up their calculation.Our PSFs models can accommodate any non-circular aberration and can easily be modified to account for a tilted plane such as a tilt in a coverslip for instance.The models are also reproducible and can be easily extended to represent different imaging modalities.An experimental validation of the PSF models presented in this work is essential to verify the experimental accuracy of these models.

Fig. 4 .
Fig. 4. Visualization of the pupils and representation of the Huygens Wavelet becoming a plane wave.

Fig. 6 .
Fig. 6.Field distribution calculated from the Fourier transform of (a) a hard aperture of size 128 × 128 pixels, (b) a jinc aperture aperture of size 128 × 128 pixels, (c) a hard aperture of size 1024 × 1024 pixels cropped to 128 × 128 pixels size for display and, (d) a jinc aperture aperture of size 1024 × 1024 pixels cropped to 128 × 128 pixels size for display.A DampEdge of 15 % is applied to the generated field (full size) using the jinc-trick in (b) and (d).Figures are displayed at tan −1 () and centered at the zero of the display,  being the field distribution and  = 20.

Fig. 7 .
Fig. 7. Wrap-around of using FFTs in PSF calculations.Profiles displayed at  = 0.05 of the PSFs calculated using the slice propagation method.(a, b, c) -plane at defocus position 7 µm.(d, e, f) -cut parallel to the optical axis at 3.5 µm away from the focus and with depth range Δ = 3.5 µm.(a, d) Directly applying the propagator in Fourier space.(b, e) By zero-padding the image window size to twice its original size.(c, f) Using the chirp Z-transform operator.The parameters are NA = 1.4,immersion medium: water of refractive index 1.33, polarization: circular; emission wavelength  em = 580 nm, voxel size 80 nm × 80 nm × 140 nm and, displayed window size: 256 × 256 × 25 pixel.

Fig. 9 .Algorithm 2
Fig. 9. Diagram of the wave propagation in real space for calculating the zoom factor .

Fig. 10 .
Fig. 10.(a) Phase shift in the ideal wave exp(2  ).(b) Phase at a sub-pixel 0.25, indicated by the horizontal white line in (a).(c) Interpolation table in Fourier space containing the interpolator coefficients at 60 different sub-pixel positions.(d) Interpolation coefficients in Fourier space along the 0.25 sub-pixel indicated by the horizontal white line in (a) and (c).

6 :
Calculate the three components of the field in the McCutchen pupil following Step 1 to 4 of Algorithm 1 7: Pre-compensate for the projection factor 1/cos  of the sinc-shell by multiplying the field with cos ,  being the angular aperture that is previously defined.Apply the factor 1/ √ cos  as described in Step 5 of Algorithm 1 8: Multiply the resulting spectrum from Step 7 with the sinc-shell from Step 5 9: Perform a 3D FFT for each of the field components to obtain the sought-after field components ì ℎ   in real space 10: Calculate the PSF intensity: ℎ = | ì ℎ   | 2 11: Extract the field within the desired window [  ,   ,   ]

Fig. 11 presentsFig. 11 .
Fig. 11 presents −slices of the PSFs at  = 3.2 µm from the focal position and in-focus ( = 0 µm) and the −cross sections of the PSFs at  = 0 µm.Severe standing waves are observed at the −depth equal to 3.2 µm for the SP-FFT* and F-Shell* methods.The smallest region closest to the centre of the image in the Sinc-R* at 3.2 µm also does not appear to be perfectly circular.

Fig. 12 .
Fig. 12. Radial mean profiles, in logarithmic scale, of the PSFs at the focal plane.

Fig. 13 .
Fig. 13.Error map of the simulated PSFs relative to the highly sampled gold standard RW-GS.

Fig. 15 .
Fig. 15.(a) Representation of the optical transfer function (OTF) calculated using RW-GS and displayed at gamma of 0.2.The missing cone is presented in yellow.(b, c, d) Observation of the violation of the missing cone by calculating and presenting the integrated intensity along  of the PSFs.The profile which corresponds to RW-GS is displayed in all figures (b), (c), and (d) in red solid line as a reference.

Fig. 16 .
Fig. 16.Standard deviation of the integrated intensity over the 20 % around the focus of the given −range.

Fig. 17 .
Fig.17.Quantification of the wrap-around effect in dependence of the size of the calculation window ( and ).The red dashed horizontal lines pass by the error  of the RW for each given window size.The reference window is calculated to be 337 × 337.

Funding.
This work was funded by the DAAD through the African Institute for Mathematical Sciences and Stellenbosch University, and Friedrich Schiller University Jena.This work was also supported by the German Research Foundation (DFG) through the Collaborative Research Center PolyTarget 1278, project number 316213987, subproject C04 and the Council for Scientific and Industrial Research (CSIR), project number LREQA03.
The Fourier-based PSF models that are presented in this work are based illustrates this effect.Fig.8.Vector field propagation from a circular aperture of diameter  through a high NA system.ì   , ì   are unit vectors of  and −polarization state, ì   represents the unit vector of the −polarization after the projection on the Gaussian reference sphere.