Fast and accurate threedimensional point spread function computation for fluorescence microscopy
Abstract
The pointspread function (PSF) plays a fundamental role in fluorescence microscopy. A realistic and accurately calculated PSF model can significantly improve the performance in 3D deconvolution microscopy and also the localization accuracy in singlemolecule microscopy. In this work, we propose a fast and accurate approximation of the GibsonLanni model, which has been shown to represent the PSF suitably under a variety of imaging conditions. We express the Kirchhoff’s integral in this model as a linear combination of rescaled Bessel functions, thus providing an integralfree way for the calculation. The explicit approximation error in terms of parameters is given numerically. Experiments demonstrate that the proposed approach results in a significantly smaller computational time compared with current stateoftheart techniques to achieve the same accuracy. This approach can also be extended to other microscopy PSF models.
Key points
 We express the integral in this model as a linear combination of rescaled Bessel functions, providing an integralfree way for the calculation.
 Experiments demonstrate that the proposed approach results in significantly smaller computational time compared with the quadrature approach (e.g. PSF Generator^{1}).
 This approach can also be extended to other microscopy PSF models.
The GibsonLanni Model^{2}
This model is based on a calculation of the optical path difference (OPD) between the design conditions and experimental conditions of the objective. It accounts for coverslips and other interfaces between the specimen and the objective.
$$ OPD(\rho,\mathrm{z}; z_p, {\boldsymbol \tau}) = \left(\mathrm{z} + t_i^* \right)\sqrt{n_i^2  (\mathrm{NA}\rho)^2}+ z_p\sqrt{n_s^2  (\mathrm{NA}\rho)^2}  t_i^\sqrt{(n_i^)^2  (\mathrm{NA}\rho)^2} + t_g\sqrt{n_g^2  (\mathrm{NA}\rho)^2}  t_g^\sqrt{(n_g^)^2  (\mathrm{NA}\rho)^2}, $$ where $\rho$ is the normalized radius in the focal plane, $\mathrm{z}$ is the axial coordinate of the focal plane, $z_p$ is the axial location of the pointsource in the specimen layer relative to the cover slip and $\mathbf{p}=(\mathrm{NA}, \mathrm{n}, \mathrm{t})$ is a parameter vector containing the physical parameters of the optical system: $\mathrm{NA}$ is the numerical aperture, $\mathrm{n}$ represents the refractive index and $\mathrm{t}$ is the thickness of individual layers.
The GibsonLanni model is expressed by $$ \mathrm{PSF}(\mathrm{r}, \mathrm{z}; z_p, \mathbf{p}) = \leftA\int_0^a e^{iW(\rho,\mathrm{z}; z_p, \mathbf{p})} J_0\left(k\mathrm{r}\mathrm{NA}\rho\right) \rho \mathrm{d}\rho\right^2 $$ where the phase term $W(\rho,\mathrm{z}; z_p, \mathbf{p})=k,OPD(\rho,\mathrm{z}; z_p, \mathbf{p})$, $k=2\pi/\lambda$ is the wave number of the emitted light. $A$ is a constant complex amplitude, and $J_0$ denotes the Bessel function of the first kind of order zero.
Bessel Series Approximation
The main idea is based on the fact that the integral $\int_{0}^{\alpha}tJ_0(ut)J_0(vt)dt$ can be explicitly computed as ^{3} $$ \int_0^a t J_0(ut)J_0(vt) dt = a\Big(\frac{uJ_1(ua)J_0(va)  v J_0(ua)J_1(va)}{u^2v^2}\Big). $$
If we expand the function $e^{iW(\rho,\mathrm{z}; z_p, \mathbf{p})}$ as a linear combination of rescaled Bessel function and fit the coefficients, then $$ \mathrm{PSF}_{\mathrm{app}}(\mathrm{r}, \mathrm{z}; z_p, \mathbf{p}) = \leftA\sum_{m=1}^{M}c_m (\mathrm{z}) R_m(\mathrm{r}; \mathbf{p})\right^2, \mathrm{where\ }R_m(\mathrm{r}; \mathbf{p})=\frac{\sigma_mJ_1(\sigma_ma)J_0(\eta a)a  \eta J_0(\sigma_ma)J_1(\eta a)a}{\sigma_m^2  \eta^2}. $$
Compared with PSFGenerator
Software
Call the PSF generation from the ImageJOps
See more here.


Extensions
Noncylindrically symmetric PSF generation (e.g. with coma, astigmatism aberrations represented in vectorial models)?
In principle, the proposed approach can be extended to those cases, provided other models are able to be expressed in the form of Kirchhoff’s integral formula. Haeberle^{4} showed that the vectorial model can also be combined with the ease of use of the GibsonLanni scalar approach, which has the advantage of introducing explicitly the known or sampledependent parameters.
Here is a sample code which supports angledependent PSF generation.
PSF Engineering (e.g. doublehelix PSF)?
It maybe interesting in some sense to include the spherical aberration in the doublehelix PSF, see^{5} for the discussion.
Applications
See citations, more examples to be added.
 3D deconvolution microscopy by PURELET and 3D microscopy PSF estimation
 Cytokit: A singlecell analysis toolkit for high dimensional fluorescent microscopy imaging
 Flowdec: TensorFlow Deconvolution for Microscopy Data
References
H. Kirshner, F. Aguet, D. Sage, and M. Unser, “3D PSF fitting for fluorescence microscopy: implementation and localization application,” J. Microsc. vol. 249, no. 1, pp. 13–25, 2013. ↩︎
S. F. Gibson and F. Lanni, “Experimental test of an analytical model of aberration in an oilimmersion objective lens used in threedimensional light microscopy”, J. Opt. Soc. Am. A, vol. 9, no. 1, pp. 154166, 1992. ↩︎
G. N. Watson, “A treatise on the theory of Bessel functions”, Cambridge University Press, 1995. ↩︎
O. Haeberle, “Focusing of light through a stratified medium: a practical approach for computing microscope point spread functions. Part I: conventional microscopy”, Opt. Commun., vol. 216, no. 1, pp. 5563, 2003. ↩︎
S. Ghosh and C. Preza, “Characterization of a threedimensional doublehelix pointspread function for fluorescence microscopy in the presence of spherical aberration”, J. Biomed. Opt., vol. 18, no. 3, pp. 036010110, 2013. ↩︎