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:

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 α′ − β′ρ². 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).

Based on the image acquisition model in Eq. (3), the objective of this work is to estimate the PSF h (namely H₀) 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₀ can be done by minimizing the expected mean squared error (MSE) between the blurred ground-truth image H₀X and a linear processing U_H of the acquired image Y [23,24]: $\frac{1}{N}ℰ\left\{{\Vert {\mathbf{U}}_{\mathbf{H}}\text{Y}-{\mathbf{H}}_0\text{X}\Vert }^2\right\}$. This oracle criterion (H₀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: U_H,μ = HW_H, where W_H = H^T(H^TH + μP)⁻¹, 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₀, through ${\mathbf{H}}^{\text{T}}\mathbf{H}={\mathbf{H}}_0^{\text{T}}{\mathbf{H}}_0$, i.e., ‖ĥ(ω)‖² = ‖ĥ₀(ω)‖² in the Fourier domain. Since ĥ and ĥ₀ satisfy a very low order parametric representation, this equation is equivalent to h(r) = h₀(r − r₀) where r₀ 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₀(ω) = |H₀(ω)|² (|H₀(ω)|² + C(ω)/A(ω))⁻¹, 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(ω) ≈ μ‖ω‖² 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.,

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]) $ℰ\left\{{\text{X}}^{\text{T}}{\mathbf{H}}_0^{\text{T}}\mathbf{U}\text{Y}\right\}=ℰ\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. $ℰ\left\{{\text{X}}^{\text{T}}{\mathbf{H}}_0^{\text{T}}\mathbf{U}\text{Y}\right\}=ℰ\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⁶) 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⁻⁴. 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.

 

Fig. 1 (a) Reliability of blur-PURE as an estimate of the blur-MSE: a typical example of the Bars image degraded by a PSF and various noise parameters (α = 0.02, σ² ∈ [0, 25]). The values of the blur-PURE for various values of σ closely match those of the blur-MSE, which confirms that the blur-PURE is an accurate substitute to the blur-MSE in practice. 3D rendering of the blurred noisy image when σ² = 16 is shown at the right bottom. See details of the experimental setting in Section 5.1. (b) and (c) show the z and x profiles of band-indicators (true frequency responses of U_H0,μ and approximated U_app), respectively.

Download Full Size | PDF

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_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⁻³(αY_mean + σ²), 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:

 

Fig. 2 The scheme diagram of the proposed PSF estimation approach. With the help of the 4-basis PSF approximation model, the band-indicator U_app is firstly estimated by minimizing the blur-PURE in Theorem 1, and then the estimated PSF H_Est is obtained by an iterative optimization algorithm in Eq. (4).

Download Full Size | PDF

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

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.

 

Fig. 3 (a) Determination of the optimal η in the sense of minimal RSE: typical examples of three different focus positions z_p = 0, 500 nm and 1000 nm. Other settings: λ = 395 nm, NA = 1.4, n_s = 1.38, n_i = 1.515. The estimated value is η = 0, 529.10 nm and 1025.07 nm, respectively. (b) The scatter plot between the true z_p and the estimated η. The ideal relation is indicated by solid red line.

Download Full Size | PDF

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.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_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.

 

Fig. 4 Evaluation of the proposed approximation model on simulated PSFs. (a) Plot of the mean RSE between the theoretical and fitted PSFs over different wavelengths λ. (b) Axial intensity profile of typical PSFs (λ = 395nm, z_p = 500 nm) generated by the complete Gibson-Lanni model (the black curve) and fitting (the dashed red curve). The focal shift can be observed. The fitted PSF closely matches the ground truth PSF (RSE = 0.48%). (c) x–z section of the theoretical PSF. Note the non-symmetric pattern caused by the refractive index mismatch. The degraded PSFs in (d) and (f) are fitted by the proposed approximation model in Eq. (2) resulting in (e) and (g), respectively. Images have been cropped and rescaled for visualization purpose.

Download Full Size | PDF

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.

 

Fig. 5 A practical example. The x–y (a) and y–z (b) sections of the experimental PSF. (c–d) The x–y and y–z sections of the fitted PSF by the proposed approximation model in Eq. (2). The focus indicator (parameter η) is estimated to be 1.855 μm. The mismatch between the refractive indices of the sample and immersion medium contributes to an axially asymmetric PSF. Scale bar: 0.5 μm.

Download Full Size | PDF

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_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.

Table 1. RSE (%) comparison of the PSF estimation accuracy with other approaches under different scenarios. The results have been averaged over 5 random initializations.

View Table

 

Fig. 6 Comparison results of PSF estimations in the high noise conditions (α = 0.2, σ = 0.02). (a) z_p = 0 μm; (b) z_p = 2 μm. Other settings: λ = 622 nm, NA = 1.4, n_i = 1.33 and n_s = 1.46. Locators (orange line) indicate the location of displayed sections (z = 0). Images have been cropped and rescaled for visualization purpose.

Download Full Size | PDF

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%).

 

Fig. 7 The PSF estimation results of experimental data from Qdots. (a–b) the x–y and z–y sections of the estimated PSF. (c) the z profile of the PSF estimated with our algorithm (described in Section 3), compared to the z profiles of the experimental and fitted PSFs (same conditions as in Fig. 5).

Download Full Size | PDF

 

Fig. 8 Deconvolution results of experimental Qdots data. (a) raw image (256 × 256 × 113 voxels); (b) and (c) restored images by the Richardson-Lucy algorithm and PURE-LET algorithm [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.

Download Full Size | PDF

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