Fitting point-spread function models

Least-squares methods

To approximate the data with a point-spread function, least-squares methods [4, 1, 5] are employed to minimize the sum of (weighted) squared residuals defined by

$\chi^{2}\left(\boldsymbol{\theta}\mid\mathcal{D}\right)=\sum_{x,y\in\mathcal{D% }}{w\left(\tilde{I}\left(x,y\right)-\mathrm{PSF}\left(x,y\mid\boldsymbol{% \theta}\right)\right)^{2}}\,.$

(1)

Here the residual value for the $\left(x,y\right)$ data point is defined as the difference between the observed image intensity $\tilde{I}\left(x,y\right)$ and the value approximated by the $\mathrm{PSF}\left(x,y\mid\boldsymbol{\theta}\right)$ , where $\boldsymbol{\theta}$ are the PSF parameters. The residual value can be further weighted by , making all measurements equally significant, or weighted by $w=1/\tilde{I}\left(x,y\right)$ , which takes into account the uncertainty in the number of detected photons.

The search for parameters $\boldsymbol{\hat{\theta}}$ which minimize $\chi^{2}\left(\boldsymbol{\theta}\mid\mathcal{D}\right)$ , leads to an optimization problem formulated as

$\boldsymbol{\hat{\theta}}=\argmin_{\boldsymbol{\theta}}\chi^{2}\left(% \boldsymbol{\theta}\mid\mathcal{D}\right)\,,$

(2)

which we solve by the Levenberg-Marquardt algorithm as implemented in the Apache Commons Math library [2]. The starting point for the optimization process is computed from the data as the difference between the maximum and the minimum intensity values for the molecular intensity $\theta_{N}$ , and as the minimum intensity value for the background offset $\theta_{b}$ . Users have to choose the starting point for the approximate molecular width $\theta_{\sigma}$ . The sub-pixel refinement of the coordinates is obtained as $\hat{x}_{0}=\hat{\theta}_{x}$ and $\hat{y}_{0}=\hat{\theta}_{y}$ , where $\boldsymbol{\hat{\theta}}=\left[\hat{\theta}_{x},\hat{\theta}_{y},\ldots\right]$ .

Maximum-likelihood estimation

This approach assumes that the number of photons collected by a single camera pixel follows the Poisson distribution. Thus, the probability of $\kappa$ photons arriving at a camera pixel, where the expected number of photons is $\lambda$ , is given by

$p\left(\kappa\mid\lambda\right)=\frac{\lambda^{\kappa}\mathrm{e}^{-\lambda}}{% \kappa!}\,.$

(3)

Suppose that samples are drawn independently from the Poisson distribution, with the expected photon count $\lambda=\mathrm{PSF}\left(x,y\mid\boldsymbol{\theta}\right)$ given by the point-spread function model, and the observed photon count $\kappa=\tilde{I}\left(x,y\right)$ given by the image intensity expressed in photons. The likelihood [4, 5, 7, 3] of the parameters $\boldsymbol{\theta}$ can be modeled as

$L(\boldsymbol{\theta}\mid\mathcal{D})=\prod_{x,y\in\mathcal{D}}\frac{\mathrm{% PSF}\left(x,y\mid\boldsymbol{\theta}\right)^{\tilde{I}\left(x,y\right)}\mathrm% {e}^{-\mathrm{PSF}\left(x,y\mid\boldsymbol{\theta}\right)}}{\tilde{I}\left(x,y% \right)!}\,.$

(4)

The maximum likelihood estimate of the parameters $\boldsymbol{\theta}$ is, by definition, the value that maximizes the likelihood $L(\boldsymbol{\theta}\mid\mathcal{D})$ . Intuitively, the estimate $\boldsymbol{\hat{\theta}}$ corresponds to the value $\boldsymbol{\theta}$ that best agrees with the data. The maximization problem has the form

$\boldsymbol{\hat{\theta}}=\argmax_{\boldsymbol{\theta}}\sum_{x,y\in\mathcal{D}% }{\left[\tilde{I}\left(x,y\right)\ln{\left(\mathrm{PSF}\left(x,y\mid% \boldsymbol{\theta}\right)\right)}-\mathrm{PSF}\left(x,y\mid\boldsymbol{\theta% }\right)\right]}\,,$

(5)

which we solve by the Nelder-Mead method [6]. The starting point for the optimization process is computed from the data as the difference between the maximum and the minimum intensity values for the molecular intensity $\theta_{N}$ , and as the minimum intensity value for the background offset $\theta_{b}$ . Users have to choose the starting point for the approximate molecular width $\theta_{\sigma}$ . The sub-pixel refinement of the coordinates is obtained as $\hat{x}_{0}=\hat{\theta}_{x}$ and $\hat{y}_{0}=\hat{\theta}_{y}$ , where $\boldsymbol{\hat{\theta}}=\left[\hat{\theta}_{x},\hat{\theta}_{y},\ldots\right]$ .

Constraining parameters of PSF models

The Levenberg-Marquardt algorithm and the Nelder-Mead method used above search for values of the parameters $\boldsymbol{\theta}$ over an infinite interval. The optimization process can therefore converge to a solution with negative values which is impossible for variables corresponding to image intensity or to the standard deviation of a Gaussian PSF. We therefore limit the interval of possible values by transforming the relevant parameters and using $\mathrm{PSF}\left(x,y\mid\boldsymbol{\tilde{\theta}}\right)$ in Equations (2) and (5) instead of $\mathrm{PSF}\left(x,y\mid\boldsymbol{\theta}\right)$ . The transformation for a 2D Gaussian PSF model is $\boldsymbol{\tilde{\theta}}=\left[\theta_{x},\theta_{y},\theta_{\sigma}^{2},% \theta_{N}^{2},\theta_{b}^{2}\right]$ . The optimization process is still unconstrained but will result in positive PSF parameters. Such a transformation also improves the stability of the fit.

Guidelines for the choice of parameters

Fitting PSF models by maximum likelihood estimation generally gives slightly better results than fitting by least square methods, particularly when the photon counts are low. Maximum likelihood estimation requires correct setup of the camera parameters (photoelectrons per A/D count and the digitizer offset level). The recommended PSF model, which generally works well, is the integrated Gaussian. The initial size of sigma can be found by running ThunderSTORM on few images of the data sequence. A histogram of the fitted sizes of sigma (in pixels) can help to find the initial value. The size of the fitting radius should be an integer number close to 3*sigma.

References

[1] P. R. Bevington and D. K. Robinson(2003) Data reduction and error analysis for the physical science, McGraw-Hill Higher Education, McGraw-Hill. External Links: ISBN 9780072472271. Cited by: Least-squares methods.
[2] Commons-Math(2013-04) The Apache Commons Mathematics Library; version 3.2, External Links: Link. Cited by: Least-squares methods.
[3] F. Huang, S. L. Schwartz, J. M. Byars and K. A. Lidke(2011) Simultaneous multiple-emitter fitting for single molecule super-resolution imaging, Biomedical Optics Express 2 (5), pp. 1377–93. External Links: Document. Cited by: Maximum-likelihood estimation.
[4] M. Kendall and A. Stuart(1979) The Advanced Theory of Statistics, London: Charles Griffin. Cited by: Least-squares methods, Maximum-likelihood estimation.
[5] K. I. Mortensen, L. S. Churchman, J. A. Spudich and H. Flyvbjerg(2010) Optimized localization analysis for single-molecule tracking and super-resolution microscopy, Nature Methods 7 (5), pp. 377–381. External Links: Document. Cited by: Least-squares methods, Maximum-likelihood estimation.
[6] R. O’Neill(1971) Algorithm AS 47–function minimization using a simplex procedure, Applied Statistics 20, pp. 338–45. External Links: Link. Cited by: Maximum-likelihood estimation.
[7] C. S. Smith, N. Joseph, B. Rieger and K. A. Lidke(2010) Fast, single-molecule localization that achieves theoretically minimum uncertainty, Nature Methods 7 (5), pp. 373–5. External Links: Document, ISSN 1548-7105. Cited by: Maximum-likelihood estimation.

Generated on Mon Oct 27 10:48:41 2014 by LaTeXML [LOGO]