Asymptotics of Bayesian error probability and source super-localization in three dimensions

S. Prasad

doi:10.1364/OE.22.016008

1. Introduction

Spatial localization and tracking of single molecules at scales well below the diffraction limit are now routinely realized in flourescence-based bioimaging, as described in two excellent recent reviews [1, 2]. The unprocessed image data typically provide initial information about the coarser ranges to which the position coordinates of a single molecule in question are confined. But a more careful post-processing of the image data involving point-spread function (PSF) fitting either directly [3] or indirectly by means of either a least-squares based centroid-fitting approach [4] or statistical approaches involving maximum-likelihood estimation [5], Bayesian inference [6] or mean-squared displacement [7,8] can locate and track [9] one or more fluorescing molecules with much finer precision. It is possible to achieve a precision 10–20 times smaller -even 100–200 times smaller [10, 11] - than the standard Abbe diffraction limit, δx ≥ 0.5λ/NA, of a microscope operating at an observing wavelength λ and with numerical aperture NA.

The typical theoretical error analysis for biomolecular localization has considered the mean-squared error (MSE) in estimating the location, either for a centroid-based PSF model fit-ting algorithm [4] or in an algorithm-independent Cramér-Rao-bound (CRB) based minimum-estimator-variance analysis [12–14]. A CRB based analysis is local and incomplete, however, since it fails to account for any prior knowledge about the location based on low-resolution image data, knowledge that may be critical to the localization accuracy at low signal strengths. A local analysis also excludes any higher-order sensitivity criteria that would permit a more accurate treatment of certain problems for which a first-order sensitivity metric like the Fisher information (FI), on which the CRB is based [15], vanishes identically. An important problem of this kind is the axial localization of a molecule at zero defocus using any defocus-symmetric PSF, such as the conventional Airy-disk PSF [14]. It seems more sensible, although theoretically less tractable, to pose the problem of source localization in a Bayesian framework where any prior knowledge may be incorporated in a statistical formalism in a non-local manner.

The point-source localization problem naturally satisfies a Bayesian protocol in which the coarse spatial ranges bounding the source coordinates provide an initial, or prior, statistical distribution of these coordinates, that without additional information can be regarded as being uniform over these ranges. Under a sufficiently high signal-to-noise ratio (SNR) of fluorescence detection, the data X contain detailed information about the source location. Their conditional statistics, specified by a probability density function (PDF), P(x | m), conditioned on the knowledge of a specific source location m, determine the degree to which the location uncertainty may be improved as a function of the source flux and background and sensor noise levels.

The problem of localization with sub-diffractive errors amounts, in our multi-hypothesis-testing (MHT) based Bayesian view [16], to bounding the source location more finely by means of the posterior distribution, p(m|X), than the coarser bound provided by the uniform prior p_m, m = 1,...,M. Indeed, the mean and mode of the posterior provide two different excellent estimates of source location. The first of these estimators, the so-called minimum MSE (MMSE) estimator, upper-bounds the best achievable precision theoretically possible with any image processing protocol. The second, the maximum a posteriori (MAP) estimator, provides a different fundamental metric of minimum error, namely the minimum probability of error (MPE) in estimating the correct source position from among a set of M a priori random possible positions,

P_{e}^{(\min)} = 1 - 𝔼 [p ({\hat{m}}_{MAP} | X)],

where m̂ is the MAP estimator,

{\hat{m}}_{MAP} = \underset{m = 1, \dots, M}{argmax} p (m | X) .

In the context of a highly sensitive Bayesian detector and provided that the decision variable (e.g., source position) is parameterizable, the two metrics are closely related, as shown in our previous work [17]. The Ziv-Zakai lower bound [18] on the MMSE, which is expressed in the terms of the MPE of a related binary-hypothesis testing problem, presents a different context in which the two metrics are more generally related. For non-parameterizable hypotheses, the MMSE is not defined, but the MPE continues, very generally, to serve as a fundamental error metric. The MPE metric also provides a useful relationship between Bayesian inference and statistical information via the Fano bound [19].

The present paper extends our previous MPE based Bayesian error analysis [17] in two important ways. First, we derive useful analytical approximations for the MPE in the asymptotic regime of high SNR for the general MHT problem under rather general noise statistics. Second, we consider the more realistic asymptotic domain of many sensor data pixels in two dimensions, unlike [17] where attention was restricted to one-dimensional localization using just a few pixels along a line. It differs, however, from the more standard asymptotic analyses of MHT [20] in which one assumes that many statistically identical data frames are present, in addressing the experimentally more realistic context of many image pixels that sample a spatially varying PSF, typically possessing a single peak so the pixels progressively farther away from the peak contain progressively less signal. We present here a comprehensive analysis of the MPE within the MHT protocol, as applicable to both imaging and non-imaging problems in the presence of photon and sensor noise.

In Sec. 2, we present a brief review of our general framework for an MHT based error analysis. In Sec. 3, we illustrate our framework by applying it to the simple purely Gaussian noise model, as applicable, e.g., to additive sensor read-out noise, for which a straightforward asymptotic analysis yields an approximate form of the MPE as a sum of complementary error functions. We address in Sec. 4, under a simplifying pseudo-Gaussian approximation, generalizations of the Gaussian additive noise to include the fluctuations of the background and the signal-dependent photon noise, both described by Poisson statistics. Specifically, we combine all three noise statistics into a single Gaussian PDF with the mean being equal to the sum of the (spatially varying) signal and (spatially uniform) background mean values, and the variance being equal to the sum of the sensor noise variance and the signal and background mean values, the latter mean values also being the corresponding variances under Poisson statistics. This simplifying approach may be justified under conditions of either a large signal mean value or a large background mean value per pixel or a large additive noise variance. A continuity correction [21], not implemented here, can improve the pseudo-Gaussian approximation of the combined noise statistics still further. But even when a pseudo-Gaussian approximation cannot be justified, as for the Poisson distribution of pure shot noise at low mean signal counts, our MPE based analysis, as we show explicitly in a few instances in Sec. 5, still reflects quantitatively similar trends. In Sec. 4 we also discuss an extension to more general exponentially peaked noise PDFs, but still more general noise models that incorporate other noise sources, such as electron multiplication noise [22] of particular interest to the EMCCD detector community, are not considered. However, if their PDFs follow an exponential form, as for combined Poisson-Gaussian noise treated via the steepest-descent approximation [23], then our more general approach may again be used to derive an asymptotic expression for the MPE.

Section 5 is devoted to a practical application of our theoretical framework to discuss MPE based upper bounds on the degree of 3D localization of a single point source by different kinds of fluorescence-based imaging microscopes under varying conditions of source strength and background levels. Here we present in graphical form our numerically obtained exact and asymptotic results for the MPE for this problem. Much attention will be devoted to a detailed comparison of 2D and 3D localizations achieved by the conventional clear-aperture imaging and a recently proposed [24] pupil-phase-engineered PSF that encodes defocus via its rotation in the image plane and is thus naturally suited to perform axial localization. Our asymptotic calculations agree with the exact evaluations of the MPE, particularly in the high SNR limit for which the MPE values drop below about 0.05. Encouraged by this fact, we also performed fast asymptotic evaluations of the MPE to compare these two imagers with the astigmatic imager [25] that achieves good axial localization near the plane of best focus. Other well established rotating PSFs [26, 27] are not discussed here since their phase-mask designs, unlike that of [24], rely on complicated iterative numerical optimization and are thus not as amenable to the error analysis presented here. More exhaustive MPE based performance comparisons of different 3D imaging modalities currently in use will be studied in a future publication. Some concluding remarks appear in Sec. 6.

2. Minimum probability of error for M-ary hypothesis testing

For Bayesian discrimination among M different hypotheses, expression (1) for the MPE in inferring the correct hypothesis from data X, drawn from the set 𝒮_X, may be reformulated by means of the Bayes rule as

P_{e}^{(\min)} = 1 - \sum_{m = 1}^{M} p_{m} \int_{ℛ_{m}} d x P (x | m),

where ℛ_m is the decision region in the data set assigned to the mth hypothesis. The MAP criterion, namely

ℛ_{m} = {x | P (x | m) p_{m} > P (x | m^{'}) p_{m^{'}}, \forall m^{'} \neq m},

defines the various decision regions for which the error probability (3) takes its minimum value. The M possible hypotheses exhaust all possible decisions from any data outcome, i.e.,

𝒮_{X} = \cup_{m = 1}^{M} ℛ_{m}

. Since

\int_{𝒮_{X}} d x P (x | m) = 1

for any m and ∑_m p_m = 1, we may express Eq. (3) more conveniently as

P_{e}^{(\min)} = \sum_{m = 1}^{M} p_{m} \sum_{m^{'} \neq m} \int_{ℛ_{m^{'}}} d x P (x | m) .

We shall assume that the data space is real and multi-dimensional, as for image data typically collected at a number of pixels, say N, in the sensor plane. Thus, 𝒮_X ⊂ ℝ^N. For N >> 1, as is typically the case, and reasonably sensitive detectors, $P_{e}^{(\min)}$ may be evaluated approximately via an asymptotic analysis. For high SNR, one expects, from the law of large numbers [19], that for each m, only one decision region ℛ_m̃, the one that is “closest” to decision region ℛ_m in the sense

\tilde{m} = \underset{m^{'} \neq m}{argmax} max_{x \in ℛ_{m^{'}}} {P (x | m)},

will make a significant contribution, with all other decision regions making exponentially smaller contributions. The MPE (5) is then accurately approximated by the asymptotic expression

P_{e}^{(\min)} = \sum_{m = 1}^{M} p_{m} \int_{ℛ_{\tilde{m}}} d x P (x | m) .

We next derive an asymptotic expression for the MPE (7) for purely Gaussian data statistics.

3. Gaussian conditional data statistics

For Gaussian conditional data statistics,

P (x | m) = \frac{1}{{(2 π σ^{2})}^{N / 2}} exp [- (1 / 2) {‖ x - x_{m} ‖}_{2}^{2} / σ^{2}],

where x_m denotes the mean value of the data vector under hypothesis m and σ² the variance of data at each pixel, the definition of m̃ may equivalently be stated as

\tilde{m} = \underset{m^{'} \neq m}{argmin} min {E_{m} (x) | x \in ℛ_{m^{'}}},

where up to an additive constant E_m(x) is simply proportional to − ln[P(x|m) p_m],

E_{m} (x) ≝ {‖ x - x_{m} ‖}_{2}^{2} - 2 σ^{2} ln p_{m} .

For a given m, we determine m̃ by first mapping out the boundary between the decision region ℛ_m and other decision regions, and then finding that decision region for which E_m(x) takes the smallest possible value at the boundary. This is an optimization problem that is easily solved by requiring that at the boundary between ℛ_m and ℛ_m′, E_m(x) = E_m′ (x), i.e.,

{‖ x - x_{m} ‖}_{2}^{2} = {‖ x - x_{m} - δ x_{m m^{'}} ‖}_{2}^{2} + 2 σ^{2} ln (p_{m} / p_{m^{'}}),

where δx_mm′ simply denotes the separation vector between the mean data values for the mth and m′th hypotheses,

δ x_{m^{'} m} ≝ x_{m^{'}} - x_{m} .

Expanding the squared norm on the right hand side (RHS) of Eq. (11) using the identity

{‖ a - b ‖}_{2}^{2} = {‖ a ‖}_{2}^{2} + {‖ b ‖}_{2}^{2} - 2 a^{T} b

, where a and b are two column vectors, we see that the boundary between ℛ_m and ℛ_m′ is described by the equation

{(x - x_{m} - γ \frac{δ x_{m^{'} m}}{2})}^{T} δ x_{m^{'} m} = 0,

with γ given by

γ = 1 + \frac{2 σ^{2} ln (p_{m} / p_{m^{'}})}{{‖ δ x_{m m^{'}} ‖}_{2}^{2}} .

Equation (13) defines a hyperplane that passes through the point x_m + γδx_mm′/2 on the mean-to-mean separation vector δx_m′m and is orthogonal to that vector. Clearly, over this hyperplane the common value of E_m and E_m′ has its minimum at this point, namely at x = x_m + γδx_m′m/2, given by

{min}_{x \in ℛ_{m^{'}}} (E_{m}) = F_{m m^{'}}^{2} - 2 σ^{2} ln p_{m}

, where

F_{m m^{'}}^{2} ≝ \frac{γ^{2}}{4} {‖ x_{m} - x_{m^{'}} ‖}_{2}^{2} .

The index of the “closest” decision region to ℛ_m, as defined by Eq. (9), is then the argument of the minimum value of min_{x∈ℛ_m′} (E_m) over all m′ ≠ m, i.e.,

\tilde{m} = \underset{m^{'} \neq m}{argmin} (F_{m m^{'}}^{2} - 2 σ^{2} ln p_{m}) .

In view of the PDF (11), the asymptotically correct expression (7) may be evaluated approximately by transforming the integral on its RHS, for each m value, to a coordinate system in the N-dimensional data space for which one of the coordinate unit vectors, say t̂, is chosen to be along the separation vector x_m̃ − x_m and the remaining (N − 1) coordinate axes are chosen to span the hyperplane, orthogonal to t̂, that separates the decision region ℛ_m from its closest neighbor ℛ_m̃. The deviation of a data vector from its mean value, (x − x_m), may then be expressed in the new coordinate basis as x − x_m = x_tt̂+ x_⊥, where x_⊥ is the projection of x − x_m in the (N − 1) dimensional hyperplane orthogonal to t̂, i.e., $x_{⊥}^{T} \hat{t} = 0$ . This transformation allows us to express the squared norm in (10) as $x_{t}^{2} + {‖ x_{⊥} ‖}_{2}^{2}$ and the PDF (8) as

P (x | m) = \frac{1}{{(2 π σ^{2})}^{N / 2}} exp [- (1 / 2) (x_{t}^{2} + {‖ x_{⊥} ‖}_{2}^{2}) / σ^{2}] .

On substitutiing this form for P(x | m) in (7) and integrating the variable x_t from F_mm̃ to ∞ and the orthogonal projection x_⊥ over the full hyperplane containing it, we obtain

P_{e}^{(\min)} = \sum_{m = 1}^{M} p_{m} \frac{1}{{(2 π σ^{2})}^{1 / 2}} \int_{F_{m \tilde{m}}}^{\infty} exp [- (1 / 2) x_{t}^{2} / σ^{2}],

in which we have used the simple integral identity,

\int_{- \infty}^{\infty} d x exp [- (1 / 2) x^{2} / σ^{2}] = {(2 π σ^{2})}^{1 / 2},

(N − 1) times to integrate over the (N − 1) nonzero, mutually orthogonal components of the vector x_⊥. The integral in expression (18) is related to the complementary error function, erfc,

P_{e}^{(\min)} = \frac{1}{2} \sum_{m = 1}^{M} p_{m} erfc (F_{m \tilde{m}} / \sqrt{2 σ^{2}}),

with erfc defined as

erfc (u) ≝ \frac{2}{\sqrt{π}} \int_{u}^{\infty} exp (- x^{2}) d x .

Since $F_{m \tilde{m}} / \sqrt{2 σ^{2}}$ is proportional to the SNR for any m, it is large compared to 1 for sensitive detectors, for which expression (19) simplifies since an asymptotically valid approximation for erfc [28] may then be used,

\begin{array}{l} P_{e}^{(\min)} & \approx \frac{1}{\sqrt{2 π}} \sum_{m = 1}^{M} p_{m} \frac{σ}{F_{m \tilde{m}}} exp [- (1 / 2) F_{m \tilde{m}}^{2} / σ^{2}] \\ = \frac{1}{\sqrt{2 π}} \sum_{m = 1}^{M} p_{m} \frac{2 σ}{γ {‖ x_{\tilde{m}} - x_{m} ‖}_{2}} exp [- γ^{2} \frac{{‖ x_{\tilde{m}} - x_{m} ‖}_{2}^{2}}{8 σ^{2}}] . \end{array}

4. Pseudo-Gaussian conditional data statistics and generalization

A similar but considerably more involved treatment of the MPE may be given for a pseudo-Gaussian conditional data PDF that accurately describes the statistics of image data acquired under combined photon-number fluctuations and sensor read-out noise, whenever the mean photon count at each pixel is much larger than 1. Let the data x, given the hypothesis m, be distributed according to the PDF

P (x | m) = \frac{1}{{(2 π)}^{N / 2} {det}^{1 / 2} (Σ_{m})} exp [- (1 / 2) (x^{T} - x_{m}^{T}) Σ_{m}^{- 1} (x - x_{m})] .

For statistically independent data pixels, the data covariance matrix is a diagonal matrix,

Σ_{m} = diag (σ^{2} + x_{m}),

where σ² and x_m denote, as before, the variance of sensor read-out noise and the mean data vector, respectively, given the hypothesis m. Here diag(v) denotes a diagonal matrix formed the elements of vector v taken in order and diag(u/v) the diagonal matrix formed from the ratios of the corresponding elements of the vectors u and v. The variance of the pseudo-Gaussian PDF (22) at a pixel is the sum of the Gaussian read-out noise variance and the variance of the shot noise corresponding to Poisson photon-number fluctuations, the latter being equal to the mean photon number at that pixel.

Under asymptotic conditions, as for Gaussian data statistics, the most significant contributions to the MPE from the mth term in the sum (7) come from the vicinity of the point, x_*, on the boundary between ℛ_m and ℛ_m̃ where P(x | m) p_m has its largest value. This point does not, however, lie on the line joining the centers of the two decision regions nor is the boundary between two decision regions a hyperplane any more. Rather one must perform a maximization of P(x | m), or equivalently a minimization of −lnP(x | m),

- ln P (x | m) = (1 / 2) {(x^{T} - x_{m}^{T}) Σ_{m}^{- 1} (x - x_{m}) + ln [{(2 π)}^{N} det Σ_{m}]},

subject to the constraint that x be on the boundary, i.e., −ln[P(x | m) p_m] = − ln[P(x | m̃) p_m̃], according to the MAP decision rule underlying the MPE expression. The details of such constrained minimization by the method of Lagrange multipliers are presented in Appendix A, where we derive the following asymptotically accurate expression for the MPE:

P_{e}^{(\min)} = \frac{1}{2} \sum_{m} p_{m} erfc ({‖ U_{m} ‖}_{2} / \sqrt{2}),

with the argument of erfc given by

{‖ U_{m} ‖}_{2} = \frac{1}{2} \frac{\sum_{i = 1}^{N} \frac{{(σ^{2} + x_{m i})}^{1 / 2}}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ x_{\tilde{m} m i})}^{2}}{{[\sum_{i = 1}^{N} \frac{1}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ x_{\tilde{m} m i})}^{2}]}^{1 / 2}} .

Expression (26) is valid asymptotically under the assumptions of many pixels of data, N >> 1, and large mean combined source-background flux per pixel for which we expect ‖δx_m̃m‖₂ << ‖σ² + x_m‖₂. Numerically a somewhat more accurate value of the asymptotic expression is achieved by including two terms, rather than one, for each value of m in the sum (25), the second term corresponding to the next nearest decision region at whose boundary P(x | m) takes its next highest boundary value. We use such an improved approximation to calculate all our numerical results.

Extreme asymptotic limit

As the number of data pixels and the sensitivity of the detector grow, the MPE sum (25) will be dominated by a single term, namely that m for which the argument ${‖ U_{m} ‖}_{2} / \sqrt{2}$ of erfc is the smallest, but still large compared to 1, and the MPE asymptotically becomes a pure exponential with a characteristic exponent, ν_∞, given by

\begin{array}{l} ν_{\infty} & ≝ - lim_{N \to \infty} \frac{ln P_{e}^{(\min)}}{N} \\ = lim \frac{1}{2 N} min_{m} {‖ U_{m} ‖}_{2}^{2} \\ = \frac{1}{8} min_{m} lim_{N \to \infty} \frac{{[\frac{1}{N} \sum_{i = 1}^{N} \frac{{(σ^{2} + x_{m i})}^{1 / 2}}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ x_{\tilde{m} m i})}^{2}]}^{2}}{[\frac{1}{N} \sum_{i = 1}^{N} \frac{1}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ x_{\tilde{m} m i})}^{2}]}, \end{array}

in which all additive terms in

ln P_{e}^{(\min)}

that scale sub-linearly with N tend to 0 in the limit.

More general exponentially peaked data statistics

Our asymptotic analysis of the MPE may be extended to any data statistics that may be expressed in the exponential form,

P (x | m) = exp [- L_{m} (x; N)],

where L_m(x; N) is the negative log-likelihood function (LLF) for the mth outcome. The exponential family of PDFs [21] is an example of such distributions. We shall assume an exponential peaking for the PDF in the sense that the LLF is approximately an extensive variable in the asymptotic limit, N → ∞, i.e. lim_N→∞ L_m(x; N)/N is finite. As more and more sensor pixels are included outward from the brightest pixel in the localized image PSF, the less signal per pixel is expected to be present on average. Depending directly on the spatial extent of the PSF, this implies, in general, a threshold for the number of data pixels above which the MPE is expected to show rapidly diminishing improvement and thus essentially to saturate with N. Any requirement of the extensivity of the LLF thus need be imposed only out to such large N values.

An analysis similar to that presented for the pseudo-Gaussian noise statistics yields similar expressions for the MPE in the asymptotic limit. But the erfc in the mth term in the MPE sum (25) must be replaced a one-dimensional integral of P(x | m) from that boundary point x_* on the boundary hypersurfaces of ℛ_m with all its neighbors where P(x | m) takes its maximum value, out to ∞ along the direction of the boundary normal at that point. The logarithm of such an integral is well approximated, up to asymptotically negligible additive corrections, by the value of logP(x_* | m) at that boundary point itself. In the asymptotic limit of a large negative LLF at such a boundary point, the following expression for the MPE is thus expected to be logarithmically accurate:

P_{e}^{(\min)} = \frac{1}{2} \sum_{m} p_{m} max_{s} P (x_{* s} | m),

where x_*s, s = 1, 2,..., are all possible boundary points of ℛ_m with its neighbors, including possibly multiple points per boundary, at which P(x | m) is stationary.

We shall now illustrate our asymptotic MPE analysis for the problem of 3D super-localization in single-molecule imaging for which the combined signal-background-sensor noise may be characterized approximately by the pseudo-Gaussian statistics. Specifically, we shall consider the MPE incurred in performing this task using a rotating-PSF imager [24] and compare its performance with the conventional clear-aperture imager for which the PSF in the plane of best focus is the Airy disk. Since our preceding analysis is completely general and makes no reference to a specific algorithm, we can compare, in our Bayesian MPE-based perspective, the best achievable performance of different protocols to achieve such localization.

5. Point-source localization using conventional and rotating-PSF-based imagers

The image of a point source produced by an imager in the absence of any sensor or photon noise is proportional to its point-spread function (PSF), H(s⃗), in which s⃗ is the position vector in the 2D sensor plane. The PSF must be discretized to reflect image-plane pixelation, requiring integration of the image intensity over the area of each pixel during the recording time to generate the mean count for that pixel. In the presence of photon-number fluctuations, typically described as a Poisson random process, the count at a pixel fluctuates correspondingly. When detected by a sensor array like a CCD sensor, the read-out process adds further noise that, after calibration, may be well approximated as zero-mean Gaussian noise with a fixed variance $σ_{s}^{2}$ . The recorded count at a pixel then exhibits combined Poisson-Gaussian fluctuations from the two sources of noise, which may, under certain conditions noted in Sec. 3, be described accurately by a pseudo-Gaussian PDF (PGP) with a mean equal to the mean count at that pixel and a variance equal to the sum of the read-out noise variance and the mean pixel count. The presence of any background fluctuations, e.g., those arising from uncontrolled, largely randomly excited emitters in the imaging volume that add to the fluorescence of a specific emitter in the bioimaging problem, can be accounted for by adding a mean background count, b̄, to the mean signal count of the PGP at a pixel and including a corresponding Poisson-noise variance, b̄, in the overall variance. Assuming that the fluctuations at the different sensor pixels are statistically uncorrelated, given a specific hypothesis m, the PGP that describes this overall random process is then given by expression (22) with the following values for the mean and variance of the data pixels:

x_{m} = s_{m} + \bar{b}, Σ_{m} = diag (σ_{s}^{2} + s_{m} + \bar{b}),

where s_m is the vector obtained by remapping the 2D array of N × N data pixels into a 1D vector of N² elements, e.g., by stacking its columns consecutively one atop the next. This permits the use of matrix-algebra methods and validates the results obtained using them in the previous sections. We further assume the mean background count to be spatially uniform.

A rotating-PSF imager is a specific imaging protocol in which a superposition of pure angular-momentum states of light can be effected through a Fresnel-zone partitioning of the imaging pupil [24]. This yields a PSF that rotates uniformly, in a nearly shape- and size-invariant manner, with defocus away from the Gaussian image plane. The amount of rotation of the PSF thus encodes the axial (z) coordinate of a point source and its transverse position the (xy) coordinates of the point source. The rotating PSF, while somewhat more extended spatially in the transverse plane than the clear-aperture diffraction-limited Airy-pattern PSF, does not, unlike the latter, spread and lose sensitivity even when the defocus-dependent phase at the edge of the pupil is many waves in magnitude. This represents a trade-off between transverse and axial encoding, which we can study and compare for both kinds of imagers under varying strengths of sensor and photon noise using our Bayesian MPE based metric.

In the formal developments of this section, we keep the PSF as being general that we may choose, as needed, to be either the rotating PSF or the conventional PSF of a clear aperture without any engineered pupil phase. The defocus adds a quadratically varying radial phase in the pupil of form π(δz/λ)(ρ/l)², where λ is the mean illumination wavelength, δz the defocus distance from the in-focus object plane a distance l from the imaging pupil, and ρ is the radial distance in the circular pupil. The phase at the edge of the pupil of radius R, namely ζ ≝ π(δz/λ)(R/l)², defines a defocus phase parameter that we shall henceforth use instead of the actual defocus distance δz. Unlike the rotating PSF, the rapid spreading of the conventional PSF with increasing defocus rapidly reduces its sensitivity to encode the axial coordinate of a defocused point source, although its tighter spatial footprint in focus means greater sensitivity and resolution to encode the transverse position of a well focused point source. It is this fundamental trade-off between the decreased source-localization sensitivity and increased depth of field for the rotating PSF and their reversal for the conventional PSF that we expect to capture with our Bayesian MPE analysis. An alternative, minimum mean-squared error (MMSE) based analysis may also be given to describe this trade-off, but the two Bayesian error analyses yield similar conclusions, at least in the highly-sensitive detection limit [17].

Let IH(s⃗ − s⃗_m; z_m) be the mean intensity in the image of a point source of total brightness I located at position (s⃗_m, z_m) in the object space. For our discrete representation, we evaluated our rotating PSF on a finer grid of subpixels than the actual sensor pixel grid, shifted it by the vector amount s⃗_m in the pixel plane, and finally summed over the subpixels constituting each sensor pixel to determine the mean count recorded by the pixel. We denote such a discrete version of the shifted PSF H(s⃗ − s⃗_m; z_m) by $h_{i j}^{(m)}$ , with ij being the 2D pixel index, so the mean count recorded by the ij pixel is $K h_{i j}^{(m)}$ , where K is the source brightness I expressed in photon count units.

To determine the decision region ℛ_m̃ “closest” to the decision region ℛ_m for a given m, we required that U_m be the smallest of all (M − 1) quantities of the same form as the RHS of (26) in which m̃ is replaced by m′ and all m′ ≠ m are considered. The contribution from the next nearest decision region was also included for each value of m, as needed for a numerically improved version of the asymptotic expression (25). The mean background level, taken to be spatially uniform, was allowed to float with the mean signal strength at the brightest pixel in the conventional in-focus Airk-disk image, the ratio of the two fixed at 0.1. For such a uniform mean background level, pixels increasingly farther from the brightest pixel see progressively weaker signal levels, with the background eventually dominating the signal. This situation only gets worse when the rotating PSF with its somewhat larger footprint but far better axial sensitivity than the conventional PSF is employed for localization. For this case, even with the most compact rotating PSF, the mean background level was roughly 0.55 of the signal level at the brightest image pixel. We therefore considered only a small 12 × 12 sub-image centered at the brightest pixel.

We divide our 3D localization error analysis into two parts. In the first part, we calculate the MPE and its asymptotic value for 2D transverse super-localization by factors 2x, 4x, 8x, and 16x at a range of focal depths, corresponding to the defocus phase ranging from 0 to 16 radians in steps of 2 radians, and the reference depth resolution set at a nominal 1 radian in that phase. The second part, by contrast, computes the MPE for achieving depth localization enhancements of 2x and 4x, corresponding to 1/2 and 1/4 radian of defocus phase, respectively, at two different defocus phases, 0 and 16 radians, for the same 4 transverse super-localization factors, 2x, 4x, 8x, and 16x. These localization error probabilities were computed for the pseudo-Gaussian statistics describing the combined fluctuations of photon number, sensor read-out, and background illumination. The nominal 3D localization cell, for our purposes, is the diffraction-limited volume over which the prior is taken to be uniformly random. Its two transverse dimensions were each taken to be 4 pixels wide, while its axial dimension, as we have stated before, was taken to be 1 radian when expressed in terms of the corresponding defocus phase. These somewhat arbitrary choices can be shown to be consistent with the diffraction limited imaging criterion with respect to the parameters chosen for our PSFs. Note that much as in actual microscopy of biological samples, all our images are thus oversampled by the sensor array.

Transverse localization in 2D

The error criterion we employed for reliable localization was that the MPE be less than 5%, corresponding to a statistical confidence limit (CL) for correct match of better than 95%. Both exact and asymptotic numerical evaluations of the MPE, according to relations (3) and (25), were performed, the former via a Monte-Carlo sampling approach in which N_s data-vector samples {x⁽¹⁾,...,x^(N_s)}, with N_s sufficiently large, were drawn according to the approximate pseudo-Gaussian PDF (22), with the mean and variance given by relation (30) for each source location index m. The numerically exact MPE evaluation simply recorded the acceptance or rejection of each data sample x based on whether P(x | m) was larger than all other P(x | m′), m′ ≠ m, or not. Subtracting from 1 the acceptance rate averaged over all M values of the hypothesis label m then yielded the exact MPE acoording to Eq. (3). The results we present here were obtained for N_s = 5000, but our evaluations changed by no more than a percent even when we increased N_s to 20000 and still larger values in certain cases, giving us excellent confidence in the accuracy of our numerical results.

We present in Fig. 1 both the exact and asymptotic values of the MPE for a number of source strengths and varying degrees of 2D sub-diffractive transverse localization, but without any axial, or depth, localization, for the conventional imager at two different depths, namely ζ = 0 and ζ = 16 radians. The sensor read-out noise variance was held fixed at $σ_{s}^{2} = 1$ . Since the conventional imager loses sensitivity quickly with increasing depth, the two depth values are attended by a large disparity in the corresponding values of the MPE. By contrast, the rotating-PSF imager has rather similar values of the MPE over a considerable depth range, necessitating two different figures, Figs. 2(a) and 2(b), to display them properly for the same two values of the depth. For the in-focus case (ζ = 0), the conventional PSF based imager yields significantly smaller error probabilities than the rotating PSF based imager for the transverse localization task, regardless of the value of the enhancement factor M_⊥ that we considered here, namely 2, 4, 8, and 16. This is consistent with the fact that the conventional PSF is more compact than the rotating PSF when in focus. However, the behavior reverses dramatically at the large defocus phase, ζ = 16 radians, since the rotating PSF maintains its shape and compactness over such large changes of defocus, and thus the MPE values for 2D localization are quite comparable both in focus and at large defocus.

Fig. 1 (a) Doubly logarithmic plots of MPE vs. source brightness K for the conventional imager for two different values of the defocus phase, ζ, namely 0 and 16 radians. The plots for the various transverse super-localization factors are indicated explicitly both by means of their different colors, marker symbols, and labels, with the curves in the bottom half referring to ζ = 0 (in-focus) and those in the top half to ζ = 16. The corresponding asymptotic values of the MPE are connected by dashed line segments. The values of the MPE for 2x super-localization in the in-focus case are both too small and at too small K to appear in the figure, while the asymptotic results for the 8x and 16x super-localizations for ζ = 16 are suppressed for clarity. (b) The same plots are now shown on a linear vertical scale for ease of comparison with Fig. 2.

Download Full Size | PDF

Fig. 2 Plots of MPE vs. source brightness K, in dB, for the rotating-PSF imager for the same two values of the defocus phase, ζ, (a) 0 and (b) 16 radians, as in Fig. 1. The plots of the MPE for the various transverse super-localization factors are indicated explicitly both by means of their different colors, marker symbols, and labels. The corresponding asymptotic values are connected by dashed line segments.

Download Full Size | PDF

Based on the comparisons presented in Figs. 1 and 2, it is clear that with increasing defocus the rotating PSF based imager outstrips the 2D localization performance of the conventional imager at some fairly low value of the defocus phase, ζ, but in a manner that depends on the specific values of K and M_⊥. In Figs. 3(a)–(c), we present the comparison of the MPE based performance of the two imagers as a function of the defocus phase for a low and a high value of K. The reversal of the behavior of the error probabilities for 2D super-localization with increasing defocus is easily seen in these plots. The cross-over of the error curves for the two kinds of imagers occurs at fairly low defocus since the conventional PSF blurs rapidly with defocus and thus fails to provide good transverse super-localization at all but the lowest values of the defocus. For both imagers, the error probabilities rise, as expected, with the degree of 2D super-localization sought, requiring increasingly larger K to keep the MPE acceptably low, but the rotating-PSF based imager continues to provide excellent localization enhancement at fairly modest values of K.

Fig. 3 Plots of MPE vs. the defocus phase, ζ, for the rotating-PSF imager (blue curves) and the conventional imager (green curves) for two different values of K. The three subfigures refer to the values (a) 2x; (b) 4x; and (c) 16x of the 2D super-localization factor. The dashed lines display the corresponding asymptotic results except for the rotating PSF imager at K = 10³.

Download Full Size | PDF

We next tested the accuracy of the pseudo-Gaussian approximation to simulate the pure Poisson statistics of signal count and background in the absence of any read-out noise, σ_s = 0. For each source location, labeled by m and characterized by the mean signal-background sum count vector x_m, we drew N_s = 5000 image frames according to Poisson statistics using the poissrnd program in Matlab. By evaluating the rate of acceptance of such Monte-Carlo draws according to the MAP rule, as before, and averaging over all possible source locations, we obtained a numerically exact value of the MPE (3). We plot in Figs. 4(a)–(c) these values versus the source brightness K along side the corresponding pseudo-Gaussian-statistics based results for the two imagers for three different values of M_⊥, namely 2, 4, and 8, for both in-focus and defocused (ζ = 16) 2D super-localization. An examination of these figures reveals that the pseudo-Gaussian approximation for purely Poisson statistics fails at low photon numbers where the MPE curves for the two probability distributions deviate at the level of 15–20% from each other. However, at intermediate to large signal counts per pixel, as for K ∼ 5 × 10² − 5 × 10³, the approximate and exact distributions produce MPE results that are within a few percent of each other. The apparent failure of the approximation at very large counts, as seen for the 8x super-localization curves around K ∼ 10⁴ for the rotating PSF imager [Figs. 4(a)–(b)], is most likely due to insufficient numerical accuracy in the Monte-Carlo evaluations of the exact expression (3) for one or both PDFs in the limit of small MPE for which the small deviations of the acceptance rate of the Monte-Carlo draws from 1 are presumably being computed inaccurately. For the conventional imager, the approximate and exact probability distributions produce virtually indistinguishable MPE results, particularly for the highly defocused case (dashed and dot-dashed lines), as seen from the doubly logarithmic plots of Fig. 4(c).

Fig. 4 Plots of MPE vs. source brightness K, in dB, for the rotating-PSF imager for (a) ζ = 0 (in-focus); (b) ζ = 16 (defocused); and for the conventional imager (c) for ζ = 0 and 16. Three different values of M_⊥, namely 2, 4, and 8, were considered here. The solid line segments refer to the approximate pseudo-Gaussian distribution, while the dashed line segments refer to the Poisson distribution, which is exact for the present case of zero sensor noise.

Download Full Size | PDF

From our graphs of the MPE in Figs. 1 and 2, we note that the asymptotic values of the MPE are quite accurate whenever they fall below about 0.05. Since this also constitutes a well accepted statistical confidence threshold (SCT) of a better than 95% error-free performance rate, we regard our asymptotic expressions as being quite useful for the superlocalization problem treated here. We demonstrate their value in deriving important scaling laws for a different but related application of optical super-resolution of a pair of closely spaced point sources in a separate paper [29].

From the results plotted in Figs. 1 and 2 and those for other defocus values not and similar results for other values of the defocus phase ζ, not shown here, we can also read off the minimum requirements on K to achieve, at the 95% statistical confidence limit, or 5% MPE, a specific value of the 2D super-localization factor M_⊥. We plot in Fig. 5 the minimum source photon number, K_min, as a function of $M_{⊥}^{2}$ for both the conventional and rotating-PSF imager for four different values of the defocus phase, namely 0, 4, 8, and 16 radians. As expected, for the rotating-PSF imager, K_min, for each value of M_⊥, increases rather modestly as it is defocused more and more over the full 16-radian range, while for the conventional imager defocused localization even for ζ = 4 requires roughly double the K_min needed for the former imager operating at ζ = 16. Only at best focus, ζ = 0, does the conventional imager deliver a better localization performance at more modest photon numbers. Its PSF spreading with increasing defocus is simply too dramatic for it to stay competitive with the rotating-PSF imager. Although our results were obtained by fixing the ratio of the mean background to the peak brightness of the conventional in-focus image at 10% and negligible sensor noise, we have verified this relative performance of the two imagers more generally. The approximately linear form of the plots in Fig. 5 for both imagers confirms the approximately quadratic dependence of the minimum source strength on the degree of super-localization sought at any value of the defocus. This conclusion is quite consistent with previous MSE based analyses [4, 12] in the signal dominated regime of operation.

Fig. 5 Plots of K_min vs. $M_{⊥}^{2}$ for the two imagers for four different values of the defocus phase, ζ, namely 0, 4, 8, and 16 radians. For the rotating-PSF imager, results shown for ζ = 0 and 16 radians bracket those for the intermediate values of ζ, not shown here. The dashed lines away from marker symbols for the conventional imager are extrapolations of our numerical results.

Download Full Size | PDF

Full 3D localization

We now address the problem of localizing the spatial position of a point source in all three dimensions by means of an imager that employs either the rotating PSF or the conventional PSF. As is well known [14], the axial localization provided by the conventional imager is rather poor for a source at the plane of best focus because of a lack of first-order sensitivity of the PSF relative to the defocus at this plane. The rotating PSF imager, however, has no such insensitivity, and is helped further by its ability to maintain a tight PSF profile over a large defocus range.

In Fig. 6, we display, for the rotating PSF imager, the MPE as a function of the source photon number for the same 10% background level used in the previous figures, but now for different degrees of transverse and axial localizations, M_⊥ and M_‖, at two different depths, ζ = 0 [Fig. 6(a)] and ζ = 16 [Fig. 6(b)]. As expected, the MPE increases with increasing demand on the degree of depth localization from 2x to 4x for each of the transverse localization enhancement factors, 2x, 4x, 8x, and 16x, at each reference depth. The rather modest differences between the overall behaviors presented in the two figures are a result of the approximate invariance of the rotating PSF shape and size across a large range of defocus. By contrast, the next two figures, Figs. 7(a) and 7(b), which present the corresponding results for the conventional imager, indicate a much larger MPE, even at zero defocus because of its first-order insensitivity in providing any depth discrimination at the in-focus plane. However, since the Bayesian error probability is a global sensitivity metric, the higher order sensitivity of conventional imagers to yield depth discrimination is fully accounted for in our results.

Fig. 6 Plots of MPE vs. K, in dB, for two different values of defocus, (a) ζ = 0 and (b) ζ = 16, for two different axial and four different transverse localization enhancement factors for the rotating-PSF imager. The latter factors are indicated by different marker symbols, namely circle for 2x, + for 4x, square for 8x, and diamond for 16x, while the former are indicated by the line type, solid for 2x and dashed for 4x.

Download Full Size | PDF

Fig. 7 Same as in Figs. 6(a) and 6(b), except here for the conventional imager.

Download Full Size | PDF

The competition between axial and transverse localizations is also evident in these figures. With increasing signal strength, the behavior of the MPE is strongly influenced by any changes in the requirement of transverse localization, rising as it does with increasing values of the latter, but ultimately at sufficiently high strengths the MPE is limited by the requirement of axial localization alone, as seen in the asymptotic behavior of the various MPE curves for fixed M_‖ but for different values of M_⊥. This behavior is quite universal for both imagers at sufficiently large signal strengths for each axial localization, but it is particularly noticeable for the in-focus (ζ = 0) conventional imager for which the 2x and 4x transverse super-localization curves are essentially indistinguishable, independent of the degree of axial super-localization demanded, over the full range of signal strengths considered here.

3D localization by the astigmatic imager

We note from Figs. 1 and 3 that the asymptotic approximation based on (25) is quite accurate whenever the MPE values become small, specifically as they fall below the CL threshold of 0.05. Robustness of performance requires such low MPE values, which means that the asymptotic approximation can be a useful tool in the practical design of imagers purposed for single-molecule localization and tracking.

We now use an asymptotic MPE evaluation to compare the full 3D localization performance of our rotating-PSF imager with another popular 3D imager, the astigmatic imager [25] that employs a cylindrical lens to achieve an anisotropic PSF. We consider here a simplified version of that imager whose aperture is square shaped with sides that are parallel to the principal sections of the cylindrical lens. The degree of lens astigmatism is described by a single parameter η = πδf R²/(2λf²), where δf is the difference between the two focal lengths, f ± δf/2, of the cylindrical lens and 2R is each side of the square lens aperture. Under the assumption δf << f, the focal astigmatism adds an optical phase of form η(x² − y²) in the pupil, where Rx, Ry are the Cartesian coordinates in the square pupil.

In Figs. 8(a) and (b), we plot the asymptotic values of the MPE versus the source photon number for an astigmatic imager with η =2 for 3D super-localization by factors M_‖ = 2, 4 and M_⊥ = 2, 4, 8, 16. Note that due to a nonzero η parameter, the 3D localization minimum error probabilities are greatly reduced at the mean focal position, which we call in-focus localization [Fig. 8(a)]. But the localization away from the mean focal position, at focal depth ζ = 8 [Fig. 8(b)], is not particularly better than that for the clear-aperture spherical lens imager, as seen from a comparison with the MPE plots of Fig. 7(b). The asymptotic values of the MPE for the rotating-PSF imager, shown in Figs. 8(c) and (d), indicate, on the other hand, that while it does not perform as well as the astigmatic imager in-focus [Fig. 8(c)], its MPE based error performance is strikingly more uniform with respect to defocus. We regard the rotating-PSF-based imager as achieving an optimal 3D super-localization performance over a large axial depth range around best focus.

Fig. 8 Plots of asymptotic values of the MPE vs. source strength for the astigmatic imager for (a) in-focus and (b) 8 rad out of focus 3D localization. The plots are labeled by the transverse super-localization factors in the same way as in Figs. 6 and 7. The corresponding rotating-PSF-imager results are shown in (c) and (d).

Download Full Size | PDF

6. Concluding remarks

We have considered in this paper an asymptotic analysis of the Bayesian problem of multi-hypothesis testing and then applied it to treat the fidelity of point-source localization with sub-diffractive error. We applied our exact and approximate analyses of the minimum probability of error (MPE) in discriminating among $M_{⊥}^{2} \times M_{‖}$ possible outcomes of the source position inside an elementary base resolution volume subdivided uniformly into $M_{⊥}^{2} \times M_{‖}$ subvolumes. The image data were drawn from a small sub-image centered around the brightest pixel in the full image. Three different imaging systems, two based on the clear-aperture spherical-lens and cylindrical-lens PSFs and a third on a phase engineered rotating PSF, were compared in detail for their MPE-based potential to achieve 3D source super-localization.

In the signal-dominated regime for which we have presented our detailed calculations here, we confirmed a number of conclusions drawn by previous researchers using mean-squared error (MSE) based analyses about the minimum source signal strength needed for super-localizing a point source spatially within the diffraction limited volume of order (λ/NA)² × λ/NA². We confirmed, in particular, a quadratic ( $M_{⊥}^{2}$ ) dependence of the minimum source strength needed to achieve a transverse localization improvement factor of M_⊥ at a statistical confidence threshold of 95% or better. The agreement in the predictions of MSE and MPE based analyses in the high-sensitivity asymptotic limit is a consequence of the equivalence of the two metrics in that limit [17].

For combined 3D transverse-axial localization, we demonstrated an interplay between the axial and transverse localization enhancements. We found that it is the axial localization enhancement that typically limits 3D localization at sufficiently high signal strengths, at which all of the transverse localization improvements from 2x to 16x entail no added error cost. The reduced sensitivity of the conventional imager to perform any axial resolution at zero defocus, noted previously in a local error-bound analysis [14], is also borne out well in our more global MPE based Bayesian analysis.

We finally applied our asymptotic MPE analysis to evaluate the full 3D localization performance of the cylindrical-lens-based astigmatic imager, and demonstrated its superiority over both the conventional and rotating PSF imagers when in focus, but, much like the conventional spherical-lens imager, a rapidly diminishing performance with increasing axial depth. By contrast, the rotating-PSF imager continues to perform uniformly well over large axial depths.

It is worth noting that the MPE is the best error performance one can achieve from the view point of error probability in a Bayesian detection protocol. Most actual localization algorithms are likely to be sub-optimal. Reasons for this sub-optimality are many, not the least of which are an incomplete identification of the statistical sources of error and their imperfect statistical characterization. Effort should be devoted primarily in mitigating these systematic sources of error before employing the MAP estimate for localization. Without such mitigation, the performance bound presented by our MPE analysis may seem overly optimistic under more realistic operating conditions.

Appendix A. Derivation of MPE asymptotics under pseudo-Gaussian noise fluctuations

In view of the form (24) for the negative log-likelihood function (LLF), finding the boundary point x_* on the boundary between ℛ_m and ℛ_m̃, amounts, via the use of a Lagrange multiplier λ, to the minimization,

min_{x} (x^{T} - x_{m}^{T}) Σ_{m}^{- 1} (x - x_{m}) - λ [(x^{T} - x_{m}^{T}) Σ_{m}^{- 1} (x - x_{m}) - (x^{T} - x_{\tilde{m}}^{T}) Σ_{\tilde{m}}^{- 1} (x - x_{\tilde{m}})],

which is performed by taking its gradient w.r.t. x and setting it to zero at x = x_*,

(1 - λ) Σ_{m}^{- 1} (x_{*} - x_{m}) + λ Σ_{\tilde{m}}^{- 1} (x_{*} - x_{\tilde{m}}) = 0 .

This equation is readily solved for (x_* − x_m) by writing (x_* − x_m̃) = (x_* − x_m) − δx_m̃m, where δx_m̃m is defined by relation (12). The solution is simplified by employing the definition (23) for the diagonal covariance matrices Σ_m and Σ_m̃ and then performing simple algebra,

(x_{*} - x_{m}) = diag {[1 + \frac{(1 - λ)}{λ} \frac{(σ^{2} + x_{\tilde{m}})}{(σ^{2} + x_{m})}]}^{- 1} δ x_{\tilde{m} m} .

A similar expression for (x_* − x_m̃) also follows from (33),

\begin{array}{l} (x_{*} - x_{\tilde{m}}) & = (x_{*} - x_{m}) - δ x_{\tilde{m} m} \\ = - diag {[1 + \frac{λ}{(1 - λ)} \frac{(σ^{2} + x_{m})}{σ^{2} + x_{\tilde{m}}}]}^{- 1} δ x_{\tilde{m} m} . \end{array}

These expressions, when substituted into the constraint, −ln[P(x | m) p_m] = − ln[P(x | m̃) p_m̃], with the negative LLF given by expression (24), yield the following equation for λ:

δ x_{\tilde{m} m}^{T} diag {\frac{λ^{2} (σ^{2} + x_{m}) - {(1 - λ)}^{2} (σ^{2} + x_{\tilde{m}})}{{[λ (σ^{2} + x_{m}) + (1 - λ) (σ^{2} + x_{\tilde{m}})]}^{2}}} δ x_{\tilde{m} m} = ln \frac{p_{m}^{2} det Σ_{\tilde{m}}}{p_{\tilde{m}}^{2} det Σ_{m}} .

Under asymptotic conditions of many pixels, N >> 1, and large combined source-background flux per pixel, we expect the inequality, ‖δx_m̃m‖₂ << ‖σ² + x_m‖₂, to hold, since the PSFs for the two most closely spaced source-location alternatives give rise to pixel counts that are correspondingly close, pixel by pixel. Since (σ² + x_m̃) = (σ² + x_m) + δx_m̃m, we may then expand the element-wise quotient, r_m, of the two vectors inside the braces in Eq. (35) in powers of δx_m̃m/(σ² + x_m) as

\begin{array}{l} r_{m} & = \frac{(2 λ - 1) (σ^{2} + x_{m}) - {(1 - λ)}^{2} δ x_{\tilde{m} m}}{{[λ (σ^{2} + x_{m}) + (1 - λ) (σ^{2} + x_{\tilde{m}})]}^{2}} \\ = \frac{(2 λ - 1)}{(σ^{2} + x_{m})} [1 - \frac{(3 λ - 1) (1 - λ)}{(2 λ - 1)} \frac{δ x_{\tilde{m} m}}{(σ^{2} + x_{m})} + O (δ^{2})], \end{array}

where δ is a shorthand notation for the order of the asymptotically small elements of the element-wise quotient vector δx_m̃m/(σ² + x_m). The RHS of Eq. (35) may be similarly expanded by writing the logarithm of the ratio of the determinants of the diagonal covariance matrices, Σ_m̃ and Σ_m, as a sum of the logarithms of the ratios of their diagonal elements, ratios that, in the asymptotic limit, differ from 1 by the small ratios of the corresponding elements of δx_m̃m and (σ² + x_m). On performing this expansion to the lowest order in these ratios and substituting expression (36) in Eq. (35), we obtain a simple equation for λ in this order, from which λ may be determined as

λ = \frac{1}{2} + \frac{1}{2} \frac{\sum_{i = 1}^{N} \frac{δ x_{\tilde{m} m i}}{(σ^{2} + x_{m i})} + 2 ln \frac{p_{m}}{p_{\tilde{m}}}}{\sum_{i = 1}^{N} \frac{δ x_{\tilde{m} m i}^{2}}{(σ^{2} + x_{m i})}} .

The mth term in the sum (7) for the MPE may be calculated in the asymptotic regime by recognizing that the integral of P(x | m) over ℛ_m̃ is dominated by its maximum value, P(x_* | m), on the boundary between ℛ_m and ℛ_m̃. More precisely, we first determine the unit normal n_* to the boundary at the point x_*, and then resolve the vector,

v ≝ Σ_{m}^{- 1 / 2} (x - x_{m}),

whose negative squared norm occurs in the exponent of the Gaussian form of P(x | m), along n_* and normal to it. The integral of P(x | m) over ℛ_m̃ may be approximately evaluated by fixing the origin of an orthonormal coordinate system at x_*, with coordinate axes along n_* and orthogonal to it.

The unit normal n_* on the boundary between regions ℛ_m and ℛ_m̃ is along the gradient of the negative LLF difference, −[lnP(x | m) − lnP(x | m̃)], evaluated at the boundary point x_*,

\begin{array}{l} n_{*} & = - \frac{\nabla_{*} ln P (x_{*} | m) - \nabla_{*} ln P (x_{*} | \tilde{m})}{{‖ \nabla_{*} ln P (x_{*} | m) - \nabla_{*} ln P (x_{*} | m) ‖}_{2}} \\ = \frac{Σ_{m}^{- 1} (x_{*} - x_{m}) - Σ_{\tilde{m}}^{- 1} (x_{*} - x_{\tilde{m}})}{{‖ Σ_{m}^{- 1} (x_{*} - x_{m}) - Σ_{\tilde{m}}^{- 1} (x_{*} - x_{\tilde{m}}) ‖}_{2}}, \end{array}

where Eq. (24) was used to reach the second inequality. From Eq. (32) and the fact that λ ≈ 1/2 from expression (37) under conditions used to derive that expression, we see that expression (39) for n_* simplifies approximately to the form

\begin{array}{l} n_{*} & = \frac{Σ_{m}^{- 1} (x_{*} - x_{m})}{{‖ Σ_{m}^{- 1} (x_{*} - x_{m}) ‖}_{2}} \\ = \frac{Σ_{m}^{- 1} (x_{*} - x_{m})}{{[(x_{*}^{T} - x_{m}^{T}) Σ_{m}^{- 2} (x_{*} - x_{m})]}^{1 / 2}}, \end{array}

where we used the definition of the norm,

{‖ v ‖}_{2}^{2} = v^{T} v

, to arrive at the second equality. We now write v = v_‖ + v_⊥, where

\begin{array}{l} v_{‖} & = n_{*} n_{*}^{T} v \\ = n_{*} n_{*}^{T} Σ_{m}^{- 1 / 2} (x - x_{*}) + n_{*} n_{*}^{T} Σ_{m}^{- 1 / 2} (x_{*} - x_{m}) \\ = u_{m} + U_{m} and \\ v_{⊥} & = v - v_{‖} \end{array}

are the projections of v along n_* and in its orthogonal complement, respectively, with the former subdivided further into two parts, u_m and U_m,

u_{m} = n_{*} n_{*}^{T} Σ_{m}^{- 1 / 2} (x - x_{*}); U_{m} = n_{*} n_{*}^{T} Σ_{m}^{- 1 / 2} (x_{*} - x_{m}),

in which U_m defines the shift vector between the mean data point x_m within ℛ_m to the origin, at x_*, of the new coordinate system whose axes are individually scaled by the diagonal elements of the diagonal matrix

Σ_{m}^{- 1 / 2}

.

In view of the definition (38) and the decomposition v = v_‖ + v_⊥, as given in Eq. (41), expression (24) may now be exponentiated to arrive at a simplified form for P(x | m),

P (x | m) = \frac{1}{{[{(2 π)}^{N} det Σ_{m}]}^{1 / 2}} exp [- (1 / 2) (v_{‖}^{T} v_{‖} + v_{⊥}^{T} v_{⊥})] .

The integral of P(x | m) over x in the decision region ℛ_m̃ can now be performed approximately as the product of the integral over the variable ‖v_‖‖₂ from ‖U_m‖₂ to ∞ and (N − 1) integrals over the remaining (N − 1) orthogonal components of v, each of the latter integrals having its limits extended to ±∞. The scaling of the coordinate axes in going from the x space to the v space, according to definition (38), exactly cancels out the determinantal factor in the denominator of expression (43), while the (N − 1) infinite integrals over the orthogonal complement of n_* produce merely an overall factor (2π)^(N−1)/2, leaving just a single Gaussian integral over v_‖ to be done. In other words, in the asymptotic limit the following approximate value may be obtained for the overall multiple integral:

\begin{matrix} \int_{ℛ_{\tilde{m}}} P (x | m) d x = \frac{1}{{(2 π)}^{1 / 2}} \int_{{‖ U_{m} ‖}_{2}}^{\infty} exp [- (1 / 2) {‖ v_{‖} ‖}_{2}^{2}] d {‖ v_{‖} ‖}_{2} \\ = \frac{1}{2} erfc ({‖ U_{m} ‖}_{2} / \sqrt{2}), \end{matrix}

which yields the following asymptotically valid expression for the MPE (3):

P_{e}^{(\min)} = \frac{1}{2} \sum_{m} p_{m} erfc ({‖ U_{m} ‖}_{2} / \sqrt{2}) .

According to Eq. (42), ‖U_m‖₂ is simply the number $| n_{*}^{T} Σ_{m}^{- 1 / 2} (x_{*} - x_{m}) |$ . Using expression (40) for n_*, we may write it as

{‖ U_{m} ‖}_{2} = \frac{(x_{*}^{T} - x_{m}^{T}) Σ_{m}^{- 3 / 2} (x_{*} - x_{m})}{{[(x_{*}^{T} - x_{m}^{T}) Σ_{m}^{- 2} (x_{*} - x_{m})]}^{1 / 2}}

In view of the relation (33) between (x_* − x_m) and δx_m̃m and since λ ≈ 1/2 according to result (37), we may express Eq. (46) approximately as

{‖ U_{m} ‖}_{2} = \frac{1}{2} \frac{δ x_{\tilde{m} m}^{T} diag [\frac{{(σ^{2} + x_{m})}^{2}}{{(σ^{2} + {\bar{x}}_{m})}^{2}}] Σ_{m}^{- 3 / 2} δ x_{\tilde{m} m}}{{δ x_{\tilde{m} m}^{T} diag [\frac{{(σ^{2} + x_{m})}^{2}}{{(σ^{2} + {\bar{x}}_{m})}^{2}}] Σ_{m}^{- 2} δ x_{\tilde{m} m}}^{1 / 2}},

where x̄_m ≝ (1/2)(x_m + x_m̃) is the arithmetic mean of the mean data vectors. Using the definition (23) for Σ_m, we may finally express ‖U_m‖₂ as the ratio

{‖ U_{m} ‖}_{2} = \frac{1}{2} \frac{\sum_{i = 1}^{N} \frac{{(σ^{2} + x_{m i})}^{1 / 2}}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ {\tilde{x}}_{\tilde{m} m i})}^{2}}{{[\sum_{i = 1}^{N} \frac{1}{{(σ^{2} + {\bar{x}}_{m i})}^{2}} {(δ x_{\tilde{m} m i})}^{2}]}^{1 / 2}} .

Acknowledgments

Helpful conversations with R. Kumar, S. Narravula, J. Antolin, Z. Yu, H. Pang, K. Lidke, R. Heintzmann, and R. Ober are gratefully acknowledged. The work reported here was supported in part by AFOSR under grant numbers FA9550-09-1-0495 and FA9550-11-1-0194.

References and links

1. H. Deschout, F. Zanacchi, M. Mlodzianoski, A. Diaspro, J. Bewersdorf, S. Hess, and K. Braeckmans, “Precisely and accurately localizing single emitters in fluorescence microscopy,” Nat. Methods 11, 253–266 (2014). [CrossRef] [PubMed]

2. A. Small and S. Stahlheber, “Fluorophore localization algorithms for super-resolution microscopy,” Nat. Methods 11, 267–279 (2014). [CrossRef] [PubMed]

3. C. Anderson, G. Georgiou, I. Morrison, G. Stevenson, and R. Cherry, “Tracking of cell surface receptors by fluorescence digital imaging microscopy using a charge-coupled device camera,” J Cell Sci. 101, 415–425 (1992).

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

5. N. Bobroff, “Position measurement with a resolution and noise-limited instrument,” Rev. Sci. Instrum. 57, 1152–1157 (1986). [CrossRef]

6. J. Yoon, A. Bruckbauer, W. Fitzgerald, and D. Klenerman, “Bayesian inference for improved single molecule fluorescence tracking,” Biophys. J. 94, 4932–4947 (2008). [CrossRef] [PubMed]

7. N. Monnier, S.-M. Guo, M. Mori, J. He, P. Lenart, and M. Bathe, “Bayesian approach to MSD-based analysis of particle motion in live cells,” Biophys. J. 103, 616–626 (2012). [CrossRef] [PubMed]

8. S. Cox, E. Rosten, J. Monypenny, T. Jovanovic-Talisman, D. Burnette, J. Lippincott-Schwartz, G. Jones, and R. Heintzmann, “Bayesian localization microscopy reveals nanoscale podosome dynamics,” Nat. Methods 9, 195–200 (2012). [CrossRef]

9. A. Gahlmann and W. Moerner, “Exploring bacterial cell biology with single-molecule tracking and super-resolution imaging,” Nat. Rev. Microbiol. 12, 9–22 (2014). [CrossRef]

10. A. Yildiz, J. Forkey, A. McKinney, T. Ha, Y. Goldman, and P. Selvin, “Myosin V walks hand-over-hand: Single fluorophore imaging with 1.5-nm localization,” Science 300, 2061–2065 (2003). [CrossRef] [PubMed]

11. A. Pertsinidis, Y. Zhang, and S. Chu, “Subnanometre single-molecule localization, registration and distance measurements,” Nature 466, 647–651 (2010). [CrossRef] [PubMed]

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

13. S. Ram, E. Sally Ward, and R. Ober, “Beyond Rayleigh’s criterion: a resolution measure with application to single-molecule microscopy,” Proc. Natl. Acad. Sci. U. S. A. 103, 4457–4462 (2006). [CrossRef] [PubMed]

14. S. Ram, E. Sally Ward, and R. Ober, “How accurately can a single molecule be localized in three dimensions using a fluorescence microscope?” Proc. SPIE 5699, 426–435 (2005). [CrossRef] [PubMed]

15. S. Kay, Fundamentals of Statistical Signal Processing: I. Estimation Theory (Prentice Hall, 1993), Chap. 3.

16. H. van Trees and K. Bell, Detection, Estimation, and Modulation Theory (Wiley, 2013), Part I.

17. S. Prasad, “New error bounds for M-testing and estimation of source location with subdiffractive error,” J. Opt. Soc. Am. A 29, 354–366 (2012). [CrossRef]

18. J. Ziv and M. Zakai, “Some lower bounds on signal parameter estimation,” IEEE Trans. Inform. Theory 15, 386–391 (1969). [CrossRef]

19. T. Cover and J. Thomas, Elements of Information Theory (Wiley, 1991). [CrossRef]

20. C. Leang and D. Johnson, “On the asymptotics of M-hypothesis Bayesian detection,” IEEE Trans. Inform. Theory 43, 280–282 (1997). [CrossRef]

21. G. Roussas, A Course in Mathematical Statistics (Academic, 1997).

22. J. Chao, E. Sally Ward, and R. Ober, “Fisher information matrix for branching processes with application to EMCCDs,” Multidimens. Syst. Signal Process. 23, 349–379 (2012). [CrossRef] [PubMed]

23. D. Snyder, C. Helstrom, A. Lanterman, M. Faisal, and R. White, “Compensation for readout noise in CCD images,” J. Opt. Soc. Am. A 12, 272–287 (1995). [CrossRef]

24. S. Prasad, “Rotating point spread function via pupil-phase engineering,” Opt. Lett. 38, 585–587 (2013). [CrossRef] [PubMed]

25. 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]

26. S. Pavani, M. Thompson, J. Biteen, S. Lord, N. Liu, R. Twieg, R. Piestun, and W. Moerner, “Three-dimensional, single-molecule fluorescence imaging beyond the diffraction limit by using a double-helix PSF,” Proc. Natl. Acad. Sci. U. S. A. 106, 2995–2999 (2009). [CrossRef] [PubMed]

27. M. Lee, S. Lew, M. Badieirostami, and W. Moerner, “Corkscrew point spread function for far-field three-dimensional nanoscale localization of pointlike objects,” Opt. Lett. 36, 202–204 (2011). [CrossRef] [PubMed]

28. M. Abramowitz and I. Stegun, eds., Handbook of Mathematical Functions (Dover, 1964), Chap. 7.

29. S. Prasad, “Asymptotics of Bayesian error prability and 2D pair superresolution,” Opt. Express 22, 16029–16047 (2014).

Asymptotics of Bayesian error probability and source super-localization in three dimensions

Abstract

1. Introduction

2. Minimum probability of error for M-ary hypothesis testing

3. Gaussian conditional data statistics

4. Pseudo-Gaussian conditional data statistics and generalization

Extreme asymptotic limit

More general exponentially peaked data statistics

5. Point-source localization using conventional and rotating-PSF-based imagers

Transverse localization in 2D

Full 3D localization

3D localization by the astigmatic imager

6. Concluding remarks

Appendix A. Derivation of MPE asymptotics under pseudo-Gaussian noise fluctuations

Acknowledgments

References and links

Cited By

Figures (8)

Equations (49)

Optics Express