Fast and Accurate Three-Dimensional Point Spread Function Computation for Fluorescence Microscopy

A realistic and accurately calculated PSF model can significantly improve the performance in deconvolution microscopy and also the localization accuracy in single-molecule microscopy. In this work [1], we propose a fast and accurate approximation to the Gibson-Lanni model [2].

Key points:

  • We express the integral in this model as a linear combination of rescaled Bessel functions, providing an integral-free 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 [3]).

  • This approach can also be extended to other microscopy PSF models.

The Gibson-Lanni Model

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 point-source 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 Gibson-Lanni model is expressed by

\[ \mathrm{PSF}(\mathrm{r}, \mathrm{z}; z_p, \mathbf{p}) = \left|A\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 [4]

\[ \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^2-v^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}) = \left|A\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


Figure 1. (a) One example PSF (\(256 \times 256 \times 128\)). (b) Comparison of computational time with PSF Generator [3] ('Best’ option) for a variety of image size. The approximation parameters \(M\) and \(K\) in the proposed approach are chosen to result in the same accuracy.


(a) Matlab
(b) ImageJ Plugin
(C) Icy Plugin
(d) Python
(by Dr. Kyle Douglass)

Call the PSF generation from the ImageJ-Ops, see more.
        import net.imagej.ops.AbstractOpTest;
        import net.imglib2.Dimensions;
        import net.imglib2.FinalDimensions;
        import net.imglib2.img.Img;
        import net.imglib2.type.numeric.real.DoubleType;

        double NA = 1.4; // numerical aperture
        double lambda = 610E-09; // wavelength
        double ns = 1.33; // specimen refractive index
        double ni = 1.5; // immersion refractive index, experimental
        double resLateral = 100E-9; // lateral pixel size
        double resAxial = 250E-9; // axial pixel size
        double pZ = 2000E-9D; // position of particle

        DoubleType type = new DoubleType(); // pixel type of created kernel
        Dimensions dims = new FinalDimensions(256, 256); //dimensions

        Img<DoubleType> kernel = ops.create().kernelDiffraction(dims, NA, lambda, ns, ni, resLateral, resAxial, pZ, type);


  • [1] J. Li, F. Xue and T. Blu, “Fast and accurate three-dimensional point spread function computation for fluorescence microscopy”, J. Opt. Soc. Am. A, vol. 34, no. 6, pp. 1029-1034, 2017. [PDF][Slides][Matlab Code][ImageJ Plugin][Icy Plugin][Python Code] (Top Downloaded Article in Oct 2017)

  • [2] 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, vol. 9, no. 1, pp. 154-166, 1992. [Link]

  • [3] 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. [Link]

  • [4] G. N. Watson, “A treatise on the theory of Bessel functions”, Cambridge University Press, 1995. [Link]