## Abstract

Fluorophores that are fixed during image acquisition produce a diffraction pattern that is characteristic of the orientation of the fluorophore’s underlying dipole. Fluorescence localization microscopy techniques such as PALM and STORM achieve super-resolution by applying Gaussian-based fitting algorithms to in-focus images of individual fluorophores; when applied to fixed dipoles, this can lead to a bias in the range of 5–20 nm. We introduce a method for the joint estimation of position and orientation of dipoles, based on the representation of a physically realistic image formation model as a 3-D steerable filter. Our approach relies on a single, defocused acquisition. We establish theoretical, localization-based resolution limits on estimation accuracy using Cramér-Rao bounds, and experimentally show that estimation accuracies of at least 5 nm for position and of at least 2 degrees for orientation can be achieved. Patterns generated by applying the image formation model to estimated position/orientation pairs closely match experimental observations.

©2009 Optical Society of America

## 1. Introduction

Fluorescence localization microscopy (FLM) has emerged as a powerful family of techniques for optically imaging biological samples at molecular resolutions, down to the nanometer scale [1]. This is achieved by employing specific fluorescent labels that can be activated [2], switched on [3], or that intrinsically blink [4], which make it possible to image sparse subsets of the fluorophores contained within the specimen in sequence. Fluorophores appear as spatially isolated spots in the resulting image frames; their center can then be computationally localized with an accuracy that far surpasses the resolution limit formulated by Abbe. A super-resolved image of the specimen is generated by imaging a large number of such fluorophore subsets, localizing every molecule within each frame, and combining the resulting positions to render a composite image at a finer scale. The resolution achieved by these techniques is thus directly dependent upon the achievable localization accuracy. This accuracy in turn depends on the number of photons collected [5] and on the image formation model used in the localization algorithm.

The primary implementations of FLM proposed to date include photoactivated localization microscopy (PALM) [6], fluorescence photoactivation localization microscopy (FPALM) [7], and stochastic optical reconstruction microscopy (STORM) [8]. Developed contemporaneously, these methods demonstrated resolutions at the 10 nm scale experimentally. Due to the long acquisition times required for the collection of a sufficient amount of frames, these methods were initially designed to image 2-D sections of thin, fixed specimens. For such samples, it is generally assumed that the image of a single molecule corresponds to the in-focus section of the system’s 3-D point spread function (PSF), which can be approximated by a 2-D Gaussian function. Under the latter assumption, it has been shown that Gaussian-based fitting algorithms do not result in a significant loss in localization accuracy [9].

Prior to the introduction and practical feasibility (due to a lack of suitable fluorophores) of FLM, localization-based approaches were generally limited to single molecule tracking applications, with inherent isolation of individual fluorophores. In a study of the progression of the molecular motor Myosin V on actin fibers, Yildiz et al. achieved a Gaussian localization-based resolution of 1.5 nm, and accordingly named their technique fluorescence imaging with one-nanometer accuracy (FIONA) [10]. In the framework of single particle tracking, theoretical limits on the achievable localization accuracy have been formulated for lateral (i.e. *x*-*y*) localization [11] and axial localization [12]; these limits translate into a measure of resolution when extended to incorporate multiple sources.

In all of the previously cited methods, it is assumed that individual fluorophores act as isotropically emitting point sources (which also implies that their image corresponds to a section of the system’s PSF). Whereas this is valid for molecules that freely rotate during image acquisition, it does not hold when imaging fixed fluorophores. Their corresponding diffraction patterns differ significantly from the PSF, and are highly specific of the orientation of the fluorophore’s underlying electromagnetic dipole [13, 14].

#### 1.1. Localization of fluorescent dipoles

Dipole diffraction patterns are generally not radially symmetric, and as a consequence, their point of maximum intensity in the image plane is shifted with respect to the position of the fluorophore. Applying Gaussian-based localization to such patterns can lead to a bias in the range of 5–20 nm, even for dipoles that are imaged in focus [15]. This is especially significant in the context of FLM, where resolutions of the same order are striven for. Due to the complexity of the dipole patterns, avoiding this bias requires localization methods based on a more accurate image formation model than a simple Gaussian.

Moreover, localization based on a physically realistic model makes it possible to estimate the dipole’s orientation in addition to its position. So far, this has been exploited in an extension of the FIONA technique named dipole orientation and position imaging (DOPI) [16]. In the study of Myosin V, the estimation of the fluorescent labels’ orientation made it possible to characterize its progression on actin much more precisely than with FIONA. A shortcoming of the current implementation of DOPI is its reliance on two separate procedures for orientation and position estimation.

The current state of the art in dipole orientation and position estimation is based on matched filtering of defocused images. Using an advanced image formation model for dipoles [17], Patra et al. proposed to estimate the orientation based on a finite number of precomputed templates corresponding to rotated versions of a dipole diffraction pattern for a given amount of defocus [18]. The accuracy of this approach is inherently limited by the angular sampling stemming from the finite number of templates, which is usually restricted due to computational cost. As a further consequence, the matched filtering approach limits position estimates to pixel-level accuracy. In DOPI, super-resolved position information is recovered by taking two images of every fluorophore: a defocused image for orientation estimation, and an in-focus image for position estimation, which uses the Gaussian-based localization method proposed in FIONA.

Other methods for estimating the orientation of fluorescent dipoles have been proposed in the literature. Often relying on specific instrumentation, most of these methods are incompatible with a joint position and orientation approach over a reasonably large field of view. Notable techniques are based on direct imaging of the emission pattern in the back focal plane of the objective lens (see, e.g., [19–21]), and on annular illumination [22].

#### 1.2. New approach based on 3-D steerable filters

In this work, we introduce an efficient method for the joint estimation of position and orientation for fluorescent dipoles from a single (defocused) image. We show that image formation for dipoles can be expressed as a linear combination of six templates weighted by trigonometric functions of the dipole angles. As a consequence, the estimation problem can be formulated as the optimization of a 3-D steerable filter [23]. The resulting algorithm serves both as an initial detection of dipoles with localization accuracy at the pixel level, and as an efficient means for iteratively updating the orientation estimates when performing the subsequent sub-pixel localization to achieve super-resolution. We formulate the theoretical limits on estimation accuracy for position and orientation using Cramér-Rao bounds (CRB), and show that our method reaches these bounds. Notably, these limits make it possible to establish the acquisition settings (e.g., defocus) that will lead to the highest estimation accuracy. Our experimental results indicate that dipoles can be localized with a position accuracy of at least 5 nm and an orientation accuracy of at least 2 degrees, at a peak signal-to-noise ratio (PSNR) of approximately 25 dB.

#### 1.3. Organization of the paper

In the next section, we introduce an image formation model for fluorescent dipoles and show how it decomposes into six non-orthogonal templates. Subsequently, in Section 3, we use this steerable decomposition to formulate a localization algorithm for the estimation of both the position and orientation of fluorescent dipoles. In Section 4, we establish the theoretical limits on the accuracy of these estimators, using Fisher information. Finally, in Section 5, we demonstrate our method on simulated and experimental data. We conclude with a discussion of potential applications and extensions of the proposed technique.

## 2. Image formation

A single fluorescent molecule can be characterized as a harmonically oscillating dipole with moment p = (sin *θ _{p}* cos

*ϕ*, sin

_{p}*θ*sin

_{p}*ϕ*, cos

_{p}*θ*) and position

_{p}*x*= (

_{p}*x*,

_{p}*y*,

_{p}*z*), where the angles

_{p}*θ*and

_{p}*ϕ*are the zenith (i.e., between the dipole and the optical axis) and azimuth angle of the dipole, respectively. In our notation, the subscript

_{p}*p*indicates a dipole parameter. The intensity radiated by such a dipole is modeled by propagating the emitted electric field through the optical system. We begin by expressing the amplitude of the electric field in object space as

where the vectors **e**
* _{p}^{s}*, and

**e**

^{s}are the

*p*- and

*s*-polarized components of

*in the sample layer [24], as illustrated in Fig. 1 (the superscript*

**ε**_{s}*s*on the

**e**

_{p}vector denotes the sample medium;

**e**

*remains constant throughout the system). A constant factor representing the magnitude (see, e.g., [24, 25]) has been neglected in this expression and will be reintroduced later. We assume a standard model for the object space, consisting of a sample layer (refractive index*

_{s}*n*), a glass layer (corresponding to the coverslip, index

_{s}*n*), and an immersion layer (index

_{g}*n*) [26]. After propagating through these layers and the objective, the field is given by

_{i}where *t _{p}*

^{(l)}and

*t*

_{s}^{(l)}are the Fresnel transmission coefficients for

*p*-polarized and

*s*-polarized light from layer

*l*to layer

*l*+ 1, i.e.,

$${t}_{s}^{\left(l\right)}=\frac{2{n}_{l}\mathrm{cos}{\theta}_{l}}{{n}_{l}\mathrm{cos}{\theta}_{l}+{n}_{l+1}\mathrm{cos}{\theta}_{l+1}},$$

where (*n*
_{1},*n*
_{2},*n*
_{3}) = (*n _{s}*,

*n*,

_{g}*n*) for the configuration considered here. Note that Eq. (2) can easily be generalized to an arbitrary number of strata. The unit vectors of the relevant polarization directions, expressed in spherical coordinates, are

_{i}$$=(\frac{1}{{n}^{s}}\sqrt{{n}_{s}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}\mathrm{cos}\varphi ,\frac{1}{{n}^{s}}\sqrt{{n}_{s}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}\mathrm{sin}\varphi ,-\frac{{n}_{i}}{{n}_{s}}\mathrm{sin}\theta )$$

$${\mathbf{e}}_{p}^{a}=\left(\mathrm{cos}\varphi ,\mathrm{sin}\varphi ,0\right)$$

$${\mathbf{e}}_{s}=\left(-\mathrm{sin}\varphi ,\mathrm{cos}\varphi ,0\right),$$

where *θ _{s}* is the zenith angle between the wavevector

**k**and the optical axis in the sample layer,

*θ*is the corresponding angle in the immersion layer, and

*ϕ*is the azimuth angle (see Fig. 1). The wavevector is thus parameterized as

where $k=\frac{2\pi}{\lambda}$ is the wavenumber. The field in the detector plane, expressed in object space coordinates, is given by the Richards-Wolf integral [27, 28]:

$$=\frac{-i{A}_{0}}{\lambda}{\int}_{0}^{2\pi}{\int}_{0}^{\alpha}{\epsilon}_{a}{e}^{\mathrm{ikr}{n}_{i}\mathrm{sin}\phantom{\rule{.2em}{0ex}}\theta \phantom{\rule{.2em}{0ex}}\mathrm{cos}(\varphi -{\varphi}_{d}){e}^{-ik{n}_{i}z\mathrm{cos}\theta}}{e}^{\mathrm{ik}\Lambda (\theta ;\mathit{\tau})}\sqrt{\mathrm{cos}\phantom{\rule{.2em}{0ex}}\theta}\mathrm{sin}\phantom{\rule{.2em}{0ex}}\theta \phantom{\rule{.2em}{0ex}}\mathrm{d\theta}\phantom{\rule{.2em}{0ex}}d\varphi ,$$

where *A*
_{0} is the scalar amplitude of the field, *λ* its wavelength, and a is the maximal angular aperture of the objective (i.e., *α* = sin^{-1} (NA/*n _{i}*)). The phase component Λ(

*θ*;

*τ*) describes the system’s aberrations. Its argument

is a vector containing the optical parameters of the setup; the subscripts on the refractive indices *n* and thicknesses *t* indicate the respective layer for each parameter (see Fig. 1), and an asterisk signals a design value (i.e., *n _{i}* is the refractive index of the immersion medium under experimental conditions, and

*n*

_{i}^{*}its corresponding design value). The field

*is a function of the point of observation*

**ε***x*= (-

*r*cos

*ϕ*, -

_{d}*r*sin

*ϕ*,

_{d}*z*) on the detector, where $r=\sqrt{{\left(x-{x}_{p}\right)}^{2}+{\left(y-{y}_{p}\right)}^{2}}$ and

*ϕ*= tan

_{d}^{-1}((

*y*-

*y*)/(

_{p}*x*-

*x*));

_{p}*z*denotes defocus. All units are expressed in object space coordinates, which allows us to place the origin of the coordinate system at the interface between the sample layer and the coverslip in order to facilitate the expression of the phase term Λ(

*θ*;

*τ*).

We describe the system’s aberrations based on the optical path difference (OPD) between actual (experimental) imaging conditions and the corresponding design values of the system, as proposed in [26]. Aberrations are assumed to arise exclusively as the result of a mismatch between the refractive indices and/or the thicknesses of the different layers in the sample setup. State-of-the-art vectorial PSF calculations all incorporate this phase term (see, e.g., [29, 30]), and although the formalism employed differs from [26], these approaches lead to the same result for the phase and are considered equivalent [31, 32]. We thus replace the aberration term Λ(*θ*; *τ*) with the OPD Λ(*θ*, *z*; *z _{p}*,

*τ*), defined as

$$+{z}_{p}\sqrt{{n}_{s}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}+{t}_{g}\sqrt{{n}_{g}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}$$

$$-{t}_{g}^{*}\sqrt{{n}_{{g}^{*}}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}-{t}_{g}^{*}\sqrt{{n}_{{i}^{*}}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta .}$$

Note that Λ(*θ*,*z*;0,*τ*) is the phase term corresponding to the standard defocus model [33].

By substituting Eq. (2) into Eq. (6) and simplifying the result, we can rewrite the vector * ε* as

where

$${I}_{1}\left(\mathit{x};{\mathit{x}}_{p},\mathit{\tau}\right){\text{}=\int}_{0}^{\alpha}{B}_{1}\left(\theta \right){t}_{p}^{\left(1\right)}{t}_{p}^{\left(2\right)}\frac{{n}_{i}}{{n}_{s}}\mathrm{sin}\theta d\theta $$

$${I}_{2}\left(\mathit{x};{\mathit{x}}_{p},\mathit{\tau}\right)={\int}_{0}^{\alpha}{B}_{2}\left(\theta \right)\left({t}_{s}^{\left(1\right)}{t}_{s}^{\left(2\right)}-{t}_{p}^{\left(1\right)}{t}_{p}^{\left(2\right)}\frac{1}{{n}_{s}}\sqrt{{n}_{s}^{2}-{n}_{i}^{2}{\mathrm{sin}}^{2}\theta}\right)d\theta $$

with

It should be noted that these integrals are standard for vectorial PSF calculations (for equivalent expressions see, e.g. [30, 32]). The amplitude factor *A*
_{0} was dropped from this formulation and will again be reintroduced later. Due to an already involved notation, we choose to omit the argument (*x*; *x _{p}*,

*τ*) from functions such as ∣

*I*

_{0}∣

^{2}from this point on.

The orientation-dependent intensity in the detector plane is accordingly given by

$$\phantom{\rule{5.2em}{0ex}}={\mathrm{sin}}^{2}{\theta}_{p}\left({\mid {I}_{0}\mid}^{2}+{\mid {I}_{2}\mid}^{2}+2\mathrm{cos}\left(2{\varphi}_{p}-2{\varphi}_{d}\right)\Re e\left\{{I}_{0}^{*}{I}_{2}\right\}\right)$$

$$\phantom{\rule{6.2em}{0ex}}-2\mathrm{sin}\left(2{\theta}_{p}\right)\mathrm{cos}\left({\varphi}_{p}-{\varphi}_{d}\right)\Im m\left\{{I}_{1}^{*}\left({I}_{0}+{I}_{2}\right)\right\}+4{\mid {I}_{1}\mid}^{2}{\mathrm{cos}}^{2}{\theta}_{p}$$

$$\phantom{\rule{5.2em}{0ex}}={p}^{T}\mathrm{Mp}.$$

An asterisk denotes complex conjugation, and *v*
^{T} stands for the Hermitian transpose of the vector *ν*; ℜ*e*{*ν*} and ℑ*m*{*ν*} represent the real and imaginary components of *ν*, respectively. The symmetric matrix M = (*m _{ij}*)1≤

*i,j*≤3 is specified by

$${m}_{12}=2\Re e\left\{{I}_{0}^{*}{I}_{2}\right\}\mathrm{sin}2{\varphi}_{d}$$

$${m}_{13}=-2\mathrm{cos}{\varphi}_{d}\Im m\left\{{I}_{1}^{*}\left({I}_{0}+{I}_{2}\right)\right\}$$

$${m}_{22}={\mid {I}_{0}\mid}^{2}+{\mid {I}_{2}\mid}^{2}+2\Re e\left\{{I}_{0}^{*}{I}_{2}\right\}\mathrm{cos}2{\varphi}_{d}$$

$${m}_{23}=-2\mathrm{sin}{\varphi}_{d}\Im m\left\{{I}_{1}^{*}\left({I}_{0}+{I}_{2}\right)\right\}$$

$${m}_{33}=4{\mid {I}_{1}\mid}^{2}.$$

The quadratic form in Eq. (12) is of particular interest, since it decouples the dipole orientation from the calculation of the field propagation. The dipole diffraction pattern can thus be modeled as the linear combination of six non-orthogonal templates, which forms the basis of the steerable filter-based algorithm presented in this work.

The result of Eq. (12) is consistent with the other vectorial PSF models cited earlier [28–32]: for a fluorophore that is freely rotating during exposure, the resulting intensity is obtained by integrating the above result over all possible orientations, i.e.,

$$=\frac{8\pi}{3}\left({\mid {I}_{0}\mid}^{2}+2{\mid {I}_{1}\mid}^{2}+{\mid {I}_{2}\mid}^{2}\right).$$

An illustration of the different diffraction patterns observed for different orientation and defocus values is provided in Fig. 11 in the appendix.

#### 2.1. Noise model

In single molecule fluorescence microscopy, shot noise is the dominant source of noise. Further potential sources include a background term due to autofluorescence, as well as residual signals from other fluorophores in the sample. Read-out noise in the detector may also contribute to observations. Whereas the former sources obey Poisson statistics, read-out noise is Gaussian distributed, which presumes an additive noise model. However, given that the Poisson distribution rapidly converges towards a Gaussian with equal mean and variance when the variance is large enough (this is usually considered the case when *σ*
^{2} > 10), we propose a general noise model consisting of a shifted Poisson formulation that incorporates a term accounting for the read-out noise factor and the background. Consequently, we formulate the expected photon count *q*̅(*x*; *x _{p}*,

*τ*) corresponding to a point x on the detector (in object-space coordinates) as

where *A* is the amplitude, *c* is a conversion factor, and *b* is the sum of the background fluorescence signal and the variance *σ _{r}*

^{2}, (in intensity) of the read-out noise. The probability of detecting

*q*photons at

*x*is then given by

#### 2.2. Pixelation

The pixelation of the detector has a non-negligible effect on the measured intensity distribution corresponding to a dipole, and must therefore be taken into account in the implementation of the model. Hereinafter, we assume that whenever a point of observation *x* represents a pixel on the CCD, functions of *x* incorporate integration over the pixel’s area (pixels are assumed to be contiguous and non-overlapping).

## 3. Dipole localization using steerable filters

A standard solution to estimating the position and orientation of arbitrarily rotated image features consists in correlating the input data with a set of filters corresponding to rotated versions of the feature template. This is the approach currently used in joint dipole orientation and position estimation methods [16, 18]. Mathematically, it is formulated as

where *θ _{p}*

^{*}(

*x*) and

*ϕ*

_{p}^{*}(

*x*) are the optimal orientations in every point

*x*,

*f*(

*x*) is the measured data, * is the convolution operator, and

*g*

_{θp,ϕp}(

*x*) =

*h*

_{θp,ϕp}(

*x*; (0,0,

*z*),

_{p}*τ*) is a feature template corresponding to the dipole diffraction pattern for the orientation (

*θ*,

_{p}*ϕ*). The accuracy on

_{p}*θ*

_{p}^{*}(

*x*) and

*ϕ*

_{p}^{*}(

*x*) depends on the sampling of

*θ*and

_{p}*ϕ*chosen for the optimization—for accuracies of a few degrees, this requires convolution with hundreds of templates, rendering the process very costly.

_{p}For specific feature template functions, rotation of the templates can be decoupled from the convolution by decomposing the template into a small number of basis templates that are interpolated by functions of the orientation. These functions, called steerable filters, were introduced by Freeman and Adelson [23]. The best known class of steerable filters are Gaussian derivatives, due to their applicability to edge and ridge detection in image processing [34].

The quadratic form of Eq. (12) shows that the 3-D rotation of a dipole can be decoupled from filtering; i.e., that the orientation and position estimation can be expressed through a 3-D steerable filter. Specifically, the dipole diffraction pattern is decomposed into six non-orthogonal templates, as illustrated in Fig. 2. Orientation estimation then amounts to filtering measured dipole patterns with these templates, followed by optimizing the trigonometric weighting functions that interpolate the dipole model. This process is illustrated in the filterbank representation of Fig. 3.

For each spatial location *x*, the dipole model is fitted to the data by minimizing the least-squares criterion

$$={\parallel A{h}_{{\theta}_{p}{,\varphi}_{p}}\left(\mathit{x};{\mathit{x}}_{p},\mathit{\tau}\right)\parallel}^{2}+{\int}_{\Omega}f{\left(\mathit{x}-\mathit{v}\right)}^{2}d\mathit{v}$$

$$-2A{h}_{{\theta}_{p}{,\varphi}_{p}}\left(\mathit{x};{\mathit{x}}_{p},\mathit{\tau}\right)*f\left(\mathit{x}\right),$$

where *f*(*x*) is the observed signal from which the average background value has been subtracted, and where Ω ⊂ ℜ^{2} is the support of *h*
_{θp,ϕp}. The correlation term is steerable and can be expressed as

$$={\mathbf{p}}^{T}{\mathbf{M}}_{f}\mathbf{p},$$

where *a _{ij}*(

*θ*,

_{p}*ϕ*) are the weighting functions given in Fig. 3, and where [

_{p}*]*

**M**_{f}_{ij}=

*m**

_{ij}*f*. Concretely, this means that we can compute Eq. (19) for each location

*x*very efficiently by first convolving the image

*f*with the six templates

*m*—which yields

_{ij}*—followed by applying the trigonometric weights*

**M**_{f}*a*(

_{ij}*θ*,

_{p}*ϕ*), which amounts to computing the quadratic form

_{p}**p**

^{T}

**M**

_{f}

**p**. Due to the symmetry of

**M**, only six basis templates are involved in the process.

In order to simplify the expression for the optimization of Eq. (18), we rewrite the model energy term, which is independent of *ϕ _{p}*, as

where **u**
*θ _{p}* = (sin

^{2}

*θ*, sin 2

_{p}*θ*, cos

_{p}^{2}

*θ*)

_{p}^{T}, and where

**E**is defined as

The notation 〈*m _{ij}*〉 stands for integration over the support Ω of

*m*(

_{ij}*x*). We can then rewrite the cost criterion as

The data term ∫_{Ω}
*f*(*x*-*ν*)^{2} d*ν* in Eq. (18) has no effect on the optimization and was eliminated from the criterion.

#### 3.1. Orientation estimation

The criterion in Eq. (22) cannot be solved in closed form for *θ _{p}*,

*ϕ*, and A. However, for a given set of templates corresponding to a position estimate, an iterative algorithm is obtained by setting the partial derivatives

_{p}$$\frac{\partial}{\partial {\varphi}_{p}}J\left(\mathit{x};{\theta}_{p},{\varphi}_{p}\right)=-4A{\mathbf{p}}^{T}{\mathbf{M}}_{f}\frac{\partial}{\partial {\varphi}_{p}}\mathbf{p}$$

to zero, and alternately solving for *θ _{p}* and

*ϕ*; we found this to be more efficient than a standard gradient-descent based approach. The quartic equations whose solution yields the optimal values for both angles are given in the appendix. The notation $\frac{\partial}{\partial t}\nu $ stands for the componentwise derivative of the vector

_{p}*ν*with respect to

*t*. Between these iterations, the amplitude is updated using the least-squares estimator

For fixed values of *x _{p}* and

*z*, the angles can be estimated to sufficient accuracy in a small number of iterations (usually less than five).

#### 3.2. Sub-pixel position estimation

Discretization of Eq. (22) yields a filter-based algorithm that returns position estimates with pixel-level accuracy. In order to recover super-resolved position estimates, an iterative fitting algorithm analogous to the Gaussian-based fitting algorithms employed by FLM can be applied to the dipole model in Eq. (12). In this work, this was achieved with the Levenberg-Marquardt algorithm. An alternative would be to use a maximum-likelihood-based cost criterion, which is optimal with respect to the Poisson-based noise model [12]. However, the localization accuracy achieved using the proposed least-squares criterion is near-optimal in the sense that it reaches the theoretical limits (see Section 4).

#### 3.3. Joint estimation algorithm

In addition to the lateral position (*x _{p}*,

*y*), amplitude

_{p}*A*, and orientation (

*θ*,

_{p}*ϕ*) of dipoles, which can be estimated using the proposed steerable filter in the minimization of Eq. (22), the axial position

_{p}*z*of the dipole and the position of the focus

_{p}*z*are further degrees of freedom in the model. Estimation of these two parameters can be challenging due to their mutual influence in the phase term of Eq. (8); possible solutions for resolving this ambiguity were discussed in [12], in the framework of 3-D localization of isotropic point sources.

However, for the distances *z _{p}* from the coverslip that are observable under TIRF (total internal reflection fluorescence) excitation, which is currently the method of choice for FLM due to its thin depth excitation (around 100 nm), the axial shift variance in the PSF is essentially negligible. This means that a dipole with axial position

*z*that is observed at the focal distance

_{p}*z*is virtually undistinguishable from an identical dipole at the origin that is observed at

*z*-

*z*. Consequently, we set

_{p}*z*= 0 for our experiments.

_{p}For a given value of defocus *z*, the minimization of Eq. (22) is achieved in two steps. The input image is first filtered with the six templates corresponding to *z*, which yields orientation and position estimates that are accurate at the pixel level. This is followed by an iterative optimization that refines the position estimates to sub-pixel accuracy. At each iteration, a new set of templates corresponding to the sub-pixel position estimates is evaluated; the corresponding orientation is computed with the method described in Section 3.1. If required, the estimation of *z* can be included between steps of this iterative procedure, both at the pixel level and the super-resolution level. A simplified flowchart representation of the algorithm is given in Fig. 4.

## 4. Dipole localization accuracy

Resolution in FLM is defined as a function of localization accuracy. Theoretical limits on localization accuracy can be established based on the Cramér-Rao bounds for the parameters of interest (see [11, 35] for an analysis of lateral localization accuracy, and [12] for axial localization accuracy). The CRB yields a lower bound on the variance of an estimator, provided that the latter is unbiased. For multi-parameter estimation, this bound is obtained by computing the inverse of the Fisher information matrix [36]. For our image formation model, this matrix is defined through

where the parameters to be estimated are *ϑ* = (*x _{p}*,

*y*,

_{p}*θ*,

_{p}*ϕ*,

_{p}*z*,

*A*). In the matrix

**F**, the cross-terms between

*ϕ*and the other parameters are zero, which slightly simplifies the inversion. The CRBs for

_{p}*ϑ*= (

*x*,

_{p}*y*,

_{p}*θ*,

_{p}*z*,

*A*) are thus given by

and the bound for *ϕ _{p}* by

The partial derivatives of Eq. (12) relevant to the computation of the CRBs are given in the appendix.

An example of the CRBs for all parameters is shown in Fig. 5, as a function of *z* and *θ _{p}*. As the different plots illustrate, the bounds are relatively flat for a large range of values. In-focus imaging leads to lower accuracy for orientation and defocus estimation. Off-focus imaging, on the other hand, preserves excellent localization accuracy in the plane, while leading to higher accuracies for orientation and defocus estimation. Note that these observations hold true for the general case; changing the system and acquisition parameters essentially amounts to a scaling of the bounds.

#### 4.1. Performance of the algorithm

The localization accuracy of the proposed algorithm was evaluated by performing the fit on simulated acquisitions generated using the image formation model in Eq. (15). In Fig. 6, we show standard deviations of the estimation results for each parameter compared to the associated CRB. The standard deviations consistently match the value of the CRB, which indicates that the performance of the algorithm is near-optimal. A fully optimal solution for the proposed Poisson-based noise model would consist in a maximum-likelihood formulation. However, this would lead to a computationally less efficient solution compared to the filter-based approach obtained through the proposed least-squares formulation.

## 5. Results

#### 5.1. Experimental setup

To demonstrate the performance of our algorithm experimentally, we imaged single Cy5 dipoles at an air/glass interface. Samples were prepared by drying a solution of Cy5 molecules bound to single-stranded DNA (ssDNA) fragments onto a coverslip (isolated Cy5 binds to the glass surface at a specific orientation, leading to limited diversity in the observed diffraction patterns). These were then imaged using TIRF excitation at a wavelength of *λ* = 635 nm. Imaging was performed on the microscope setup described in [37] (*α* Plan-Fluar 100 × 1.45 NA objective with Immersol^{™} 518F (Carl Zeiss Jena, Jena, Germany)), where we added a highly sensitive CCD camera (Luca^{EM} S 658M (Andor, Belfast, Northern Ireland), 10 × 10 *μ*m^{2} pixels) in the image plane. Images were acquired with an EM gain setting of 100 (the shot noise variance distortion resulting from the electron multiplying process [38] can be compensated in the background term of Eq. (15) in order to conserve the Poisson statistics). Despite weak photostability in air, individual molecules could be imaged for several tens of seconds before bleaching.

#### 5.2. Evaluation

For an experimental assessment of the accuracy of our algorithm, we applied the estimation to image sequences of individual molecules taken under identical acquisition conditions. Standard deviations on position, orientation, and defocus for such an experiment are given in Table 1; some frames from the sequence used are shown in Fig. 7. The results are in good agreement with the values predicted by the CRB. Specifically, it is the relative magnitude between the different standard deviations that precisely matches the behavior predicted by the CRB (for comparison, see Fig. 6). Although the average PSNR was relatively high for these acquisitions (26.43 dB, computed with respect to the average of all acquisitions), they provide an indication of experimentally achievable accuracy.

In a similar experiment, we compared the estimation accuracy over a series of frames of the same dipole taken at different levels of defocus. The resulting values, given in Table 2, indicate that the algorithm and image formation model are consistent and correspond well to the observed measurements. The frames from this experiment are shown in Fig. 8.

The performance of the steerable filter-based estimation of orientation and position is illustrated in Fig. 9 for simulated data and in Fig. 10 for an experimental acquisition. The image generated using the estimated positions and orientations matches the experimental measurement well, despite some residual aberrations that slightly modify the observed diffraction patterns.

## 6. Discussion

The method proposed in this paper should be of immediate benefit to optical imaging-based studies of molecular motors, as demonstrated in the DOPI paper by Toprak et al. [16]. These authors showed that the diffraction pattern of dipolar quantum dots [39] attached to myosin V can be clearly detected and tracked over time. Due to the computational cost of the matched filtering algorithm employed [18], the angular sampling of the templates was limited to 10° for *ϕ _{p}* and to 15° for

*θ*, which also required the selection of dipoles oriented almost orthogonal to the optical axis for better accuracy (see supporting text of [16]). Due to a relatively small change in the observed diffraction pattern for values of

_{p}*θ*between 0° and 45°, it was assumed that

_{p}*θ*could only be determined with low accuracy using a matched filtering approach [18]. Our theoretical analysis shows that both angles in the parameterization of the dipole’s orientation can be recovered with high accuracy.

_{p}DOPI requires acquisition of a pair of images: a defocused image for orientation estimation, and an in-focus image for Gaussian-based localization. As shown by Enderlein et al. [15], Gaussian-based localization applied to dipole diffraction patterns can introduce a significant bias in the estimated position, even for molecules that are imaged in focus. Our method removes these limitations: it performs well over a wide range of dipole orientations (only the estimation of *φ _{p}* for values of

*θ*that are close to zero remains inherently inaccurate), and is relatively insensitive to the amount of defocus used, as shown in Fig. 5.

_{p}Recent advances in the development of FLM include expansions to 3-D, using either a combination of cylindrical optics and bivariate Gaussian fitting [40], or an interferometric estimation of fluorophore positions in direction of the optical axis [41], as well as orientation sensitive imaging using polarization-based techniques [42]. All of these approaches require some modification of the imaging system; the fitting algorithm described in this paper can be readily adapted to include 3-D localization, and is compatible with standard TIRF setups.

## 7. Conclusion

We have introduced an efficient algorithmic framework based on 3-D steerable filters for the estimation of fluorescent dipole positions and orientations. Image formation for fluorescent dipoles can be expressed as a function of six templates weighted by functions of the dipole’s orientation, leading to an effective filter-based estimation at the pixel-level, followed by an iterative refinement that yields super-resolved estimates. Experimental results on Cy5 dipoles demonstrate the potential of the proposed approach; estimation accuracies of the order of 5 nm for position and 2° for orientation are reported.

## A. Appendix

#### A.1. Quartic equations for orientation estimation

The orientation parameters *θ _{p}* and

*ϕ*are estimated by iteratively solving the two sets of quartic equations given below. The solution for

_{p}*θ*is obtained by solving

_{p}$$+{\mathrm{tan}}^{3}{\theta}_{p}\left({f}_{33}-({f}_{11}\mathrm{cos}{\left({\varphi}_{p}\right)}^{2}+{f}_{12}\mathrm{sin}\left(2{\varphi}_{p}\right)+{f}_{22}\mathrm{sin}{\left({\varphi}_{p}\right)}^{2})+A\left({e}_{11}-{e}_{13}-2{e}_{22}\right)\right)$$

$$+3{\mathrm{tan}}^{2}{\theta}_{p}A\left({e}_{12}-{e}_{23}\right)$$

$$+\mathrm{tan}{\theta}_{p}\left({f}_{33}-({f}_{11}\mathrm{cos}{\left({\varphi}_{p}\right)}^{2}+{f}_{12}\mathrm{sin}\left(2{\varphi}_{p}\right)+{f}_{22}\mathrm{sin}{\left({\varphi}_{p}\right)}^{2})+A\left({e}_{13}-{e}_{33}+2{e}_{22}\right)\right)$$

$$-\left({f}_{13}\mathrm{cos}\left({\varphi}_{p}\right)+{f}_{23}\mathrm{sin}\left({\varphi}_{p}\right)-A{e}_{23}\right)$$

$$=0$$

for tan *θ _{p}*, where [

**E**]

_{ij}=

*e*. Similarly, the solution for

_{ij}*ϕ*is obtained by solving

_{p}$$+{\mathrm{tan}}^{3}{\varphi}_{p}\left({f}_{13}{f}_{23}{\mathrm{cos}}^{2}{\theta}_{p}+{f}_{12}\left({m}_{11}-{f}_{22}\right){\mathrm{sin}}^{2}{\theta}_{p}\right)$$

$$+{\mathrm{tan}}^{2}{\varphi}_{p}\left(\left({\left({f}_{11}-{f}_{22}\right)}^{2}-2{f}_{12}^{2}\right){\mathrm{sin}}^{2}{\theta}_{p}-\left({f}_{13}^{2}+{f}_{23}^{2}\right){\mathrm{cos}}^{2}{\theta}_{p}\right)$$

$$+\mathrm{tan}{\varphi}_{p}\left({f}_{13}{f}_{23}{\mathrm{cos}}^{2}{\theta}_{p}-{f}_{12}\left({f}_{11}-{f}_{22}\right){\mathrm{sin}}^{2}{\theta}_{p}\right)$$

$$+{f}_{12}^{2}{\mathrm{sin}}^{2}{\theta}_{p}-{f}_{23}^{2}{\mathrm{cos}}^{2}{\theta}_{p}$$

$$=0$$

for tan *ϕ _{p}*, where [

**M**

_{f}]

_{ij}=

*f*.

_{ij}#### A.2. Derivatives of the dipole model

The partial derivatives required in the evaluation of the Fisher information and in the fitting algorithm are

$$\phantom{\rule{3.2em}{0ex}}-2\mathrm{sin}\left(2{\theta}_{p}\right)\mathrm{cos}\left({\varphi}_{p}-{\varphi}_{d}\right)\Im m\left\{{I}_{1}^{*}\left(\frac{\partial {I}_{0}}{\partial {x}_{p}}+\frac{\partial {I}_{2}}{\partial {x}_{p}}\right)-\frac{\partial {I}_{1}}{\partial {x}_{p}}\left({I}_{0}^{*}+{I}_{2}^{*}\right)\right\}$$

$$\phantom{\rule{3.2em}{0ex}}+8{\phantom{\rule{.2em}{0ex}}\mathrm{cos}}^{2}{\theta}_{p}\Re e\left\{{I}_{1}^{*}\frac{\partial {I}_{1}}{\partial {x}_{p}}\right\}$$

$$\frac{\partial {h}_{{\theta}_{p},{\varphi}_{p}}}{\partial {\theta}_{p}}=\mathrm{sin}2{\theta}_{p}({\mid {I}_{0}\mid}^{2}+{\mid {I}_{2}\mid}^{2}+2\mathrm{cos}\left(2{\varphi}_{p}-2{\varphi}_{d}\right)\Re e\left\{{I}_{0}^{*}{I}_{2}\right\}-4{\mid {I}_{1}\mid}^{2})\text{}$$

$$\phantom{\rule{3.2em}{0ex}}-4\mathrm{cos}\left(2{\theta}_{p}\right)\mathrm{cos}\left({\varphi}_{p}-{\varphi}_{d}\right)\Im m\left\{{I}_{1}^{*}\left({I}_{0}+{I}_{2}\right)\right\}$$

$$\frac{\partial {h}_{{\theta}_{p},{\varphi}_{p}}}{\partial {\varphi}_{p}}=-4{\mathrm{sin}}^{2}{\theta}_{p}\mathrm{sin}\left(2{\varphi}_{p}-2{\varphi}_{d}\right)\Re e\left\{{I}_{0}^{*}{I}_{2}\right\}$$

$$\phantom{\rule{3.2em}{0ex}}+2\mathrm{sin}\left(2{\theta}_{p}\right)\mathrm{sin}\left({\varphi}_{p}-{\varphi}_{d}\right)\Im m\left\{{I}_{1}^{*}\left({I}_{0}+{I}_{2}\right)\right\}.$$

## Acknowledgments

This work was supported by the Swiss National Science Foundation under grant 200020-109415, as well as by the Hasler foundation under grant 2033.

## References and links

**1. **S. W. Hell, “Microscopy and its focal switch,” Nature Methods **6**, 24–32 (2009). [CrossRef] [PubMed]

**2. **G. H. Patterson and J. Lippincott-Schwartz, “A photoactivatable GFP for selective photolabeling of proteins and cells,” Science **297**, 1873–1877 (2002). [CrossRef] [PubMed]

**3. **M. Bates, T. R. Blosser, and X. Zhuang ,“Short-range spectroscopic ruler based on a single-molecule optical switch,” Phys. Rev. Lett. **94**, 108101(1–4) (2005). [CrossRef] [PubMed]

**4. **K. A. Lidke, B. Rieger, T. M. Jovin, and R. Heintzmann, “Superresolution by localization of quantum dots using blinking statistics,” Opt. Express **13**, 7052–7062 (2005). [CrossRef] [PubMed]

**5. **R. E. Thompson, D. R. Larson, and W. W. Webb, “Precise nanometer localization analysis for individual fluorescent probes,” Biophys. J. **82**, 2775–2783 (2002). [CrossRef] [PubMed]

**6. **E. Betzig, G. H. Patterson, R. Sougrat, O. W. Lindwasser, S. Olenych, J. S. Bonifacio, M. W. Davidson, J. Lippincott-Schwartz, and H. F. Hess, “Imaging intracellular fluorescent proteins at nanometer resolution,” Science **313**, 1642–1645 (2006). [CrossRef] [PubMed]

**7. **S. T. Hess, T. P. K. Girirajan, and M. D. Mason, “Ultra-high resolution imaging by fluorescence photoactivation localization microscopy,” Biophys. J. **91**, 4258–4272 (2006). [CrossRef] [PubMed]

**8. **M. J. Rust, M. Bates, and X. Zhuang, “Sub-diffraction-limit imaging by stochastic optical reconstruction microscopy (STORM),” Nature Methods **3**, 793–795 (2006). [CrossRef] [PubMed]

**9. **M. K. Cheezum, W. F. Walker, and W. H. Guilford, “Quantitative comparison of algorithms for tracking single fluorescent particles,” Biophys. J. **81**, 2378–2388 (2001). [CrossRef] [PubMed]

**10. **A. Yildiz, J. N. Forkey, S. A. McKinney, H. Taekjip, Y. E. Goldman, and P. R. Selvin, “Myosin V walks handover-hand: single fluorophore imaging with 1.5-nm localization,” Science **300**, 2061–2065 (2003). [CrossRef] [PubMed]

**11. **R. J. Ober, S. Ram, and S. Ward, “Localization accuracy in single-molecule microscopy,” Biophys. J. **86**, 1185–1200 (2004). [CrossRef] [PubMed]

**12. **F. Aguet, D. Van De Ville, and M. Unser, “A maximum-likelihood formalism for sub-resolution axial localization of fluorescent nanoparticles,” Opt. Express **13**, 10,503–10,522 (2005). [CrossRef]

**13. **A. P. Bartko and R. M. Dickson, “Imaging three-dimensional single molecule orientations,” J. Phys. Chem. B **103**, 11,237–11,241 (1999).

**14. **R. Schuster, M. Barth, A. Gruber, and F. Cichos, “Defocused wide field fluorescence imaging of single CdSe/ZnS quantum dots,” Chem. Phys. Lett. **413**, 280–283 (2005). [CrossRef]

**15. **J. Enderlein, E. Toprak, and P. R. Selvin, “Polarization effect on position accuracy of fluorophore localization,” Opt. Express **14**, 8111–8120 (2006). [CrossRef] [PubMed]

**16. **E. Toprak, J. Enderlein, S. Syed, S. A. McKinney, R. G. Petschek, T. Ha, Y. E. Goldman, and P. R. Selvin, “Defocused orientation and position imaging (DOPI) of myosin V,” Proc. Natl. Acad. Sci. USA **103**, 6495–6499 (2006). [CrossRef] [PubMed]

**17. **M. Böhmer and J. Enderlein, “Orientation imaging of single molecules by wide-field epifluorescence microscopy,” J. Opt. Soc. Am. A **20**, 554–559 (2003). [CrossRef]

**18. **D. P. Patra, I. Gregor, and J. Enderlein, “Image analysis of defocused single-molecule images for three-dimensional molecule orientation studies,” J. Phys. Chem. A **108**, 6836–6841 (2004). [CrossRef]

**19. **M. A. Lieb, J. M. Zavislan, and L. Novotny, “Single-molecule orientations determined by emission pattern imaging,” J. Opt. Soc. Am. B **21**, 1210–1215 (2004). [CrossRef]

**20. **Z. Sikorski and L. M. Davis, “Engineering the collected field for single-molecule orientation determination,” Opt. Express **16**, 3660–3673 (2008). [CrossRef] [PubMed]

**21. **M. R. Foreman, C. M. Romero, and P. Török, “Determination of the three-dimensional orientation of single molecules,” Opt. Lett. **33**, 1020–1022 (2008). [CrossRef] [PubMed]

**22. **B. Sick, B. Hecht, and L. Novotny, “Orientational imaging of single molecules by annular illumination,” Phys. Rev. Lett. **85**, 4482–4485 (2000). [CrossRef] [PubMed]

**23. **W. T. Freeman and E. H. Adelson, “The design and use of steerable filters,” IEEE Trans. Pattern Anal. Mach. Intell. **13**, 891–906 (1991). [CrossRef]

**24. **E. H. Hellen and D. Axelrod, “Fluorescence emission at dielectric and metal-film interfaces,” J. Opt. Soc. Am. B **4**, 337–350 (1987). [CrossRef]

**25. **G. W. Ford and W. H. Weber, “Electromagnetic interactions of molecules with metal surfaces,” Phys. Rep. **113**, 195–287 (1984). [CrossRef]

**26. **S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oil-immersion objective lens used in three-dimensional light microscopy,” J. Opt. Soc. Am. A **8**, 1601–1613 (1991). [CrossRef]

**27. **E. Wolf, “Electromagnetic diffraction in optical systems—I. An integral representation of the image field,” Proc. R. Soc. London A **253**, 349–357 (1959). [CrossRef]

**28. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems—II. Structure of the image field in an aplanatic system,” Proc. R. Soc. London A **253**, 358–379 (1959). [CrossRef]

**29. **S. W. Hell, G. Reiner, C. Cremer, and E. H. K. Stelzer, “Aberrations in confocal fluorescence microscopy induced by mismatches in refractive index,” J. Microsc. **169**, 391–405 (1993). [CrossRef]

**30. **P. Török and R. Varga, “Electromagnetic diffraction of light focused through a stratified medium,” Appl. Opt. **36**, 2305–2312 (1997). [CrossRef] [PubMed]

**31. **O. Haeberlé, “Focusing of light through a stratified medium: a practical approach for computing microscope point spread functions. Part I: Conventional microscopy,” Opt. Commun. **216**, 55–63 (2003). [CrossRef]

**32. **A. Egner and S. W. Hell, “Equivalence of the Huygens-Fresnel and Debye approach for the calculation of high aperture point-spread functions in the presence of refractive index mismatch,” J. Microsc. **193**, 244–249 (1999). [CrossRef]

**33. **M. Born and E. Wolf, *Principles of Optics*, 7th ed. (Cambridge University Press, 1959).

**34. **M. Jacob and M. Unser, “Design of steerable filters for feature detection using Canny-like criteria,” IEEE Trans. Pattern Anal. Mach. Intell. **26**, 1007–1019 (2004). [CrossRef]

**35. **K. A. Winick, “Cramér-Rao lower bounds on the performance of charge-coupled-device optical position estimators,” J. Opt. Soc. Am. A **3**, 1809–1815 (1986). [CrossRef]

**36. **D. L. Snyder and M. I. Miller, *Random point processes in time and space*, 2nd ed. (Springer, 1991). [CrossRef]

**37. **M. Leutenegger, H. Blom, J. Widengren, C. Eggeling, M. Gösch, R. A. Leitgeb, and T. Lasser, “Dual-color total internal reflection fluorescence cross-correlation spectroscopy,” J. Biomed. Opt. **11**, 040502(1–3) (2006). [CrossRef] [PubMed]

**38. **M. S. Robbins and B. J. Hadwen, “The noise performance of electron multiplying charge-coupled devices,” IEEE Trans. Electron Devices **50**, 1227–1232 (2003). [CrossRef]

**39. **J. Hu, L.-s. Li, W. Yang, L. Manna, L.-w. Wang, and A. P. Alivisatos, “Linearly polarized emission from colloidal semiconductor quantum rods,” Science **292**, 2060–2063 (2001). [CrossRef] [PubMed]

**40. **B. Huang, W. Wang, M. Bates, and X. Zhuang, “Three-dimensional super-resolution imaging by stochastic optical reconstruction microscopy,” Science **319**, 810–813 (2008). [CrossRef] [PubMed]

**41. **G. Shtengel, J. A. Galbraith, C. G. Galbraith, J. Lippincott-Schwartz, J. M. Gillette, S. Manley, R. Sougrat, C. M. Waterman, P. Kanchanawong, M. W. Davidson, R. D. Fetter, and H. F. Hess, “Interferometric fluorescent super-resolution microscopy resolves 3D cellular ultrastructure,” Proc. Natl. Acad. Sci. USA **106**, 3125–3130 (2009). [CrossRef] [PubMed]

**42. **T. J. Gould, M. S. Gunewardene, M. V. Gudheti, V. V. Verkhusha, S.-R. Yin, J. A. Gosse, and S. T. Hess, “Nanoscale imaging of molecular positions and anisotropies,” Nature Methods **5**, 1027–1030 (2008). [CrossRef] [PubMed]