## Abstract

A proper estimation of realistic point-spread function (PSF) in optical microscopy can significantly improve the deconvolution performance and assist the microscope calibration process. In this work, by exemplifying 3D wide-field fluorescence microscopy, we propose an approach for estimating the spherically aberrated PSF of a microscope, directly from the observed samples. The PSF, expressed as a linear combination of 4 basis functions, is obtained directly from the acquired image by minimizing a novel criterion, which is derived from the noise statistics in the microscope. We demonstrate the effectiveness of the PSF approximation model and of our estimation method using both simulations and real experiments that were carried out on quantum dots. The principle of our PSF estimation approach is sufficiently flexible to be generalized non-spherical aberrations and other microscope modalities.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The point-spread function (PSF) describes the response of the imaging system to a point source or object. It plays a fundamental role in understanding imaging performance, such as the theoretical resolution limit and the optical sectioning capacity, and to identify any problems with the microscope settings for calibration [1, 2]. Thus a realistic and accurate knowledge of the PSF is crucial to optimize the performance of a microscope and for successful deconvolution, allowing for highly accurate reconstruction of biological structures [3–5].

Depending on the imaging modalities (wide-field, confocal, light-sheet etc), PSFs have different shapes. As the exact PSF in the microscope is rarely known, one has to rely on an approximation, which can be accessed using an experimental or an analytical approach. The classical experimental method is to image sub-diffraction sized fluorescent microspheres or quantum dots (<200 nm diameter) [1,6], which contains both the intrinsic and extrinsic aberrations [7]. However, this process often gets tedious, since several measurements need to be averaged over user selected ROIs due to noise degradations. Directly denoising them may cause loss in essential information, and their resolution are tied to the resolution of the acquisition. Moreover, in the case of deconvolution microscopy, the imaging conditions of experimental PSFs may be different from those in actual imaging, which makes the restoration performance suboptimal. For instance, the PSF in Stimulated Emission Depletion (STED) microscopy depends strictly on the properties of the fluorescence dye, while in most of the cases the dye for microspheres is different from the one used for labeling the sample.

Several approaches have been proposed to address the question how to accurately model analytically the PSF in accordance with the laws of optics using diffraction theory and knowledge of the optical components [8–10]. The advantage is that the generated model can be evaluated for the entire object space. Because real optical systems are prey to a variety of aberrations, the PSF does not obey some simplified models. This is especially apparent when using lenses with high numerical aperture. Extensive PSF models, such as the Richards-Wolf model [11] and the Gibson-Lanni model [8], are characterized by a large number of approximations with varying degrees of accuracy. In addition, despite the widespread PSF approximation by a Gaussian in single molecule localization microscopy, it has been argued that a more realistic model can significantly improve the localization accuracy [12–15].

Unfortunately, some of the model parameters such as the refractive index of the specimen, are difficult to obtain or might change due to heating of live samples during the course of experiments. Thus the resulting PSF is unlikely to match experimental conditions. Moreover, in some microscopy techniques, for example Selective Plane Illumination Microscopy (SPIM), the PSF is object-dependent (the product of the widefield detection PSF and the intensity profile of the excitation light sheet along the axial direction) [16], which makes it impossible to compute the PSF.

It is thus preferable to use a more global (i.e. object-dependent) approach for estimating a PSF, directly from the acquired data. Most algorithms that choose this option perform the simultaneous estimation of the PSF and of the object, based on *prior* hypotheses on them under a maximum a posteriori (MAP) framework (blind deconvolution) [17–19]. Alternatively, based on an analytic model and a carefully designed optimization criterion, one can estimate the PSF by finding its parameters first [20–23], then perform deconvolution to restore the image of the sample. It is this approach that we are going to follow, because it allows to apply higher quality non-blind image deconvolution techniques [23,24].

Aiming at 3D wide-field fluorescence microscopy, we propose a blind PSF estimation algorithm which does not require to measure any additional materials (e.g. fluorescence beads or quantum dots). Inspired by the commonly-used Gibson-Lanni model [8], our PSF model which considers only spherical aberrations consists of a linear combination of 4 basis functions. These basis functions are also derived from Kirchhoff’s integral formula, thus automatically satisfy all the constraints imposed by the optical model. Then the PSF is obtained by minimizing a novel optimization criterion, the “blur-PURE” (Poisson Unbiased Risk Estimate), which is based on the realistic noise statistics in fluorescence microscopy [25]. Compared with a non-parametric representation (i.e., from an implicit regularization framework), a parametric formulation requires much less values to estimate, and dramatically reduces the degrees of freedom in the problem. In addition, the estimated PSF is defined in the continuous space, thus can be used to restore different-size acquisitions under the same imaging settings.

The paper is organized as follows. We firstly describe the 4-basis PSF approximation model for 3D wide-field fluorescence microscopy in Section 2. Section 3 introduces the blur-PURE criterion and our PSF estimation method. Section 4 and 5 demonstrate the effectiveness of the proposed PSF approximation model and of our estimation method on both synthetic and experimental acquired data.

## 2. 4-basis approximation of a wide-field microscopy PSF

Wide-field microscopy is still the foundation of some advanced techniques, such as multi-photon microscopy and super resolution localization microscopies. It is widely used to capture the 3D structures of living biological samples by collecting a stack of 2D images. While this wide-field recording has the advantage of fast acquisition and low light exposure, the 3D image is always blurred by the contribution of light from out-of-focus planes. Meanwhile, the PSF of wide-field microscopy is axially asymmetric due to the mismatch between refractive indices of the immersion medium and of the specimen [8].

Some theoretically derived models of 3D wide-field microscopy that assume the use of an aberration-free objective lens have been shown experimentally to be inaccurate [26]. Thus the aberrations introduced when the microscope is used under non-design optical conditions must be incorporated into the theoretical model. The Gibson-Lanni model [8] is widely used since it can predict the non-symmetric patterns in the axial direction, which reflects the spherical aberration and often appears in realistic imaging conditions. This model is based on a calculation of the optical path difference (OPD) between experimental conditions and the design conditions of the objective. It has been shown to be very useful for deconvolution microscopy [5,27,28] and also for particle localization [29–33]. However, for the PSF estimation purpose, the Gibson-Lanni model depends on a large number of non-linear parameters, which renders it inadequate for a global optimization approach.

Inspired by the Gibson-Lanni model, we propose to approximate the spherically aberrated PSF by a linear combination of 4 basis functions. For each basis function, the key idea is to reduce Gibson-Lanni’s OPD expression (made of a sum of five square roots) to a sum of two square roots only:

where*z*is the axial coordinate of the focal plane, and

*ρ*is the normalized radius in the focal plane. The (

*p*,

_{m}*q*) are akin to refractive indices, and their values are fixed once for all (see below). The parameter

_{m}*η*is an indication of the focus position

*z*. More specifically, the PSF

_{p}*h*is written as

*x*,

*y*,

*z*) is the 3D spatial coordinate,

*c*,

_{m}*m*= 1, 2, 3, 4} are the coefficients of the basis functions, (

*p*

_{1},

*q*

_{1}) = (1.35, 1.35), (

*p*

_{2},

*q*

_{2}) = (1.35, 1.45), (

*p*

_{3},

*q*

_{3}) = (1.45, 1.35), (

*p*

_{4},

*q*

_{4}) = (1.45, 1.45),

*A*is a constant complex amplitude,

*K*= 2

*π*/

*λ*is the wave number,

*r*= ‖r‖ and

*J*

_{0}denotes the Bessel function of the first kind of order zero. Each basis function

*h*(r;

_{m}*η*) is computed efficiently using a fast approach based on Bessel series approximation in [34]. We demonstrate (see Section 4) that the 4-basis model in Eq. (2) approximates any Gibson-Lanni PSF very accurately, provided that the range of refractive index of the specimen are within [1.345, 1.5].

Intuitively, the OPD expression in the original Gibson-Lanni model seems to be over parametrized: it is essentially made of a sum of terms like $\alpha \sqrt{{\beta}^{2}-{\rho}^{2}}$ which, when *ρ* ≪ *β*, are equivalent to *α′* − *β′ρ*^{2}. Hence, a sum of two square roots as in Eq. (1) should be sufficient to approximate accurately the OPD in the Gibson-Lanni model.

For the proposed parametrization, when the wavelength *λ* and NA are provided (which is usually the case), only the focus indicator *η* and 4 linear coefficients c have to be determined. The degrees of freedom is greatly reduced compared with the original Gibson-Lanni model. Note that the linear approximation in Eq. (2) is different from the approaches in Markham et al. [20] and Soulez et al. [22], where the phase term of the original model is approximated by a linear combination of several (typically larger than 6) polynomials (either power or Zernike polynomials). The corresponding coefficients are thus highly non-linear and more difficult to estimate accurately since they are involved in an integral. The proposed parametrization is also distinct from the approach in [18], where the basis functions are learned from a training set of PSFs with only different focus positions.

By using a different set of basis functions, our strategy can be generalized to non-spherical aberrations (e.g. coma, astigmatism) that are correctly modelled in vectorial models [10,11,35]. In fact, Haeberlé [36] showed that the vectorial model can also be combined with the ease of use of the Gibson-Lanni scalar approach, which has the advantage of introducing explicitly the known or sample-dependent parameters [37].

## 3. PSF estimation method

#### 3.1. Imaging model

The optical microscope can be modeled as a linear shift-invariant system, due to the incoherent nature of the emitted light [3,4]. Such a system is characterized by a PSF *h*(r), where r = (*x*, *y*, *z*) is the 3D spatial coordinate. We do not consider space-varying PSFs [38–40] in the current work. The noise degradation during image acquisition is considered to be mixed Poisson-Gaussian statistics [41,42], which accounts for a broad range of imaging situations (from photon-limited to sensor-limited imaging).

*denotes the distorted observation of the unknown true original image X ∈ ℝ*

^{N}*,*

^{N}*N*=

*N*×

_{x}*N*×

_{y}*N*is the size of the acquired 3D image. We assume independence of the individual components of the random variable Y and that of the photon-counting process and the read-out noise. The Scale

_{z}*α*represents the gain of the acquisition device, which controls the noise during the acquisition process.

**H**

_{0}∈ ℝ

^{N×N}is a block-circulant matrix, which implements a discrete convolution with the PSF h.

**I**denotes the identity matrix, and

*𝒫*(·) and

*𝒩*(·, ·) represent the effect of the Poisson noise and the additive Gaussian noise (variance

*σ*

^{2}), respectively. The values of

*α*and

*σ*

^{2}can be estimated by a robust linear regression performed on a collection of local estimates of the sample mean and variance [42,43].

Based on the image acquisition model in Eq. (3), the objective of this work is to estimate the PSF *h* (namely **H**_{0}) directly from the measurement Y. This is a well-known difficult inverse problem since both the original image X and the PSF are unknown. However, as stated in the introduction part, the PSF approximation with few involved parameters can greatly reduces the degrees of freedom [20,22–24].

#### 3.2. Blur-PURE as an optimization criterion

When the original image X is known, the estimation of **H**_{0} can be done by minimizing the expected mean squared error (MSE) between the blurred ground-truth image **H**_{0}X and a linear processing **U _{H}** of the acquired image Y [23,24]: $\frac{1}{N}\mathcal{E}\left\{{\Vert {\mathbf{U}}_{\mathbf{H}}\text{Y}-{\mathbf{H}}_{0}\text{X}\Vert}^{2}\right\}$. This oracle criterion (

**H**

_{0}X is assumed to be known) is named “

*blur-MSE*”, since it measures the closeness to the blurred ground-truth image. A good choice for

**U**is Tikhonov regularized deconvolution:

_{H}**U**

_{H,μ}=

**HW**, where

_{H}**W**=

_{H}**H**

^{T}(

**H**

^{T}

**H**+

*μ*

**P)**

^{−1},

**P**is an approximation of the power density spectrum of the origin image X and

*μ*is some positive scalar. Similar to [23], it can be proven that for all

**U**

_{H,μ}, the solution

**H**minimizing the

*blur-MSE*is related to the true matrix

**H**

_{0}, through ${\mathbf{H}}^{\text{T}}\mathbf{H}={\mathbf{H}}_{0}^{\text{T}}{\mathbf{H}}_{0}$, i.e., ‖

*ĥ*(

*ω*)‖

^{2}= ‖

*ĥ*

_{0}(

*ω*)‖

^{2}in the Fourier domain. Since

*ĥ*and

*ĥ*

_{0}satisfy a very low order parametric representation, this equation is equivalent to

*h*(r) =

*h*

_{0}(r − r

_{0}) where r

_{0}is some arbitrary shift.

The matrix **U**_{H,μ} corresponds to a filter whose frequency response can be thought as a kind of band-indicator since it marks a certain frequency band as 0 or 1 with a narrow transition between the two values. The exact band-indicator **U**_{0}(*ω*) = |**H**_{0}(*ω*)|^{2} (|**H**_{0}(*ω*)|^{2} + *C*(*ω*)/*A*(*ω*))^{−1}, where *A*(*ω*) and *C*(*ω*) are the power spectral densities of signal X and noise, respectively. The **P** in **U**_{H,μ} is often expressed by the discrete Laplacian operator (${\Vert \omega \Vert}^{2}={\omega}_{x}^{2}+{\omega}_{y}^{2}+\gamma {\omega}_{z}^{2}$) so that *C*(*ω*)/*A*(*ω*) ≈ *μ*‖*ω*‖^{2} and **U**_{H,μ} thus serves as an approximation to the exact band-indicator, where *γ* is the ratio between lateral and axial resolutions.

In practice, the *blur-MSE* cannot be minimized directly since X is unknown. However, without any assumptions on the noise-free data, the quantity of *blur-MSE* can be replaced by an statistical estimate, blur-PURE (Poisson Unbiased Risk Estimate), which involves the measurement Y only. For the linear degradation model in Eq. (3), we have the following theorem (similar to [23,24,44]).

**Theorem 1** *Consider the degradation model in* *Eq.* (3) *and* **U** *an arbitrary matrix, the random variable blur-PURE is an unbiased estimate of the blur-MSE; i.e.,*

*where*e

_{n}is the N-dimensional vector with components δ_{k−n},

*k*= 1, 2, ...,

*N, and N is the number of pixels of the image.*

This theorem is the consequence of probabilistic relation that satisfied for Poisson-Gaussian random variables: with the notation of Eq. (3), we have that (see Eq. (7) in [42]) $\mathcal{E}\left\{{\text{X}}^{\text{T}}{\mathbf{H}}_{0}^{\text{T}}\mathbf{U}\text{Y}\right\}=\mathcal{E}\left\{{\text{Y}}^{\text{T}}\mathbf{U}\left(\text{Y}-\alpha {\text{e}}_{n}\right)-{\sigma}^{2}\text{Tr}(\mathbf{U})\right\}$, where *ℰ*{·} denotes the mathematical expectation operator. The proof of Theorem 1 is very similar to the one in [23], the main difference being that, in [23], only Gaussian statistics was considered (i.e. $\mathcal{E}\left\{{\text{X}}^{\text{T}}{\mathbf{H}}_{0}^{\text{T}}\mathbf{U}\text{Y}\right\}=\mathcal{E}\left\{{\text{Y}}^{\text{T}}\mathbf{U}\text{Y}-{\sigma}^{2}\text{Tr}(\mathbf{U})\right\}$).

This criterion solely depends on the observed image Y thus is computable. The statistical unbiasedness with the *blur-MSE* and the fact that the pixel number of the image *N* is very large (typically for a 3D image 256 × 256 × 32, *N* > 2 × 10^{6}) indicates that the blur-PURE can be used as a reliable subsitute of the *blur-MSE* (law of large numbers). See Fig. 1(a) for a typical example of the closeness between the *blur-MSE* and blur-PURE. The maximum difference for all cases is 3.74 × 10^{−4}. Note that Theorem 1 is valid for any linear distortion **H**; i.e., is not limited to convolutions. Here, since **U**_{H,μ} is a convolution, all matrices involved are diagonalized by the 3D discrete Fourier transformation, and so, the blur-PURE can be efficiently computed in the Fourier domain.

#### 3.3. PSF estimation by band-indicator approximation

The optimization criterion, blur-PURE, provides explicit control over the PSF estimation problem. However, directly minimizing the blur-PURE cannot guarantee the optimal solution, which is mainly because this is a highly non-convex optimization problem that has many local minima. Instead, we propose to firstly approximate the true band-indicator **U**_{H0,μ} by a linear combination of basis functions ${\mathbf{U}}_{\text{app}}={\sum}_{n}{a}_{n}{\mathbf{U}}_{{\mathbf{H}}_{n},\mu}$. Here **H*** _{n}*’s are constructed by the corresponding PSF basis

*h*(r;

_{n}*η*), which are described in Section 2. The accuracy of this approximation is shown in Fig. 1(b). Thanks to the quadratic nature of the blur-PURE, the search for the optimal coefficients a = {

*a*,

_{n}*n*= 1, 2, 3, 4} boils down to the solution of a linear system of equations. Auxiliary parameters

*η*and

*λ*are estimated during this stage by a derivative-free optimization technique. Specifically, we used the Nelder-Mead Simplex method with bound constraints to find the optimal values. The maximum iteration number is set to be 400, and the starting points for

*η*and

*λ*are 0 and 10

^{−3}(

*α*Y

_{mean}+

*σ*

^{2}), respectively.

Once **U**_{app} has been found, finding the PSF **H** amounts to solving the following unconstrained minimization problem: $\underset{\mathbf{H}}{\text{min}}{\Vert {\mathbf{HW}}_{\mathbf{H}}-{\mathbf{U}}_{\text{app}}\Vert}_{2}^{2}$, where **H** assumes the parametric expression in Eq. (2). We propose an iterative optimization algorithm:

**H**

^{(0)}is randomly initialized. The solution

**H**

^{(k)}involves finding the coefficients ${c}_{m}^{(k)}$’s by solving a linear system of equations at each iteration until ‖

**H**

^{(k)}−

**H**

^{(k−1)}‖/‖

**H**

^{(k−1)}‖ ≤ 10

^{−3}. The convergence to stationary points of the objective function is achieved generally within 50 iterations. The scheme diagram describing the proposed approach is shown in Fig. 2.

## 4. Evaluation of the PSF approximation model

#### 4.1. 4-basis model validation

The PSF approximation error is measured in terms of the relative squared error (RSE) [34,45] calculated as

*h*

_{0}is generated based on the complete Gibson-Lanni model under different acquisition conditions and normalized to have a maximum value of one [34]. The approximated PSF

*h*

_{app}is obtained by fitting the 4-basis model in Eq. (2) to Gibson-Lanni “ground-truth” from various settings. Basically, the estimation is thought to be accurate when RSE < 9% [45]. Calculations are carried out on a Macbook Pro with a 2.8 GHz Intel Core i7, with 16 GB of RAM.

### Fitting procedure

To find the optimal parameters (*η* and c in Eq. (2)) such that the approximated PSF can be as close as possible to the ground truth, we perform an exhaustive search over the possible values of *η* within the thickness of the sample. For each candidate of *η*, those 4 coefficients c are determined from the least square fitting by solving a linear system of equations. The one corresponding to the smallest value of RSE between the reconstructed and ground truth PSFs is selected. Fig. 3(a) shows a typical example of the fitting process corresponding to different focus positions (*z _{p}* = 0, 500 nm, 1000 nm). The estimated

*η*is 0, 529.10 nm, and 1025.07 nm respectively.

In this work, we focus on the PSF estimation instead of some specific parameters. As stated before, this *η* is an indicator of focus position but is not equal to *z _{p}*. Fig. 3(b) shows the scatter plot between the ground truth

*z*and the estimated

_{p}*η*. The mean and max difference between them are 25.29 nm and 46.00 nm respectively. As we can see from Fig. 3(b) and other experiments (not shown here), the relation between

*η*and

*z*is deterministic (function variables may depend on wavelength

_{p}*λ*, NA etc). This may be caused by the focal shift occurred in optical microscopy due to the refractive index mismatch [46]. This suggests that it is possible to tabulate the relation between

*η*and

*z*, which may have some interest beyond the scope of this paper.

_{p}### Noise-free conditions

We consider the wavelength *λ* used in the range from 340 nm to 750 nm with a step of 10 nm, as usual in real experiments. The refractive indices *n _{s}*’s of typical cellular components are in the range from 1.354 to 1.5 with a step of 0.02. Other parameters are consistent with the setting of a 100× magnification, 1.4 NA oil immersion objective (the refractive index

*n*= 1.515). The focus position

_{i}*z*varies from −1.24

_{p}*μ*m to 1.24

*μ*m with a step of 0.05

*μ*m, which covers the range of expected spherical aberrations. There are totally 20, 496 PSFs of size 127 × 127 × 63 generated.

The mean and standard deviation of the RSE are 0.26% and 0.42% respectively. Fig. 4(a) shows the mean RSEs with respect to the wavelength. Note that PSFs with longer wavelength have lower frequency (less rings) thus are more accurately approximated. For a typical PSF (*z _{p}* = 500 nm,

*n*= 1.394 and

_{s}*λ*= 395 nm), the

*z*-profiles of the true and fitted PSFs are shown in Fig. 4(b). Apart from the low fitting errors, the non-symmetric pattern caused by the refractive index mismatch can be well predicted by the proposed model. These results indicate that the 4-basis approximation model fits the Gibson-Lanni model very well, thus can be used as a good alternative with much fewer parameters, albeit with less straightforward physical interpretation.

### Noisy conditions

To evaluate the fitting performance in noisy conditions, we further degrade the above typical PSF with two noise levels (low noise: *α* = 0.02, *σ* = 0 and high noise: *α* = 0.2, *σ* = 0.02). Then the approximation model in Eq. (2) is again used to fit these noisy PSFs. The ground truth, noisy and fitted PSFs are shown in Figs. 4(b)–4(f), which shows that the model in Eq. (2) provides a robust approximation to the spherically aberrated PSFs even in noisy conditions.

#### 4.2. Experimental PSF fitting

For experimental validation, we prepared samples by drying dilutions of Qdot^{TM} 605 ITK streptavidin conjugated quantum dots (Qdots) on to a slide. These Qdots including the coating have an approximate diameter of 20 nm. A small drop of the mounting medium was added and a coverslip was placed to cover the dried Qdots on the slide. Fluorescence was collected by a 100×, 1.49 NA oil-immersion objective (Nikon) with an sCMOS camera (Andor Zyla VSC-03278). The excitation and emission peaks of these Qdots are 375 nm and 605 nm, respectively. The physical pixel size of the camera mounted on the microscope is 65.29 nm, and the measurement of the 3D stack is performed with a step size of 89.83nm, whose value is suggested in [1] to ensure the proper sampling frequency. Complete 3D stacks of the sample are acquired, and PSFs are extracted by the PSFj software [2] from the dataset. The average PSF is taken as the measured PSF.

The *x*–*y* and *y*–*z* sections of such a PSF (55 × 55 × 113) and its fitted result is shown in Fig. 5. It shows that the proposed approximation model fits the data well over the complete sample space, even for current relatively difficult condition (high NA, obvious asymmetric pattern). Fitting a Gaussian function will certainly fail in this case to represent the long tails of the measured PSF. The estimated focus indicator *η* in Eq. (2) is 1.855 *μ*m. We can also see a substantial suppression of the extraneous noise in the fitted PSF. These results illustrate the appropriateness of this 4-basis model in Eq. (2) for microscopy PSF in wide-field fluorescence microscopy.

## 5. PSF estimation results

#### 5.1. Validation on simulated data

We adapt the confocal image simulator in [47] to the wide-field settings. Two PSFs with different focus positions (*z _{p}* = 0

*μ*m, 2

*μ*m respectively) are generated [34] according to the Gibson-Lanni model with typical values for the parameters (

*λ*= 622nm, NA=1.4,

*n*= 1.33,

_{i}*n*= 1.46). They are corresponding to the cases without and with the spherical aberration respectively. Other parameters are consistent with the setting of a Nikon Apo microscope. The pixel size is 100 nm in the

_{s}*x*–

*y*plane and 250 nm along the

*z*-axis. The Bars image (available at http://bigwww.epfl.ch/deconvolution) is convolved with these PSFs. The blurred images are subsequently contaminated by mixed Poisson-Gaussian noise with different noise levels (corresponding to different

*α*and

*σ*values). Typically, we obtain a low noise image when

*α*= 0.02,

*σ*= 0.02 and a high noise image when

*α*= 0.2,

*σ*= 0.2.

We compare the proposed approach with a popular blind deconvolution method (the blind Richardson-Lucy algorithm [48]) and the EpiDEMIC plugin in Icy (http://icy.bioimageanalysis.org/plugin/EpiDEMIC). The latter uses a PSF parametrization by means of decomposition of the pupil on Zernike basis, and a continuous optimization method to iteratively estimate both PSF parameters and the object. Similar to the present method, the PSF parametrization requires only some general parameters (wavelength *λ*, NA and the refractive index of the immersion medium *n _{i}*).

Table 1 presents the RSEs between the estimated and ground-truth PSFs over different scenarios. Note that all estimated PSFs are optimally shifted in the *z*-axis to best match the ground-truth PSF. It can be seen that the proposed approach generally yields significantly more accurate and consistent PSF estimation than other methods. Two typical cases are shown in Fig. 6. The proposed approach succeeds in estimating the spherical aberration of the PSF. The iterative algorithm in Eq. (4) generally converges within less than 50 iterations. The proposed method takes around 60 seconds generally, the Blind-RL takes around 70 seconds and EpiDEMIC takes more than 120 seconds for the estimation. Fast computation is particularly beneficial in localization microscopy and, generally, for high-throughput image restoration. Note that we do not estimate the physical parameters (refractive indices, working distance etc), but only the PSF that they generate. It may be interesting, though, to investigate how to retrieve some of these physical parameters from the estimated PSF.

#### 5.2. Validation on real data

In practice, Qdots in the acquired images are not always well-separated due to aggregation and sedimentation. In this case, we can estimate the PSF from the acquired image using our algorithm (described in Section 3), and then compare the estimated PSF with the measured one from isolated Qdots under the same setting.

The *x*–*y* and *z*–*y* sections of the estimated PSF from the real wide-field image shown in Fig. 8(a) are shown in Figs. 7(a) and 7(b) respectively. Fig. 7(c) shows the *z* profile of the estimated PSF as well as the Qdot experimental (“measured”) PSF and its approximation (“fitted”) using the parametric expression in Eq. (2). As can be seen from this figure, these three profiles show similar signal distributions, including their sidelobes, which indicates that the proposed method can make an accurate estimation of the PSF in real conditions. Note that the PSF fitting is performed on a single measured PSF (described in Section 4.2), while the estimation is done for the whole image. As expected, under the same imaging setting, the difference between estimated and fitted PSFs is very small (RSE = 2.12%).

By using the experimental or estimated PSF with our approach, we can perform deconvolution on the image to retrieve the contrast of Qdots. We firstly apply the popular Richardson-Lucy algorithm and our recent 3D PURE-LET algorithm [5,42] with the experimental PSF. Figs. 8(a)–8(c) show the blurred noisy observation as well as the deconvolved images with these two approaches. With the experimental PSF, the Richardson-Lucy algorithm fails to reconcentrate the spreading pattern. In addition, both approaches produce some unpredictable structures or ring artifacts, which is mainly due to the inevitable noise of the measured PSF.

We further obtained the restored image using the PURE-LET algorithm with the PSF estimated by the proposed estimation method (described in Section 3), as shown in Fig. 8(f). As baselines, the blind-RL algorithm and the EpiDEMIC plugin are applied and the restored images are shown in Figs. 8(d) and 8(e). Compared with these two methods, the PURE-LET significantly enhances image resolution as many dots details can be clearly resolved in the deconvolved image (see zoomed-in regions). The restored image has clearly removed the blur and suppressed the noise. In particular, the elongation phenomenon and asymmetric pattern along the optical axis due to the PSF are effectively eliminated. These results demonstrate that a properly estimated PSF by the proposed approach is an essential ingredient, in addition to an efficient restoration algorithm, to effectively improve the spatial resolution in wide-field fluorescence microscopy.

## 6. Conclusion

We have proposed a blind PSF estimation approach dedicated to wide-field fluorescence microscopy. The method consists in expressing the PSF as a linear combination of four basis functions, which makes it depend only on five parameters (4 linear and 1 non-linear). In order to estimate these parameters, we have introduced a novel criterion based on the mixed Poisson-Gaussian noise statistics (namely blur-PURE) and an efficient iterative algorithm to minimize this criterion. We illustrate the effectiveness of proposed approach using both simulations and experimental measured data, and demonstrate that the estimated PSF can be used to effectively deconvolve 3D wide-field fluorescence microscopy data. The flexibility of our approach makes it suitable for generalization to other microscope modalities.

Although our parametrization takes care of spherical aberrations only, preliminary tests show that our algorithm is quite robust to other types of aberration (coma and astigmatism), as long as they remain mild (e.g. 30m*λ*): the result obtained is a close, spherically aberrated approximation of the true PSF. Moreover, using this approximation in a state-of-the-art restoration algorithm, results in negligible quality loss. Our future plans are to incorporate more aberrations into the PSF approximation model, and apply the PSF estimation technique to microscope calibration [49,50].

## Funding

Research Grants Council of Hong Kong (CUHK14210617, CUHK14201317, AoE/M-05/12); National Natural Science Foundation of China (61401013).

## References

**1. **R. W. Cole, T. Jinadasa, and C. M. Brown, “Measuring and interpreting point spread functions to determine confocal microscope resolution and ensure quality control,” Nat. Protoc. **6**, 1929–1941 (2011). [CrossRef] [PubMed]

**2. **P. Theer, C. Mongis, and M. Knop, “PSFj: know your fluorescence microscope,” Nat. Methods **11**, 981–982 (2014). [CrossRef] [PubMed]

**3. **P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Process. Mag. **23**, 32–45 (2006). [CrossRef]

**4. **D. Sage, L. Donati, F. Soulez, D. Fortun, G. Schmit, A. Seitz, R. Guiet, C. Vonesch, and M. Unser, “Deconvolutionlab2: An open-source software for deconvolution microscopy,” Methods **115**, 28–41 (2017). [CrossRef] [PubMed]

**5. **J. Li, F. Luisier, and T. Blu, “PURE-LET deconvolution of 3D fluorescence microscopy images,” in *Proceedings of IEEE International Symposium on Biomedical Imaging* (IEEE, 2017), pp. 723–727.

**6. **W. Wallace, L. H. Schaefer, and J. R. Swedlow, “A workingperson’s guide to deconvolution in light microscopy,” Biotechniques **31**, 1076–1097 (2001). [CrossRef]

**7. **J. W. Shaevitz and D. A. Fletcher, “Enhanced three-dimensional deconvolution microscopy using a measured depth-varying point-spread function,” J. Opt. Soc. Am. A **24**, 2622–2627 (2007). [CrossRef]

**8. **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 **9**, 154–166 (1992). [CrossRef] [PubMed]

**9. **M. Born and E. Wolf, *Principles of optics: electromagnetic theory of propagation, interference and diffraction of light* (Cambridge University, 1999). [CrossRef]

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

**11. **B. Richards and E. Wolf, “Electromagnetic diffraction in optical systems II. structure of the image field in an aplanatic system,” P. Roy. Soc. Lond. A Mat. **253**, 358–379 (1959). [CrossRef]

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

**13. **D. Sage, H. Kirshner, T. Pengo, N. Stuurman, J. Min, S. Manley, and M. Unser, “Quantitative evaluation of software packages for single-molecule localization microscopy,” Nat. Methods **12**, 717–724 (2015). [CrossRef] [PubMed]

**14. **Y. Li, M. Mund, P. Hoess, U. Matti, B. Nijmeijer, V. J. Sabinina, J. Ellenberg, I. Schoen, and J. Ries, “Fast, robust and precise 3d localization for arbitrary point spread functions,” bioRxiv 172643; doi: https://doi.org/10.1101/172643.

**15. **H. P. Babcock and X. Zhuang, “Analyzing single molecule localization microscopy data using cubic splines,” Sci. Rep. **7**, 552 (2017). [CrossRef] [PubMed]

**16. **L. Gao, L. Shao, B.-C. Chen, and E. Betzig, “3D live fluorescence imaging of cellular dynamics using Bessel beam plane illumination microscopy,” Nat. Protoc. **9**, 1083–1101 (2014). [CrossRef] [PubMed]

**17. **E. F. Y. Hom, F. Marchis, T. K. Lee, S. Haase, D. A. Agard, and J. W. Sedat, “AIDA: an adaptive image deconvolution algorithm with application to multi-frame and three-dimensional data,” J. Opt. Soc. Am. A **24**, 1580–1600 (2007). [CrossRef]

**18. **T. Kenig, Z. Kam, and A. Feuer, “Blind image deconvolution using machine learning for three-dimensional microscopy,” IEEE Trans. Pattern Anal. Mach. Intell. **32**, 2191–2204 (2010). [CrossRef] [PubMed]

**19. **M. Keuper, T. Schmidt, M. Temerinac-Ott, J. Padeken, P. Heun, O. Ronneberger, and T. Brox, “Blind deconvolution of widefield fluorescence microscopic data by regularization of the optical transfer function (OTF),” in *Proceedings of IEEE Conference on Computer Vision and Pattern Recognition* (IEEE, 2013), pp. 2179–2186.

**20. **J. Markham and J.-A. Conchello, “Parametric blind deconvolution: a robust method for the simultaneous estimation of image and blur,” J. Opt. Soc. Am. A **16**, 2377–2391 (1999). [CrossRef]

**21. **P. Pankajakshan, B. Zhang, L. Blanc-Féraud, Z. Kam, J.-C. Olivo-Marin, and J. Zerubia, “Blind deconvolution for thin-layered confocal imaging,” Appl. Opt. **48**, 4437–4448 (2009). [CrossRef] [PubMed]

**22. **F. Soulez, L. Denis, Y. Tourneur, and É. Thiébaut, “Blind deconvolution of 3D data in wide field fluorescence microscopy,” in *Proceedings of IEEE International Symposium on Biomedical Imaging* (IEEE, 2012), pp. 1735–1738.

**23. **F. Xue and T. Blu, “A novel SURE-based criterion for parametric PSF estimation,” IEEE Trans. Image Process. **24**, 595–607 (2015). [CrossRef]

**24. **J. Li, F. Xue, and T. Blu, “Gaussian blur estimation for photon-limited images,” in *Proceedings of IEEE International Conference on Image Processing* (IEEE, 2017), pp. 495–499.

**25. **F. Luisier, T. Blu, and M. Unser, “Image denoising in mixed Poisson-Gaussian noise,” IEEE Trans. Image Process. **20**, 696–708 (2011). [CrossRef]

**26. **Y. Hiraoka, J. W. Sedat, and D. A. Agard, “Determination of three-dimensional imaging properties of a light microscope system. partial confocal behavior in epifluorescence microscopy,” Biophys. J. **57**, 325–333 (1990). [CrossRef] [PubMed]

**27. **C. Preza and J.-A. Conchello, “Depth-variant maximum-likelihood restoration for three-dimensional fluorescence microscopy,” J. Opt. Soc. Am. A **21**, 1593–1601 (2004). [CrossRef]

**28. **J. Kim, S. An, S. Ahn, and B. Kim, “Depth-variant deconvolution of 3D widefield fluorescence microscopy using the penalized maximum likelihood estimation method,” Opt. Express **21**, 27668–27681 (2013). [CrossRef]

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

**30. **H. Kirshner, F. Aguet, D. Sage, and M. Unser, “3-D PSF fitting for fluorescence microscopy: implementation and localization application,” J. Microsc. **249**, 13–25 (2013). [CrossRef]

**31. **S. Liu, E. B. Kromann, W. D. Krueger, J. Bewersdorf, and K. A. Lidke, “Three dimensional single molecule localization using a phase retrieved pupil function,” Opt. Express **21**, 29462–29487 (2013). [CrossRef]

**32. **J. Huang, M. Sun, K. Gumpper, Y. Chi, and J. Ma, “3D multifocus astigmatism and compressed sensing (3D MACS) based superresolution reconstruction,” Biomed. Opt. Express **6**, 902–917 (2015). [CrossRef] [PubMed]

**33. **M. Štefko, B. Ottino, K. M. Douglass, and S. Manley, “Design principles for autonomous illumination control in localization microscopy,” bioRxiv 295519; doi:https://doi.org/10.1101/295519.

**34. **J. Li, F. Xue, and T. Blu, “Fast and accurate three-dimensional point spread function computation for fluorescence microscopy,” J. Opt. Soc. Am. A **34**, 1029–1034 (2017). [CrossRef]

**35. **C. Smith, M. Huisman, M. Siemons, D. Grünwald, and S. Stallinga, “Simultaneous measurement of emission color and 3d position of single molecules,” Opt. Express **24**, 4996–5013 (2016). [CrossRef] [PubMed]

**36. **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]

**37. **F. Aguet, “Super-resolution fluorescence microscopy based on physical models,” Ph.D. thesis, Ecole Polytechnique Fédérale de Lausanne (2009).

**38. **M. Arigovindan, J. Shaevitz, J. McGowan, J. W. Sedat, and D. A. Agard, “A parallel product-convolution approach for representing depth varying point spread functions in 3D widefield microscopy based on principal component analysis,” Opt. Express **18**, 6461–6476 (2010). [CrossRef] [PubMed]

**39. **B. Kim and T. Naemura, “Blind depth-variant deconvolution of 3D data in wide-field fluorescence microscopy,” Sci. Rep. **5**, 9894 (2015). [CrossRef] [PubMed]

**40. **N. Patwary and C. Preza, “Image restoration for three-dimensional fluorescence microscopy using an orthonormal basis for efficient representation of depth-variant point-spread functions,” Biomed. Opt. Express **6**, 3826–3841 (2015). [CrossRef] [PubMed]

**41. **C. Vonesch, F. Aguet, J.-L. Vonesch, and M. Unser, “The colored revolution of bioimaging,” IEEE Signal Process. Mag. **23**, 20–31 (2006). [CrossRef]

**42. **J. Li, F. Luisier, and T. Blu, “PURE-LET Image Deconvolution,” IEEE Trans. Image Process. **27**, 92–105 (2018). [CrossRef]

**43. **J. Boulanger, C. Kervrann, P. Bouthemy, P. Elbau, J.-B. Sibarita, and J. Salamero, “Patch-based nonlocal functional for denoising fluorescence microscopy image sequences,” IEEE Trans. Med. Imag. **29**, 442–454 (2010). [CrossRef]

**44. **J. Li, F. Xue, and T. Blu, “Accurate 3D PSF estimation from a wide-field microscopy image,” in *Proceedings of IEEE International Symposium on Biomedical Imaging*, (IEEE, 2018), pp. 501–504.

**45. **B. Zhang, J. Zerubia, and J.-C. Olivo-Marin, “Gaussian approximations of fluorescence microscope point-spread function models,” Appl. Opt. **46**, 1819–1829 (2007). [CrossRef] [PubMed]

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

**47. **S. Dmitrieff and F. Nédélec, “ConfocalGN: A minimalistic confocal image generator,” SoftwareX **6**, 243–247 (2017). [CrossRef]

**48. **D. A. Fish, A. M. Brinicombe, E. R. Pike, J. G. Walker, and R. L. Algorithm, “Blind deconvolution by means of the Richardson-Lucy algorithm,” Appl. Opt. **12**, 58–65 (1995).

**49. **J.-S. Lee, T.-L. E. Wee, and C. M. Brown, “Calibration of wide-field deconvolution microscopy for quantitative fluorescence imaging,” J. Biomol. Tech. **25**, 31 (2014). [CrossRef] [PubMed]

**50. **M. Siemons, C. Hulleman, R. Thorsen, C. Smith, and S. Stallinga, “High precision wavefront control in point spread function engineering for single emitter localization,” Opt. Express **26**, 8397–8416 (2018). [CrossRef] [PubMed]