1. Introduction

Polarimetry is the measurement of the polarization state of light and/or the polarization properties of materials. Such measurements are usually characterized in terms of the Stokes parameters and the Mueller matrix, respectively, which are directly accessible from measurements of the intensity. Imaging polarimetry, in which the polarization is measured as a function of position, is particularly important in applications ranging from microscopy to remote sensing.

Conventional techniques for Stokes polarimetry require multiple intensity measurements, either through time-sequencing or by splitting different polarization components into several separate detection channels [1,2]. For example, one common method uses a rotating quarter-wave plate (QWP) followed by a fixed linear polarizer, in which the Stokes parameters are deduced from successive intensity measurements with the QWP oriented at different angles [3]. While these techniques can produce highly accurate measurements, they can be relatively complicated and/or time-consuming, generally involving moving parts or multiple beam paths.

When a short acquisition time is desirable, rapid polarization measurements may be taken using single-shot polarimetry, in which the Stokes parameters are estimated from a single intensity measurement. A variety of methods exist for single-shot polarimetry, involving gradient index lenses [4], patterned nanoscale gratings [4], or a split aperture composed of multiple polarizers [5]. One particularly simple method, referred to as star test polarimetry, uses a spatially-varying birefringent mask (BM) followed by a uniform polarization analyzer in the pupil plane of an exit-telecentric imaging system. With an appropriately chosen birefringence distribution, the inserted elements can encode full polarization information into the shape of the PSF in the rear focal plane of the lens. This method has been demonstrated experimentally using a stress-engineered optic (SEO), which is a BK7 glass window subjected to stress with trigonal symmetry at its periphery, followed by a circular analyzer [6,7]. Natural applications for this approach are those in which the object is a sparse set of discrete points, such as in astronomy and confocal microscopy. A recent specific application was the real-time monitoring of the output polarization states of each core within a multicore fiber bundle for applications in medical endoscopes [8]. An extension of this technique is being implemented within the context of superresolution microscopy for determining the 3D position, orientation, and vibration of independent fluorescent molecules [9]. Star test polarimetry can also be applied for imaging continuous objects if these are discretized by placing a pinhole array in the object plane of an afocal $4f$ relay system with unit magnification [10,11]. In all these applications, the polarization state of each object point can then be deduced from the corresponding PSF in the image plane, provided that the PSFs from different points do not overlap significantly. The spatial resolution of the measurement is thus limited by the size of the PSF. Encoding polarization information in the PSF necessarily implies increasing its size, since more information is being included in it.

The goal of this work is to find the birefringence distribution of the mask in the pupil plane that maximizes spatial resolution by producing the smallest possible PSF while imposing reasonable restrictions to ensure that the object’s polarization is effectively encoded in the shape of the PSF. The optimal solution is found to give very similar results to those of the birefringence distribution of the SEO. The solution’s statistical performance (i.e., the expected error in the retrieved Stokes parameters) is assessed by calculating the Fisher information matrix for the measurement.

2. System layout and notation

Consider the system shown in Fig. 1, in which a spatially-varying, transparent, thin birefringent mask with Jones matrix $\mathbb {J}(\mathbf {u})$ is placed at the pupil plane of an exit-telecentric imaging system, where $\mathbf {u}=(u\cos \phi ,u\sin \phi )$ is the pupil coordinate normalized such that $u\leq \mathrm {NA}$, with NA being the system’s numerical aperture. Collimated light (e.g., from a distant localized source) with uniform polarization $\mathbf {E}_0$ is then incident on this BM in the plane of an aperture with pupil function $A(\mathbf {u})$ (binary or apodized) and is focused by a lens, producing a polarization-dependent PSF. The input polarization may then be deduced from the shape of the PSF by separating two orthogonal polarization components of the field, denoted by $\mathbf {e}_1$ and $\mathbf {e}_2$, before forming an image of one (or both) of them. For applications where photons are not scarce, only one of the two images is often necessary, so the extraction of the component in question (say, $\mathbf {e}_1$) can be performed with a polarizer. As will be discussed in Section 9, however, there are applications in which it is best to use all photons, so instead the two components ($\mathbf {e}_1$ and $\mathbf {e}_2$) are separated by using a polarizing beamsplitter or a Nomarski or Wollaston prism and then they are imaged separately [9].

 

Fig. 1. System layout for a star test polarimetry measurement of an unknown input field $\mathbf {E}_0$, illustrated for the case in which a single circular polarization component is imaged. The pupil plane is described by a radial coordinate $u=\sin \theta$ and an azimuthal angle $\phi$.

Download Full Size | PDF

In what follows, all calculations are performed in the circular polarization basis, so that $\mathbf {e}_1=(1,0)$ and $\mathbf {e}_2=(0,1)$, corresponding to right-hand circular (RHC) and left-hand circular (LHC) polarizations, respectively. This choice facilitates the comparison with the SEO (which is used with a circular analyzer) in Section 8. If other output polarization components were analyzed, the results would remain valid after minor modifications. It is assumed that the imaging lens is slow (around 0.05 NA) so that the PSFs are imaged onto several pixels. All calculations then assume the paraxial limit, and angle-dependent polarization effects at the lens can be neglected.

3. Birefringence distribution and point-spread function

We use here the notation introduced in [12] to describe spatially-varying birefringence, where the Jones matrix of the BM in the circular polarization basis may be written as

After transmission through the BM, the $\mathbf {e}_1$ and $\mathbf {e}_2$ polarization components of the field are given by $\mathbf {e}_1^{\dagger }\hspace {1pt} \mathbb {J}(\mathbf {u})\hspace {1pt}\mathbf {E}_0$ and $\mathbf {e}_2^{\dagger }\hspace {1pt} \mathbb {J}(\mathbf {u})\hspace {1pt}\mathbf {E}_0$, respectively, where the dagger denotes a conjugate transpose. The paraxial PSF $I^{(j)}(\mathbf {x})$ of each polarization component $\mathbf {e}_j$ is the squared modulus of the Fourier transform of the pupil distribution:

The PSF may be rewritten succinctly by introducing the pupil functions

4. Constraints on the BM distribution

As mentioned previously, the primary performance metric for a candidate BM distribution is its impact on the width of the resulting PSF. However, it is necessary to impose some constraints in order to avoid solutions that would be unsuitable for polarimetry. For practical purposes, it is useful for the PSFs corresponding to each polarization component to have the same total power $\Psi ^{(j)}=\iint \! I^{(j)}(\mathbf {x})\mathrm {d}^{2} x$. This makes the signal-to-noise ratio roughly independent of the input polarization, and it ensures that all polarization information encoded by the BM is fully contained within the shape of each intensity distribution rather than the overall power. Consequently, the input polarization can be deduced by analyzing the PSF of a single component, which is guaranteed to contain roughly half of the incident photons when the signal is sufficiently large.

The total power of each polarization component can be found by integrating Eq. (6) over $\mathbf {x}$ and applying Parseval’s theorem to each term, leading to

5. Derivation of differential equation and boundary conditions

When the BM is introduced into the system, the PSF must increase in size since it now contains polarization information. Because the PSF would provide information about four values (the Stokes parameters) rather than one, the width of $I^{(j)}(\mathbf {x})$ can be expected to be roughly twice as large as that of a system with zero (or uniform) birefringence. Of course, there are many measures for the size of the PSF, each of which would give a slightly different result. The metric chosen here is the RMS irradiance width, which allows a simple variational calculation. It is expected that the optimal birefringence distribution for this measure will be nearly optimal for other measures, such as the full width at half maximum (FWHM) or the Strehl ratio.

Generally speaking, the width of the PSF does not strongly depend on the incident polarization state since the normalized intensity contributions $\mathcal {I}^{(1)}_n(\mathbf {x})$ ($n=0,1,2,3$) have similar compositions. Therefore, for simplicity, the performance may be evaluated in terms of the spread of the average PSF over all possible input polarizations, i.e., over all values of $S_1$, $S_2$, and $S_3$ within the range $[-S_0,S_0]$. This average produces $I^{(1)}(\mathbf {x})=\frac {1}{2}S_0\mathcal {I}_0^{(1)}(\mathbf {x})$, which is the PSF for unpolarized light. (Since $\mathcal {I}_0^{(1)}(\mathbf {x})=\mathcal {I}_0^{(2)}(-\mathbf {x}$), the analysis is also valid when the PSFs of both polarization components are used, as discussed in Section 9.) As shown in [12], the squared RMS width $r^{2}$ of the PSF for unpolarized light is

Since $\vec {q}$ is a unit vector, it contains only three free components for optimization. A natural choice is to use the three-dimensional vector $\mathbf {q}$ and let $q_0=(1-\lvert {\mathbf {q}}\rvert ^{2})^{1/2}$. To find the solutions we use a variational approach: we consider a linear combination of the PSF’s width increase and the three constraints, namely $M=\kappa \,\Delta r^{2}+{\boldsymbol \Lambda }\cdot \langle {\boldsymbol \beta }\rangle _A$, where ${\boldsymbol \Lambda }=(\Lambda _1,\Lambda _2,\Lambda _3)$ is a vector of Lagrange multipliers and the constant prefactor $\kappa$ was introduced for future convenience. The variational procedure consists of finding the conditions under which $M$ is stationary under infinitesimal changes in $\mathbf {q}$. For this purpose, consider a variation $\vec {q}\to \vec {q}+\vec {\epsilon }$, where $\vec {\epsilon }=(\epsilon _0,\epsilon _1,\epsilon _2,\epsilon _3)$ is an infinitesimal change. Since after this change the vector must remain a unit vector, $\vec {\epsilon }$ must be perpendicular to $\vec {q}$, so only three of its components are independent parameters, e.g., the zeroth component can be chosen to be a function of the remaining ones according to

In Eqs. (14) and (15), $\epsilon _0$ can be expressed in terms of the remaining $\epsilon _n$ by using Eq. (13). The substitution of the resulting variations into $M$ causes a variation of the form $M\to M+\iint (M_1\epsilon _1+M_2\epsilon _2+M_3\epsilon _3)\,\mathrm {d}^{2}u$. The total variation of $M$ is then insensitive to the infinitesimal changes $\epsilon _n$ at any point inside the pupil only if $M_1=M_2=M_3=0$ everywhere inside the pupil, that is, if the following three constraints hold:

6. Solution ignoring the boundary conditions

In order to solve these equations, it is illustrative to first solve the simpler problem in which the constraints are ignored in the variational calculation (that is, the Lagrange multipliers are set to zero), as are the boundary conditions. In this case, Eqs. (16a) through (16c) can be written compactly as

Consider the standard case of a circular hard aperture with pupil function

The solutions in Eq. (21) turn out to automatically satisfy the constraints $\langle {\beta _1} \rangle _A=\langle {\beta _2} \rangle _A=0$. For the third constraint ($\langle {\beta _3} \rangle _A=0$), we calculate

 

Fig. 2. Value of Eq. (22) as a function of $b_m$ (times a scaling factor) for several values of $|m|$. Note that this quantity can vanish only for $1\leq |m|\leq 4$. The inset shows the corresponding values of $\Delta r$ for the eight values of $|m|$ and $b_m$ for which the condition $\langle \beta _3\rangle _A=0$ is satisfied.

Download Full Size | PDF

While the solutions found above satisfy the constraints $\langle \boldsymbol {\beta }\rangle _A=\boldsymbol {0}$, the fact that these constraints were ignored at the outset in the variational calculation means that the solutions cannot simultaneously satisfy the constraints and the boundary condition in Eq. (17), which for the hard circular aperture considered here (where the normal derivative $\partial _\perp$ reduces to a radial derivative) is given by

7. Optimal birefringence distribution

Let us now go back to the general Eqs. (16) while again considering a hard circular pupil. Drawing inspiration from the previous case, we propose a separable solution constrained to the $q_1q_2$ plane, namely $\mathbf {q}(u,\phi )=\sin \delta (u)\,[\cos (m\phi ),\sin (m\phi ),0]$ with $m\ne 0$, which automatically satisfies the constraints $\langle {\beta _1} \rangle _A=\langle {\beta _2} \rangle _A=0$. Notice from Eq. (2) that this ansatz corresponds to $q_0(u)=\cos \delta (u)$, $\Theta =0$, and $\Phi =m\phi$. The substitution of this ansatz into Eq. (16c) implies that $\Lambda _1=\Lambda _2=0$, and the resulting form of both Eqs. (16a) and (16b) gives a differential equation for the half-retardance $\delta$, which after changing variables to $v=u/\mathrm {NA}$ becomes

 

Fig. 3. Radial retardance distribution $\bar {\delta }(v)$ for the optimal BM solution ignoring the boundary conditions, the true optimum, and an SEO with stress coefficient $c=1.166/\mathrm {NA}$, plotted as functions of the normalized radial pupil coordinate $v$.

Download Full Size | PDF

8. Performance evaluation and comparison to SEO

As discussed in the introduction, star test polarimetry is implemented experimentally by using an SEO having three points of symmetric stress. The birefringence distribution near the center of an SEO turns out to be approximately given by the same general form as the previous two solutions, namely $\mathbf {q}(u,\phi )=\sin \delta (u)\,[\cos (m\phi ),\sin (m\phi ),0]$, with $m=-1$ and where in this case the half-retardance is given by $\delta (u) = c\hspace {.5pt}u$, with $c$ being a measure of the stress in the SEO [6,7]. Given the similarity of this birefringence distribution to the optimal one, the constraints $\langle \beta _1\rangle _A=0$ and $\langle \beta _2\rangle _A=0$ are also automatically satisfied. The third constraint reduces to

Table 1 provides a comparison of the performance metrics of the BM solutions found in the previous two sections and the SEO. The half-retardance distributions are also summarized. All three solutions correspond to $|m|=1$. (The SEO corresponds specifically to $m=-1$, but the results are the same for $m=1$, which corresponds to an SEO followed a uniform half-wave plate.) The optimal solution does outperform the two other solutions not only in RMS width increase $\Delta r$ (calculated through Eq. (23)) for which it was optimized, but also for two other performance metrics: the FWHM and the Strehl ratio. However, the differences in these metrics between the optimal solution and the SEO are of only between 5 and 10%. Also provided for comparison are corresponding metrics for a diffraction-limited imaging system where no BM is used (or equivalently where the BM is uniform). We can see that the cost of encoding polarization information into the PSF (on average over all possible polarizations) is an increase in its FWHM of about 60% and a reduction of the Strehl ratio by about 50%.

Table 1. Performance metrics for the PSFs resulting from using three BM distributions characterized by the half-retardance functions in the second column, as well as for a diffraction-limited system, when the incident field is unpolarized.

View Table

The functions $G(\mathbf {x})$ and $H(\mathbf {x})$ and the corresponding PSF contributions $\mathcal {I}_n^{(j)}(\mathbf {x})$ for the optimal BM are shown in Fig. 4. Note that the PSF contributions involve negative contributions for $n\ge 1$, but the total measured PSF, which corresponds to the superposition in Eq. (6), is always non-negative as a result of the constraint $S_0^{2}\ge S_1^{2}+S_2^{2}+S_3^{2}$. For the PSF contributions, the top/bottom row corresponds to the case where the right/left circularly polarized component emerging from the BM is focused on the detector. Notice that, for these solutions, $\mathcal {I}_n^{(j)}(\mathbf {x})$ for $n=0,1,2$ are independent of $j$ (the circular component being analyzed) while $\mathcal {I}_3^{(j)}(\mathbf {x})$ changes by a global sign. That is, the symmetries of the functions $G(\mathbf {x})$ and $H(\mathbf {x})$ are such that the only sign change ($\pm$) that has an effect in Eqs. (7) is that in Eq. (7d). The equivalent plots for an SEO are virtually indistinguishable and are therefore not shown. Instead, to appreciate the very subtle differences, plots of the horizontal cross-sections $\mathcal {I}_n^{(1)}(x,0)$ for both the optimal BM solution and the SEO are shown in Fig. 5. In practice, the effect of these differences is probably insignificant compared with variations that would arrive, for example, from the pixelization of the detector.

 

Fig. 4. For the optimal BM solution: complex fields $G$ and $H$ (left, with amplitude cross-sections through the center shown on the right), and PSF contributions $\mathcal {I}_n^{(j)}(\mathbf {x})$ to the PSFs of the $\mathbf {e}_1$ (top row) and $\mathbf {e}_2$ (bottom row) output polarization components. The plots are shown over a square region with half-width $1.25\lambda /\mathrm {NA}$.

Download Full Size | PDF

 

Fig. 5. Horizontal slices through $y=0$ of the intensity contributions $\mathcal {I}_n^{(1)}$. The solid and dashed curves correspond to the optimal birefringence distribution and an SEO, respectively.

Download Full Size | PDF

The accuracy of the measurement can be assessed by using the Fisher information matrix in Eq. (8), which given the constraints imposed on the solutions reduces to

 

Fig. 6. Two-dimensional cross-sections of the error ellipsoids for incident polarization states in the Poincaré sphere’s $s_1$-$s_2$ and $s_1$-$s_3$ cross-sections (shown only for $s_1\ge 0$ due to symmetry) for (a-b) the optimal BM and (c-d) an SEO, when only the $\mathbf {e}_1$ output polarization component is measured with $\mathcal {N}=1500$ detected photons (red) and when both polarization components are measured with $\mathcal {N}=3000$ detected photons (blue). The error ellipses for unpolarized light are magnified in the inset of each plot.

Download Full Size | PDF

The fact that all solutions described so far have the form $\mathbf {q}(u,\phi )=\sin \delta (u)\,[\cos (\phi ),\pm \sin (\phi ),0]$, where $\delta (u)$ is linear near the center, can be understood in terms of the two circular components of the incident light. Recall that the BM is followed by an RHC polarizer and a Fourier-transforming lens before reaching the detector. Since the birefringence is smallest at the center of the BM and grows away from it, the BM/RHC polarizer combination acts as a rotationally-symmetric apodizer for the incident RHC-polarized light. The incident LHC-polarized component, on the other hand, gets converted into RHC light away from the BM’s center, and acquires a phase vortex due to a geometric phase effect. After Fourier transformation by the lens, the field distributions at the detector plane for these components are precisely the distributions $G$ and $H$ (shown in Fig. 5 for the optimal BM), whose forms are a localized central lobe with uniform phase and a slightly larger “donut” with a phase vortex of unit charge at its center, respectively. The radial extent of the latter would increase with $|m|$, the magnitude of the vortex charge written by the BM on $H$, so $|m|=1$ guarantees the most compact vortex distribution. The superposition of these two fields then gives access to combinations with all relative phases between them, as well as with all relative amplitudes due to their different radial dependence.

The accuracy of a polarimetric measurement is known to increase with the volume within the Poincaré sphere enclosed by a simplex whose corners correspond to the polarization components being measured [14]. The analysis above tells us that, in theory, each point of the PSF being measured corresponds to the measurement of a different polarization component, and that these cover the complete surface of the Poincaré sphere in a fairly uniform way, so that all the interior of the sphere is enclosed. The uniformity of the coverage of the Poincaré sphere in the $s_3$ direction following from the variation in the radial direction of the relative weights of the two field components is examined in Fig. 7(a) for both the optimal solution and the SEO. While the uniformity is not perfect, all points are represented in similar amounts and the centroid of the distribution is zero, as guaranteed by the constraint $\overline {\mu _3}=0$. In practice, of course, the sampling of these different polarization components is compromised by detector pixelation, and even if the pixels were to be made very small, this would come at the cost of each detecting less photons. However, it is clear that this type of measurement provides a more “democratic” coverage of the sphere when compared, say, to a polarimeter where six polarization components are measured.

 

Fig. 7. Histograms of the power coverage of different values of $s_3$ (characterized by $\mu _3$) over the extent of the PSFs, for the cases when (a) only the PSF for the $\mathbf {e}_1$ component is used, and (b) the PSFs for both components are used.

Download Full Size | PDF

9. Results when both polarization components are measured

So far we have assumed that only one circular polarization component ($\mathbf {e}_1$) emerging from the BM is focused to form an image. However, as mentioned earlier, it is possible to instead separately image both circular components and analyze the two resulting PSFs jointly. While doing so requires a more complex system, there are applications in which it is advantageous [9]. One clear advantage is that the number of detected photons is essentially doubled, reducing the uncertainty in the estimation of the Stokes parameters by roughly a factor of $1/\sqrt {2}$. Note, however, that when both images are used, $\overline {\boldsymbol {\mu }}$ vanishes automatically (where the bar now indicates weighted integration over both images) because the total number of detected photons is independent of the polarization state. The Fisher information matrix is then given by the simpler form in Eq. (27) without having to impose the constraints introduced in Section 4. While removing these constraints would clearly affect the variational derivation presented earlier, we now show that the results found earlier remain nearly optimal for the case when both components are imaged.

Consider again the representative case of nearly unpolarized light ($\mathbf {s}\approx \textbf {0}$), for which the uncertainty in the estimation of the Stokes parameters is largest, and for which the Fisher information matrix in Eq. (27) reduces to $\mathcal {N}$ times $\overline {\mu _m\mu _n}$. Given the forms found for $\mathcal {I}_n^{(j)}(\mathbf {x})$ (antisymmetric in $x$ and symmetric in $y$ for $n=1$, symmetric in $x$ and antisymmetric in $y$ for $n=2$, and rotationally symmetric for $n=3$), the $3\times 3$ matrix with components $\overline {\mu _m\mu _n}$ is automatically diagonal. Further, because $\boldsymbol {\mu }(\mathbf {x})$ is a unit vector, the trace of this matrix is unity. The uncertainty in the measurement is then minimized if the three diagonal elements have equal magnitude, that is, if $\overline {\mu _n^{2}}=1/3$ for $n=1,2,3$. It turns out that the optimal BM solution derived in Section 7 nearly achieves this, giving $\overline {\mu _1^{2}}=\overline {\mu _2^{2}}=0.32$ and $\overline {\mu _3^{2}}=0.35$, while the SEO gives $\overline {\mu _1^{2}}=\overline {\mu _2^{2}}=0.35$ and $\overline {\mu _3^{2}}=0.30$. The cross-sections of the uncertainty ellipsoids for both the optimal BM solution and the SEO are shown (in blue) in Fig. 6 for the case in which $1500$ photons are detected in each PSF, giving a total number of photons of $\mathcal {N}=3000$. The first thing one notices is that the uncertainties are indeed smaller by a factor of about $1/\sqrt {2}$ with respect to those where only the PSF for the $\mathbf {e}_1$ component was used (red). Also, the uncertainties are now exactly symmetric in $s_3$, as is the coverage of $s_3$ shown in Fig. 7(b).

Finally, note that the relations $\overline {\mu _m\mu _n}=\overline {\mu _n^{2}}\,\delta _{mn}$ and $\overline {\boldsymbol {\mu }_n}=\textbf {0}$, valid for the optimal and SEO solutions (and whether only $\mathbf {e}_1$ or both $\mathbf {e}_1$ and $\mathbf {e}_2$ components are imaged), imply that the four basic PSFs, ${\cal I}_n^{(j)}(\mathbf {x})$, are orthogonal (and almost orthonormal since $\overline {\mu _n^{2}}\approx 1/3$) under the weight $1/{\cal I}_0^{(j)}(\mathbf {x})$. This orthogonality is useful in the retrieval of the parameters from the measured PSFs.

10. Concluding remarks

A variational calculation was performed to determine the spatial birefringence distribution to be placed at the pupil plane of an imaging system that optimizes the encoding of polarization in the PSF’s shape while keeping the PSF as small as possible. This optimal birefringence pattern was found to be similar (both in distribution and performance) to that found naturally near the equilibrium points of a transparent window under stress: the mask is not birefringent at its center, but its retardance grows with the radial pupil coordinate $u$, while the orientation of the slow and fast axes rotates with the azimuthal pupil coordinate $\phi$ as $\pm \phi /2$. (As mentioned earlier, other solutions with equal performance correspond to birefringence distributions that result from the cascading of these solutions with uniform wave plates.) The only differences were in how the retardance grows with $u$ away from the center, but these cause only small variations in the resulting PSFs. If the variational derivation had been based on a different measure of PSF size and/or requirements for polarimetric optimality (e.g. a measure of the uniformity of $\mu _3$ in Fig. 7), the details of the optimal $\delta (u)$ would change slightly, these variations being comparable to those between any of these optimal solutions and the SEO. That is, the mechanics of stress-induced birefringence naturally provide a nearly optimal birefingence distribution for imaging polarimetry, without the need for nano-fabrication or combinations of spatial light modulators: by applying a pressure distribution with trigonal symmetry to the edges of a a glass window, a nearly-optimal birefringence mask (the SEO) results. Note that the SEO was also found to be nearly optimal for generating focused beams with a highly localized zero of intensity [15].

Since the aim of this work was to show that the optimal BM distribution has very similar characteristics and performance to the SEO, we considered the continuous functional form of the PSFs. The performance of these devices in real applications, however, will also depend on the pixelation and general characteristics of the detector, as well as on the algorithms used for the retrieval of the Stokes parameters from the measured PSFs. Preliminary work on this topic can be found in [13]. These aspects are particularly important in applications such as fluorescence microscopy [9], where photons are scarce and the PSFs are used to extract information not only about the two-dimensional polarization studied here, but of the three-dimensional polarization (as well as the three-dimensional position) of the light emitted by a fluorophore. For such applications, the near-optimality of the SEO is of central importance.