On-the-fly estimation of a microscopy point spread function

: A proper estimation of realistic point-spread function (PSF) in optical microscopy can signiﬁcantly improve the deconvolution performance and assist the microscope calibration process. In this work, by exemplifying 3D wide-ﬁeld ﬂuorescence 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 eﬀectiveness 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 suﬃciently ﬂexible to be generalized non-spherical aberrations and other microscope modalities.

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][13][14][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][18][19]. Alternatively, based on an analytic model and a carefully designed optimization criterion, one can estimate the PSF by finding its parameters first [20][21][22][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.

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][30][31][32][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, q) are akin to refractive indices, and their values are fixed once for all (see below). The parameter η is an indication of the focus position z p . More specifically, the PSF h is written as where r = (x, y, z) is the 3D spatial coordinate, 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 α β 2 − ρ 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].

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][39][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).
where Y ∈ R N denotes the distorted observation of the unknown true original image X ∈ R N , N = N x × N y × N z 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 α represents the gain of the acquisition device, which controls the noise during the acquisition process. H 0 ∈ R N ×N is a block-circulant matrix, which implements a discrete convolution with the PSF h. I denotes the identity matrix, and P(·) and N(·, ·) 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][23][24].  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]

Blur-PURE as an optimization criterion
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 H is Tikhonov regularized deconvolution: , 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 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- 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 ( ω 2 = ω 2 x + ω 2 y + γ ω 2 z ) 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.

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
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] where E{·} 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. E X T H T 0 UY = E Y T UY − σ 2 Tr(U) . 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.

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 H 0 ,µ by a linear combination  of basis functions U app = n a n U H n ,µ . Here H n 's are constructed by the corresponding PSF basis h n (r; η), 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: min H HW H − U app 2 2 , where H assumes the parametric expression in Eq. (2). We propose an iterative optimization algorithm: where and H (0) is randomly initialized. The solution H (k) involves finding the coefficients c (k) m '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-basis model validation
The PSF approximation error is measured in terms of the relative squared error (RSE) [34,45] calculated as where the ground-truth PSF 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 nm, and 1025 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 p and the estimated η. The mean and max difference between them are 25 nm and 46 nm respectively. As we can see from Fig. 3(b) and other experiments (not shown here), the relation between η and z p is deterministic (function variables may depend on wavelength λ, 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 p , which may have some interest beyond the scope of this paper.

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 i = 1.515). The focus position z p varies from −1.24 µ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 s = 1.394 and λ = 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.

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 nm, and the measurement of the 3D stack is performed with a step size of 89 nm, 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.

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 i = 1.33, n s = 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 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.

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.  [5] with the measured PSF, respectively; (d) restored image by the blind-RL algorithm; (e) restored image by the EpiDEMIC plugin; (f) restored image by the PURE-LET algorithm with the PSF estimated by our approach (described in Section 3). Locators (yellow line) indicate the location of displayed sections (z = 0). Zoomed-in regions (×2) show the resolution improvement of our method in (f). Also note the significant reduction of out-of-focus blur in the x-z plane. Scale bar: 1µm.

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. Various source codes are available at http://www.ee.cuhk.edu.hk/~tblu/demos.
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].