Vectorial phase retrieval in super-resolution polarization microscopy Open Access

Fluorescence microscopy is a widely used imaging modality in biological research¹ given its strong signal, selective labeling within complex systems,² and compatibility with super-resolution methods.³ Moreover, this technique also allows access to the sample’s structural properties,^4,5 making it very useful for studying biomechanics at the molecular level. For example, in single-molecule orientation localization microscopy (SMOLM), 3D spatial localization can reach a precision of a few nanometers while simultaneously allowing the characterization of the molecule’s 3D orientational behavior (e.g., mean orientation and degree of wobble).⁵ Common SMOLM techniques include polarization channel splitting^6–9 and point spread function (PSF) engineering,^10–14 which can be used together or separately. The shape of the PSF can change considerably with the emitter’s orientation and longitudinal position, so it is crucial to take this into account to enable a full estimation of the parameters^13,15,16 and to avoid localization biases.¹⁷ In that respect, PSF engineering techniques have the aim of enhancing these shape changes. Nevertheless, any optical aberration, polarization distortion, or misalignment in the imaging system can affect the final shape of the PSFs and, thus, lead to an inaccurate estimation of the parameters. In particular, polarization aberrations are delicate to correct as they require additional adaptive strategies that account for the vectorial nature of light propagation.^18,19

A common solution to this problem is to perform a set of calibration measurements^20,21 and use them in a phase retrieval optimization algorithm to determine the aberrations present in the system. For this approach to work, it is important to have an accurate model for light’s propagation from a known source to the camera and to incorporate phase diversity, e.g., through measurements at varying focal planes, both to avoid falling into local minima and to accelerate convergence.²² For single-molecule localization microscopes, the standard approach is to measure the PSFs generated by fluorescent nanobeads in varying focal planes and use this Z-stack of images in a phase retrieval algorithm.²³ Initial approaches relied on scalar models, assuming a point source emitting a spherical wavefront along with a scalar pupil representing the aberrations. In this simplified case, Gerchberg–Saxton iterative algorithms^24,25 can be used, thus reducing the complexity of the implementation. However, these algorithms are less flexible since the pupil model cannot be chosen freely, and they are not directly applicable to full vectorial models that consider the dipolar emission pattern and polarization distortions.^20,21,26 A more flexible approach is offered by casting the phase retrieval problem as a nonlinear optimization routine where most parameters in the model can be found during the optimization process, although this requires providing analytical formulas for the gradients in these parameters, hence complicating implementation.²⁵ More accurate models that take into account the vector nature of the emitted light have also been proposed.^27,28 They incorporate the effect of the interface between the medium embedding the fluorescent particles and the coverslip, which causes extra aberrations, polarization-dependent transmission, and supercritical angle fluorescence (SAF) radiation.^28–30 However, these approaches assume a scalar mask to characterize all the remaining aberrations, thus preventing them from correcting any residual polarization-dependent effects. Moreover, the polarized components of the emitted fields are eventually summed in order to model an unpolarized measurement, which is not directly applicable to polarized PSFs in SMOLM.

Recently, new SMOLM techniques have been proposed that use birefringent elements either to efficiently encode the 3D orientation and 3D localization information of the emitting dipole into the shape of the two polarization components of the PSF^10,12,13,31 or to understand the intensity and/or shape of their different polarization projections.^8,9,32 For such approaches, it is essential to take into account the emission pattern of the dipole source, its interaction with the interface between the embedded medium and coverslip, and any polarization-dependent aberrations.³³ To address these issues, a model is used here where all vector aspects of the propagation of the light emitted by the source to the back focal plane (BFP) are taken into account and where the aberrations and polarization distortions are represented by a birefringent distribution at the pupil plane (BDPP) modeled with a spatially varying Jones matrix.³⁴ It is shown that in order to properly characterize this BDPP, it is necessary to introduce polarization projection diversity, obtained by projecting the PSFs onto various polarization states, in addition to the phase diversity given by changing the location of the focal plane. The PSF images generated with these two diversities form a PZ-stack that is fed into a nonlinear optimization algorithm that allows retrieving the unknown BDPP. This approach allows for the inclusion of many parameters that are necessary for proper characterization, such as photo-bleaching amplitudes, background illumination, or diversity-dependent tilts. Accurate and computationally amenable models for the light produced by fluorescent beads (commonly used for calibration measurements) are also included. These models take into account the unpolarized nature of the emitted light and the blurring due to bead size.³⁵ The software package pyPSFstack used for the modeling of PZ-stacks and the phase retrieval process can be found in Ref. 36. The retrieval algorithm was implemented with the neural network framework PyTorch,³⁷ which greatly simplifies its implementation and flexibility due to the automatic computation of all gradients and its integration with GPUs if available. While the emphasis of this work is on the characterization of BDPPs, both the theory and code are equally applicable to scalar-only aberration pupil distributions, as shown in the supplementary material.

The goal of this work is to provide a method for determining from a set of calibration PSFs the system’s BDPP, represented by a Jones matrix of the form in Eq. (4). The various functions in this expression must first be expanded in terms of a basis, whose expansion coefficients are then determined through nonlinear optimization. Common basis choices for optimization are the Zernike polynomials⁴⁹ and a direct pixel representation,²⁸ both of which result in comparable speeds since they employ the same number of discrete Fourier transforms, the most costly operation even when implemented through the fast Fourier transform (FFT) algorithm.

It is common practice in phase retrieval algorithms for optical microscopes to assume access to a stack of intensity images for varying focal distances Δ_ζ from the location of the best focus (a Z-stack). The phase diversity provided by the varying focal distances is taken into account by the phase factor $D(Δζ)(u)$ defined in Eq. (2) within the Green tensor in Eq. (1). This additional information, referred to as phase diversity, helps the algorithm converge to an appropriate solution without falling into local minima, as well as discriminate between vortices with opposite topological charges.

Before comparing the PZ-stack computed by the model presented thus far with the one measured experimentally, it is necessary to account for other effects. First, depending on the size of the fluorescent bead, it might be necessary to include a blurring effect. As shown in Ref. 35, the exact three-dimensional blurring corresponds to a superposition of two-dimensional convolutions between the PSFs generated by point sources located along the longitudinal diameter of the bead with a kernel that weights more heavily the contributions near the center of the bead and vanishes for those at the poles. This exact blurring cannot be rewritten as a three-dimensional convolution due to the lack of translation invariance along the axis of propagation, but it can be approximated through a semi-analytic method based on a Taylor expansion. By keeping the first term of the expansion, we obtain a two-dimensional blurring model based on a two-dimensional convolution of the PSFs. If instead we keep two terms in the expansion, we get a model for the blurring along both the transverse and longitudinal directions, which is computed via two-dimensional convolutions of the PSFs and their second-derivatives with respect to the bead’s axial position, thus giving a three-dimensional blurring model. Both of these approximate models are computationally more efficient than the exact one and have been implemented in pyPSFstack. Moreover, they can be used during the BDPP retrieval process, as shown in Sec. V B, albeit at a computational cost due to the need for supplementary Fourier transforms. Nonetheless, they allow for the analytic computation of the gradients, hence significantly limiting the computational slowdown; using readily available functions to blur images, on the other hand, such as the one used in Ref. 28 with a Gaussian kernel, would require the gradients to be computed via finite differences.

The goal of the nonlinear optimization routine is to find the set of optimization parameters in the forward model that minimizes the cost function, assessing the differences between the measured and modeled PZ-stacks. This is achieved by supplying an optimization algorithm, such as Adam⁵⁶ or L-BFGS,⁵⁷ with a function that uses all the current values of the parameters to compute the forward model all the way to the value of the cost function. In addition, one should supply the optimization algorithm with another function that performs a backward computation to obtain the gradients of the cost function with respect to all the optimization parameters. These gradients are used to change the parameter values until a minimum is reached. The gradient computation is straightforward but tedious and can be achieved by following the rules outlined in Ref. 58. However, an advantage of implementing nonlinear optimization with the neural network framework PyTorch is that only the forward model must be implemented explicitly since the backward model for computing the gradients of the various parameters is computed automatically. This framework also offers the most common optimization algorithms.

The implementation using PyTorch greatly simplifies the vectorial phase retrieval algorithm, as shown in Algorithm 1. First, the system parameters and the initial values for the optimization parameters are used to initialize the forward model used to compute the estimated PSFs and evaluate the current value of the cost function. For all the examples presented in this work, the BDPP was initialized as the identity matrix, i.e., $c0(0)=1$ and $cl(W)=cl(j)$ for all j, l ≠ 0. However, this can be changed if there is prior knowledge about the BDPP. After initialization, the phase retrieval algorithm enters a for loop in which the forward model is used to compute the cost function using the current values of the parameters. Then, automatic differentiation is called for to obtain the gradients. Finally, the hyperparameters of the chosen optimization algorithm and the values of the gradients are used to update the values for all the optimization parameters. This process is iterated N_iter times, after which the process is completed. Note that other criteria can be used to terminate the process; for example, the process could be terminated when the change in the cost function falls below a set limit. All the results presented in this work were obtained using the Adam optimizer since it is more robust than simple gradient descent and faster and less memory intensive than L-BFGS. While this optimization algorithm has several parameters that can be tuned, in practice, it was found that choosing the appropriate learning rate was sufficient to obtain satisfactory solutions. The convergence was assessed by analyzing the evolution of the cost function. A more detailed explanation can be found in the examples in the repository.³⁶ 

To implement the forward model, all calculations are performed numerically using FFTs, requiring all spatial quantities to be sampled consistently with the camera’s pixel pitch p_cam, so that the size at the BFP is fixed to L_pupil = λM/(p_camn_f). The total size at the BFP is divided into N × N points used for the computation, and the resulting two-dimensional sampling is labeled with the lexicographic index ℓ, which takes over all spatial dependence in terms of u and ρ, for instance, $G0(u)→G0(ℓ)$⁠. All the parameters and steps necessary to compute the forward model are summarized in Algorithm 2.

The results of the procedure are shown in Fig. 3, where the Jones matrix for the ground truth is compared to the ones retrieved from the PZ-stacks. This figure also shows the PSFs formed by dipoles oriented along the three Cartesian axes and by an unpolarized dipole (constructed as an incoherent mixture of the three orthogonal dipoles). The PSFs shown are modeled using the standard projectors $P1$ and $P2$ [Eqs. (7) and (8)] for each birefringent mask, for the SEO, the output is projected onto the left and right circular polarizations, while for the q-plate, the output is projected onto the horizontal and vertical polarizations. To assess the accuracy of the retrieval, both the correlation between the retrieved and true BDPPs as well as the root-mean square (rms) error, ϵ_rms, between the retrieved and true PSFs are shown in Fig. 3. For the SEO, the retrieved BDPP and the corresponding PSFs are almost indistinguishable from those corresponding to the ground truth. However, for the q-plate, there are appreciable differences between the original and retrieved BDPPs. The deviation at the center is due to the chosen model, since the Zernike polynomials struggle to reproduce the singularity at the center of the q-plate. In general, this type of singularity does not arise as an aberration but rather by design of the mask, so that it can be incorporated as prior information. For instance, on top of the Zernike decomposition, a q-plate could be added to the model with its center and orientation as optimization parameters. Another noticeable difference in the case of the q-plate is that the phase between the two rows of the Jones matrix is not correct. This error happens only for birefringent masks that generate rotationally symmetric PSFs for unpolarized emitters for all polarization projections. For these masks, the algorithm cannot determine the global phase between the rows of the corresponding Jones matrix. This problem can be solved by placing the device that introduces the polarization diversity before the birefringent mask. Even when this solution is not used, the correlation (when neglecting the phase difference between the rows) is quite high, and the reproduced PSFs for any dipole orientation are again indistinguishable from the true ones. Therefore, these differences are inconsequential for the purposes of single molecule fluorescence microscopy. In addition, the error introduced for z₀ and α has no impact on the retrieval as long as the defocus Δ is included as an optimization parameter. Similar results using a pixel-based model are shown in the supplementary material, Sec. VI, as is a comparison of the convergence of the cost function between the two models.

For comparison, the retrieval of the same birefringent mask is also performed without using polarization diversity. To make this comparison consistent and fair, the number of photons reaching the camera and those in the background illumination is doubled since there is no polarization separation. Moreover, the phase diversity images are now taken from −250 to 250 nm of the RFP with a reduced step size of 5 nm, so that the total number of PSFs used is larger than in the previous case. The results are shown in Fig. 3, where the algorithm is seen to fail to retrieve the appropriate BDPP, and thus, it cannot generate the correct PSFs when they are projected onto a given polarization state. These results show the need for polarization diversity when characterizing a BDPP. A more detailed comparison between the scalar and vector models is provided in the supplementary material, Sec. IV. Note that the runtime is much longer due to the increase in the number of FFTs required. (This runtime should be comparable with that of a model that includes tilts but no blurring; see Table I.)

As an additional test of the present implementation, the retrieval of a BDPP from a highly blurred PZ-stack such as the one shown in Fig. 4 is considered. Using the SEO as an example, a PZ-stack is constructed as in Sec. III B and III C with the significant difference of increasing the radius of the nanobead, which is chosen randomly from a uniform distribution between 150 and 170 nm. Here again, the distance to the coverslip is taken to be equal to the radius of the bead. As shown in Fig. 4, this bead size causes a significant blur in the resulting PSFs. The blurred PSFs were computed with the exact blurring presented in Ref. 35, but this approach turned out to be computationally expensive and not necessary for the retrieval procedure, as shown in what follows.

A first attempt to retrieve the BDPP is performed by completely neglecting the blurring effect. However, as shown in Fig. 4, this approach fails to retrieve the appropriate BDPP, showing that bead size effects cannot be neglected in this case. Next, as discussed in the supplementary material, Sec. SIII, the two approximate models presented in Ref. 35 were used for the retrieval, where the radius of the bead is used as an optimization parameter for the blurring with an initial value of 150 nm. The first of these approaches is a 2D convolution with a kernel given by a spherical Bessel function. This approach produces a satisfactory result, whose retrieved BDPP has a correlation of 98.6% with the ground truth. The second is the semi-analytic model described in Sec. III C for reproducing the effects of the three-dimensional blurring, based on a Taylor expansion around the center of the bead. In this case, a BDPP indistinguishable from the true one is again retrieved with a slightly better correlation of 99.5%. As shown in Table I, both blurring models significantly increase the time per iteration due to the required supplementary Fourier transforms. This increase is particularly important for the 3D blurring model, which requires the computation of the first and second derivatives of the Green tensor with respect to the bead’s axial position, hence tripling the number of Fourier transforms in addition to two convolutions at the image plane. Therefore, the small gain provided by the 3D model might not justify the extra computational resources needed for it. Nonetheless, both models are able to perform the full retrieval in well under a minute with the CPU and in a few seconds with the GPU.

After validating the proposed retrieval procedure on simulated data, we apply it to retrieve a BDPP distribution from a PZ-stack measured experimentally. We use a sparse sample of fluorescent nanobeads with a diameter of 20 nm (orange carboxylate-modified FluoSpheres), immobilized on the surface of a poly-L-lysine-coated coverslip, and embedded in water. The sample is mounted on a XYZ piezo stage (Physik Instrumente) and is excited by a continuous wave laser emitting at 561 nm (Oxxius L4Cc) in a wide-field illumination configuration using an oil immersion objective lens (APO TIRF ×100, NA = 1.49, Nikon). The emitted fluorescence is collected by the same objective lens and then passes through two multiband dichroics (Semrock, R405/488/561/635-t1-25x36) and a fluorescence filter (Semrock, 605/40). A telescopic relay system (composed of two achromatic doublets with f = 250 mm so that the magnification is unity) is used to access the BFP of the objective. The different polarization projections are taken by using a QWP (AQWP05M-600, Thorlabs) placed on a rotating mount (Newport, PR50CC), followed by a quartz Wollaston polarizing 2.2° beamsplitter (Edmunds, 68-820). The final images are measured using an ORCA Fusion-Digital CMOS C14440-20 UP (1024 × 1024 pixels, 6.5 × 6.5 μm² pixel size, Hamamatsu). PZ-stacks were acquired with a step size of 50 nm and a rotation step for the QWP of 30°. Figure 5 shows part of the experimental PZ-stack.

Before launching the retrieval on the experimental PZ-stack, it should be noted that there are several factors that lead to the introduction of a diversity-dependent phase tilt at the BFP. First, the use of a Wollaston prism to spatially separate the two polarization components onto different sections of the camera might make it difficult to have the same center for the PSFs for each polarization component. Second, any slight wedge on the rotating QWP introduces a tilt that rotates with it. Finally, any slight misalignment between the stage moving the sample and the optical axis defined by the microscope objective leads to a defocus-dependent tilt. Therefore, it is best to introduce in the forward model outlined in Sec. IV extra optimization parameters to independently adjust these tilts at the BFP for each combination of diversities. In particular, steps 4 and 5 in Algorithm 2 should be reversed in order to apply the polarization diversity before propagating to the image plane, and the following additional step should be added after applying the polarization diversity:

• Apply a phase tilt $T(ℓ)=expi2πt(ζ,p)⋅ℓ$ to each diversity. Optimization parameters: $t(ζ,p)=(tx(ζ,p),ty(ζ,p))$⁠.

The downside of including these extra parameters is that the number of FFTs needed for each iteration increases by a factor equal to the number of polarization projection steps (see Table I). However, the total run time remains quite manageable even if the 2D blurring model is used.

All the optimization parameters for the retrieval procedure were used, except for the bead radius (for the blurring), since it can be neglected due to the small bead size (R_b = 10 nm). The results of the retrieval process are shown in Fig. 5, where the strong agreement between the measured PZ-stack and the one modeled with the retrieved BDPP (see the supplementary material, Sec. SVIII, for a full comparison) can be appreciated. Moreover, the presence of the SEO is visible in the retrieved BDPP. In this case, since it is known that there is an SEO at the BFP, it is worth trying to separate the total retrieved BDPP into a misaligned ideal SEO and a BDPP containing the scalar and polarization aberrations of the system. It should be noted that scalar tilts and defocus are not shown as part of the scalar aberrations. From this decomposition, it can be seen that the largest birefringent contribution comes from the SEO, but the aberrations are not negligible and must be taken into account (see the supplementary material, Sec. SVIII, for more details). Note also that the polarization aberrations show larger variations than the scalar ones, which are almost flat, showing the need to consider polarization aberrations when using polarization-dependent systems. Given that for retrieval from experimental data, the ground truth is not known a priori, the same retrieval was performed for other beads within the field of view to corroborate the results. The retrieved BDPPs are indeed highly correlated with the one shown in Fig. 5 (see the supplementary material, Sec. SVIII, for full details).

A methodology and phase retrieval algorithm were presented for the characterization of birefringent distributions at the pupil plane (BDPP) from stacks of PSFs. In particular, it was shown that, for polarization-dependent systems, aberrations should be modeled as a BDPP encoding both scalar and polarization deviations from the ideal system, and that the use of polarization diversity is essential for its proper characterization. The software program pyPSFstack created and used for the modeling and retrieval presented in this work is freely available. In particular, the BDPP retrieval implementation based on PyTorch makes this software flexible and customizable. This is shown in this work by incorporating several optimization parameters apart from those used to describe the BDPP, as well as the blurring models presented in Ref. 35. Even when considering these extra parameters and models, the runtime remains manageable. In addition, this software can also be used for the retrieval of scalar BDPPs, as shown in the supplementary material in Sec. SVII for a mask known as the tetrapod.⁵⁹ 

While it was assumed that the sources used for the characterization were unpolarized, the retrieval model and algorithm presented here can be easily adapted to sources that are dipolar, partially polarized, or with a fixed orientation, such as molecules fixed on a surface¹⁴ or in DNA origami.^60,61 Similarly, this model can also be used to find optimal scalar or BDPPs minimizing a particular combination of the Cramèr–Rao bounds for the parameters to be extracted from the shape of the PSFs. Finally, note that there are some modifications that could be implemented in the retrieval algorithm that might allow faster convergence, such as the use of stochastic gradient descent, where only a given set of measurements is used during each iteration or spectral initialization.⁶² 

See the supplementary material for the exact expression of the Green tensor, an analysis of the various criteria used to define the position for the best focus, the blurring models, a detailed comparison between the scalar and vector models, the effect of noise and the number of measurements composing the diversitiy in the performance of the retrieval, examples using a pixel-based model for the retrieval of birefringent masks, the retrieval of a scalar mask, and further corroboration of the retrieval from experimental data.

R.G.-C. acknowledges S. W. Paine and S. M. Popoff for useful discussions. L.A.A.-C. acknowledges M. Sison for the experimental advice. The authors also thank T. G. Brown for supplying the SEO. R.G.-C. acknowledges funding from the Labex WIFI (Grant Nos. ANR-10-LABX-24 and ANR-10-IDEX-0001-02 PSL*). L.A.A.-C., S.B., and M.A.A. acknowledge funding from Grant No. ANR-21-CE24-0014 and S.B. from Grant No. ANR-20-CE42-0003.

The authors have no conflicts to disclose.

R. Gutiérrez-Cuevas: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). L. A. Alemán-Castañeda: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Supervision (equal); Validation (equal); Writing – review & editing (equal). I. Herrera: Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Validation (equal); Writing – review & editing (equal). S. Brasselet: Conceptualization (supporting); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). M. A. Alonso: Conceptualization (supporting); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Software (supporting); Supervision (equal); Validation (equal); Writing – review & editing (equal).

Data underlying the results presented in this paper are available in Ref. 36.