## Abstract

Conventional deconvolution methods assume that the microscopy system is spatially invariant, introducing considerable errors. We developed a method to more precisely estimate space-variant point-spread functions from sparse measurements. To this end, a space-variant version of deblurring algorithm was developed and combined with a total-variation regularization. Validation with both simulation and real data showed that our PSF model is more accurate than the piecewise-invariant model and the blending model. Comparing with the orthogonal basis decomposition based PSF model, our proposed model also performed with a considerable improvement. We also evaluated the proposed deblurring algorithm. Our new deblurring algorithm showed a significantly better signal-to-noise ratio and higher image quality than those of the conventional space-invariant algorithm.

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

## 1. Introduction

3-D deconvolution is a powerful computational method for abrogating the effects of out-of-focus light originating from different depth planes in the specimen for any microscope system [1], such as wide-field microscopy (WFM) [2–4], laser scanning confocal microscopy (LSCM) [5–7] and light sheet fluorescence microscopy (LSFM) [8, 9]. This technique has been proven to be highly effective in removing blur and noise introduced during the imaging process, thus improving both the spatial resolution and the signal-to-noise ratio [10, 11]. As a linear system, the fluorescence microscope can be characterized by its 3-D point-spread function (PSF), which plays an important role in the 3-D deconvolution.

Conventional deconvolution methods accept the assumption of a space-invariant PSF (SI-PSF) [6, 12]. Most published deconvolution softwares, such as the commercial software *Huygens* and the open-source one *DeconvolutionLab2* [13], also assume the system PSF invariant. However, microscope systems are not spatially invariant [14, 15]. The assumption of SI-PSF limits the algorithm’s performance and introduces significant errors into the results. Many factors may lead to spatial variation in the PSF. Generally, due to the differences in the optical path and the refractive index mismatch between the medium and sample, the PSF varies significantly along the optical axis [3, 16]. In a thick specimen, this type of axial variation caused by the depth change of object point is the dominant component of the spatial variation. In addition, non-uniformity of illumination may also generate spatial variation in the PSF. For example, in a LSFM system, variation in the thickness of the light sheet, which is either a Gaussian or Bessel [8,9,15] beam, leads to radial variation in the PSF. As a result, the LSFM system achieves satisfactory resolution only in the central area of the light sheet, where its PSFs are most compact. Furthermore, stage-drift may also cause spatial variation in the PSF.

To improve the 3-D deconvolution performance, recent research studies have considered variation of the PSF [14, 17]. In [18, 19], the wavefront encoding were used to convert the SV system to a space invariant one so that a single PSF can be used in the restoration process. Because of the difficulty of measuring PSF, some studies [20–22] estimated the restored image and the system PSFs simultaneously. Although various theoretical models of the space-variant PSF (SV-PSF), such as the Gibson and Lanni model [23] and Richards and Wolf model [24], have been proposed over the past several decades, PSFs generated by these models are still inaccurate for practical systems due to the inaccuracy of the model itself and the inevitable errors of the critical parameters in the models. Therefore, the SV-PSF is usually determined using a series of sparse measurements. A simple method that uses a piecewise-invariant PSF to approximate the SV-PSF was proposed in previous studies [17], with PSF treated as invariant within a sub-region. In this model, the error of the estimated PSF increases as its distance to the measured PSF grows larger. The piecewise-invariant PSF model introduces considerable artifacts in contiguous areas of the sub-regions. Other methods involve numerical interpolation, such as the blending PSF model. For example, in some studies [25], the PSF was determined by two neighbor measurements using a linear blending function, whereas in another study [15], PSF at specific positions was interpolating using Nadaraya-Watson kernel regression with a Gaussian kernel. However, the blending PSF model only provides slightly better accuracy than that of the piecewise-invariant model when the PSF measurements are sparse. A method that decomposed the PSF to a series of orthonormal basis was also presented in study [4, 26, 27]. And in [14], linear combination and Zernike interpolation were combined for SV-PSF estimation. Recently, PSF fitting methods were introduced to this problem. In these methods, the 3-D voxel data of measurements were used to fit a theoretical model such as the Gaussian model [28] or the Gibson and Lanni model [29]. Nevertheless, fitting with a theoretical model has the following limitations: 1) The theoretical model itself may be inaccurate. For example, the Gaussian model is acceptable for LSCM but may introduce a substantial error into WFM or LSFM systems. 2) Sophisticated models such as the Gibson and Lanni model contain too many parameters. In a fitting problem, the number of parameters is restricted to avoid overfitting or reaching a local optimum. In this study, we aim to establish a more accurate method for SV-PSF estimation by considering the formation of spatial variation to the PSF, which is the fundamental problem of establishing a SV deblurring method for WFM, LSFM or other microscopies.

State-of-the-art deconvolution algorithms for 3-D fluorescent images assume that the PSF is shift variant. Algorithms based on the possibility model have been developed, such as the maximum-likelihood estimation method (MLE) [25,30] and maximum a posteriori algorithm [7,31,32]. Due to the ability of inferring out of band frequencies and the more accurate assumption of noise model, MLE has the superior capability to remove noise and is widely used in fluorescence image restoration. However, to the best of our knowledge, none of these methods proposed an effective solution to estimate the PSF with a matching optical path length (OPL). The variation of OPL from the observed object point to the back focal plane is known to be a dominant cause of PSF variation. Therefore, PSF estimation should proceed according to its OPL, which is difficult to acquire in measurements. Conventional methods directly scale the axial distance with the paraxial formula according to the reading of the stage position [33,34], which may reduce the accuracy of the PSFs.

In addition, there is a risk that the optimization process is ill-posed and that iteration may produce a noisy image. Therefore, a regularization term based on prior information must be incorporated to improve the algorithm performance. For example, the energy-based regularization presented by Conchello et al. leaves oscillations in the homogeneous area [35,36]. The Tikhonov-Miller regularization results in smoothed edges [37]. A total variation (TV) regularization was proposed to prevent noise amplification and simultaneously preserve the edges of the images [38,39].

In the current study, we aim to develop a practical 3-D space-variant (SV) image deblurring method with a measuring model of SV-PSF. We will first describe the formation of PSF variation and the main stages of our restoration framework. We then introduce the derivation of our proposed SV-PSF model with the similarity transformation and estimate the PSF located at a certain depth with this PSF measuring model. Additionally, we consider the stochastic error introduced by stage drift to achieve higher precision. In addition, we propose a space-variant maximum-likelihood (ML) method with TV regularization and present one solution to determine the actual position of the image region where the PSF should be densely estimated. Finally, we demonstrate that our proposed PSF model and restoration method provide higher resolution and more accurate restoration in experiments.

## 2. General framework of space-variant deblurring

In a SV system, the PSF varies due to the change in the position of the point source. In general, spatial variation can be categorized into radial and axial variation, which have different causes. In a microscope system, the response of a point source is a diffraction pattern image produced by the limiting aperture and spherical aberration. From the object points located at different depths in the specimen, light must pass through multiple layers of the sample, coverslip, and immersion oil, each of which may have significantly different refractive indices and thicknesses, resulting in an optical path difference (OPD, Fig. 1). Consequently, the diffusion pattern varies along the optical axis, and the OPD leads to an approximately periodic change in the concentric ring of the airy disk. In general, axial variation in the PSF that is caused by differences in depth is the most significant component in a microscope system. Therefore, a proper model capable of expressing the axial variation of the PSF is important for modeling and estimating the SV-PSF. Furthermore, the effects of stage drift may also produce stochastic axial variation in the PSF. Stage drift refers to random bias from the designed path due to the mechanical tolerance of the stage guide rail. See section 3.2 for the effects of stage drift on the PSF measurements. Nevertheless, the drift path commonly has poor reproducibility and differs among different imaging procedures. Stage drift should be estimated and corrected for both the PSF measurement and the sample image to effectively improve the accuracy of deconvolution. Radial variation in the PSF is commonly caused by objective-lens aberrations such as curvature of the field and coma aberration, which is correlated with the change in the radial distance between the object point and the optical axis. In the high-quality flat-field lenses of a commercial microscope system, these off-axis aberrations are negligible under most conditions. However, in certain microscopy techniques wherein illumination varies within the focal plane, non-uniform illumination may generate variation of the PSF in the radial plane. For example, the PSF varies within the light sheet (either Gaussian or Bessel beam) in a LSFM system [8, 15]. The thickness change of the light sheet and the non-uniform intensity distribution lead to planar variation in the PSF.

To improve the deblurring performance for imaging thick specimens, we proposed a practical framework of 3-D SV deconvolution (Fig. 2) together with a novel method for modeling the SV-PSF. The proposed framework was divided into two stages: dense estimation for the SV-PSF and SV restoration. In the first stage, dense SV-PSFs were estimated according to the sparse measurements of randomly distributed fluorescent beads. To accurately estimate the SV-PSF, the formation of spatial variation must be considered into the model. In the second stage, we developed a 3-D SV restoration algorithm optimized using the ML method and based on the dense estimation of the SV-PSF. To avoid noise amplification and to preserve object edges, a proper regularization term should be designed. For a SV deblurring problem, it is important that the optical paths of the estimated SV-PSFs and that of the object points in the sample image are estimated and matched. In an optimal SV system, the object point should be deconvolved only by the corresponding PSF such that the optical path is exactly the same as that of the object point. Nonetheless, the optical path cannot be simply calculated according to the z-position readings of the stage during sample imaging because the reading reflects only the position of the stage itself but not the position of the object point. Therefore, estimating and matching the optical paths of the SV-PSFs and the sample image are critical for developing an improved deconvolution method.

## 3. Estimating the SV-PSF from sparse measurements

To measure the SV-PSF, images of sparsely distributed fluorophore beads of known diameter were recorded. We assumed that the PSF of the system smoothly varies over the image both axially and radially. To precisely estimate the SV-PSF, the axial and radial variations were correspondingly considered in the model. To extract the PSFs, the bead images were smoothed using a Gaussian filter, and local maxima were computed as coarse bead locations. Image of a single bead was then cropped from the original image. Because the fluorophore bead is a sphere with a known diameter, the PSFs could be computed by deconvolving the measured images with the prior knowledge of the bead shape [15]. The intensity distribution of the sphere could be modeled as a Gaussian function, wherein the standard deviation is chosen according to its diameter and the estimated coarse center. We then obtained a PSF measurements set $\tilde{H}:=\left\{\tilde{{h}_{i}}(\mathbf{s})|i=1,\dots ,N\right\}$, where $\tilde{{h}_{i}}(\mathbf{s})$ refers to the discrete matrix representing one of the measured PSF and the coordinate **s** = (*u, v, w*)* ^{T}*, and a position set $\tilde{P}:=\left\{\tilde{{p}_{i}}=\left(\tilde{{x}_{i}},\tilde{{y}_{i}},\tilde{{z}_{i}}\right)|i=1,\dots ,N\right\}$ where ${\tilde{p}}_{n}$ refers to the fine center position of $\tilde{{h}_{i}}(\mathbf{s})$ calculating by the 3-D radial symmetry localization algorithm in the study [40]. The obtained PSF measurements were interpolated to form a center-adjusted PSF set

*H*:= {

*h*(

_{i}**s**) |

*i*= 1,…,

*N*} and its corresponding position set is

*P*:= {

*p*= (

_{i}*x*|

_{i}, y_{i}, z_{i}*i*= 1, …,

*N*}, where

*p*refer to the center position of the center-adjusted PSF

_{i}*h*(

_{i}**s**), to prepare for the PSF modeling. See Fig. 3 for the diagram of the PSF estimation process.

#### 3.1. Modeling and estimating the axial variation of PSF

In general, the axial shift variation in the PSF caused by depth change is the dominant component of the spatial variation in a microscope system. In this work, we considered establishing a method to generate a dense estimation of the SV-PSF from sparse measurements considering the depth variation in the PSF. Assuming *h _{i}* (

**s**) and

*h*(

_{j}**s**) are two PSF measurements from set

*H*, we assume that there exists a function $\mathcal{T}(\cdot )$ that allows a measured PSF

*h*(

_{i}**s**) to express the PSF

*h*(

_{j}**s**) at any specific position.

*h*(

_{i}**s**) refers to the intensity of light at

**s**= (

*u, v, w*)

*of the measured PSF.*

^{T}For object points located at different depths in the specimen, light from those points must pass through multiple layers with significantly varying thicknesses and refractive indices, resulting in OPD. The image of a point-source object is known to be the diffraction pattern produced when a spherical wave converges with the image point on the detector plane. The OPD between two object points leads to a phase shift in their wavefront, resulting in an approximately periodic change in the concentric ring of the airy disk. Therefore, variation of the diffraction pattern can be well approximated using a similarity transformation process. Function $\mathcal{T}(\cdot )$ that allows *h _{i}* (

**s**) to approximate

*h*(

_{j}**s**) was defined by a similarity transformation,

*θ*is the similarity transformation matrix, where (

_{ij}*α*,

_{ij}*β*) refer to the scaling parameters which transform

_{ij}*h*(

_{i}**s**) to approximate

*h*(

_{j}**s**) in the radial and axial direction, respectively. Then, the residual error is

*·*‖ refers to the L2-norm and the parameter

*θ*can be computed as

_{ij}*θ*and $\frac{\partial {h}_{i}\left({\theta}_{ij}\mathbf{s}\right)}{\partial \left(\alpha u\right)}$, $\frac{\partial {h}_{i}\left({\theta}_{ij}\mathbf{s}\right)}{\partial \left(\alpha v\right)}$, $\frac{\partial {h}_{i}\left({\theta}_{ij}\mathbf{s}\right)}{\partial \left(\beta w\right)}$ refer to the gradient of

_{ij}*h*(

_{i}*θ*

_{ij}**s**) along the x-axis, y-axis and z-axis, respectively.

We evaluated our PSF approximation model, which is based on the similarity transformation, for both synthetic and real data. The synthetic PSFs were generated using an ImageJ plugin called the PSF Generator [29] (http://bigwww.epfl.ch/algorithms/psfgenerator/) with Gibson and Lanni model. The real data was recorded from sparsely distributed fluorophore beads. The similarity transform *θ _{ij}* between pairwise PSFs was calculated using equation (5) and Table 1. The results from both synthetic and real data clearly implied that the parameters of similarity transform between PSF pairs are linear functions of their depth variation within a small range, which is

*α*=

*k*Δ

_{r}*z*,

*β*=

*k*Δ

_{z}*z*and Δ

*z*refers to the relative depth (Figs. 4(a) and 4(d)). Fig. 4(g) shows that

*k*and

_{r}*k*decrease when beads go deeper. We also compared the approximate difference between PSF pairs with and without similarity transformation (Figs. 4(b) and 4(e)). The results indicated that our model could significantly reduce the approximation error. Therefore, we established a method to estimate the SV-PSF from sparse measurements based on the similarity transformation.

_{z}To reduce error, the proposed method uses all of the PSF measurements efficiently. The PSF estimation $\widehat{h}(\mathbf{s};{\mathbf{s}}_{\mathbf{0}})$ whose point source is located at the coordinate **s _{0}** = (

*u*

_{0},

*v*

_{0},

*w*

_{0})

*in the image domain, is defined as a weighted average of all its approximations $\mathcal{T}({h}_{i}(\mathbf{s});{\theta}_{i})$ from measurements using the similarity transform.*

^{T}*θ*refers to the similarity transformation matrix that allows

_{i}*h*(

_{i}**s**) to approximates $\widehat{h}(\mathbf{s};{\mathbf{s}}_{\mathbf{0}})$,

*A*is the constraint of normalization, and the weighting function

*z*refers to the axial position of the object point

**s**,

_{0}*z*refers to the axial center position of

_{i}*h*(

_{i}**s**), and Δ

*z*refers to the threshold of the relative depth Δ

_{thresh}*z*.

To evaluate the accuracy of the proposed SV-PSF estimation method, we compared our method with the piecewise-invariant model [17] and the blending model [25] on both the synthetic and real measured PSF data. In the synthetic data tests, the MSE between the estimated PSF and the ground truth were compared (Fig. 4(c)). In the real data tests, some of the data were used as measurements to estimate the PSF and to compare with the remaining portion, which was treated as ground truth (Fig. 4(f)). The results indicated that the PSF from our model is more accurate than that from the piecewise-invariant and the blending PSF models for both synthetic and real data. In addition, as for the real data, our model remarkably reduced the MSE in the region of sparse measurements, while in the region where the PSFs were densely measured our method is still slightly better. See Fig. 5 for one example of PSF estimation with the piecewise-invariant model, the blending model and our proposed model. Besides, our PSF model was compared with PCA [4, 26] under different conditions (Figs. 4(h) and 4(i)) and the results showed that our proposed model was more precise in sparse measurements.

#### 3.2. Stage drift estimation and correction

Stage drift may result in both axial variation in the PSF and a stochastic lateral shift in the image layer. As above discussed, the drift path has poor reproducibility (Figs. 6(a) and 6(c)). Different distortions in the main axis between the sample imaging and PSF measurement cause significant degradation of the restoration results. Therefore, stage drift should be corrected in both the sample image and the PSF measurement by aligning the image layers.

Stage drift correction was performed by estimating the drift paths and aligning them to a vertical line. The stage drift of the PSF measurements was estimated by computing the symmetric centers of the PSF in each layer. Drift path estimation for a sample image may use inserted fluorophore beads or small structures on the sample. For one PSF measuring image, which contained multiple beads, the path of the drift was calculated by averaging the distortion estimated from all beads. Therefore, the PSFs in the dataset *H*:= {*h _{i}* (

**s**) |

*i*= 1,…,

*N*} were corrected by numerical interpolation for the voxel data of the PSF measurements according to the drift path of the stage. Nevertheless, the localization of the PSF center was also fine-tuned. A 3-D radial-symmetry-based localization method [40] was employed to calculate the center of the PSF at sub-pixel resolution. Figs. 6(b) and 6(d) show the drifted and corrected PSF.

#### 3.3. The radial variation of the PSF in LSFM

In WFM and LSCM systems, the radial shift variation is mainly caused by objective-lens aberrations such as curvature of field, coma aberration, and astigmatism. In high-quality flat-field lenses, these off-axis aberrations are sufficiently minimized to be negligible under most conditions over the central image field. However, in LSFM, samples are excited by Gaussian or Bessel beams, and such settings can lead to blurring with a PSF that varies within the non-uniform beam-excited LSFM [8, 15]. Considering the non-uniformity of the excitation light sheet, a system PSF of the LSFM is computed as the pointwise product of the intrinsic PSF and the intensity distribution of its illumination field [43] (Fig. 7). Here, the system between the object point and detector is characterized as the intrinsic PSF. Assuming that the illumination beam is uniform along the y-axis, the system PSF can be expressed as follows

where*SH*(

**s**;

*x*

_{0},

*z*

_{0}) is the system PSF located at

*x*=

*x*

_{0},

*z*=

*z*

_{0},

*IH*(

**s**;

*z*

_{0}) refers to the intrinsic PSF at depth

*z*=

*z*

_{0}, and

*ID*(

**s**;

*x*

_{0}) refers to the intensity distribution of its excitation beam located at the cross-section

*x*=

*x*

_{0}.

## 4. 3-D space-variant restoration

Comparing with the axial variation, the radial variation of PSF is minimized to be negligible in commercial microscope systems after stage-drift correction. Therefore, we simplified the SV system to a depth-variant (DV) one. Bello noted that the linear time-variant (LTV) system can be described via an integral transform with a time-dependent kernel [44]. Similarly, as the microscope is linearly DV in practice, we defined an operator ⊗ as the discrete version of the DV quasi-convolution operator in the image domain.

where*M*is the number of strata and

*f*(

_{m}**s**) contains the

*m*th stratum of the object.

*h*(

_{m}**s**) is the system response to a point source placed at the

*m*th stratum. As

*h*(

_{m}**s**) can be approximated by a linear combination of several transforming PSFs from the PSF measurements set, the operator ⊗ can be simplified as

*f*(

_{m}**s**) and

*h*(

_{i}**s**) only need to computed once that avoiding considerable repeated computation. Based on the depth-variant ML restoration method proposed in [25], it is easy to derive the iteration of DV deconvolution as follows

*m*th stratum of

*f*

^{(}

^{k}^{)}(

**s**), and

*f*

^{(}

^{k}^{)}(

**s**) is the estimate of the raw image after

*k*iterations.

For a DV system, the PSF is a function of its actual depth. In other words, the object point should be convolved only by the PSF for that depth matched. Unfortunately, it is difficult, even impossible to determine the axial position *z* according to readings from the stage or image coordinate *m* during sample imaging, even if the same thickness of the coverslip and the refractive index of the immersion oil and coverslip were given. Because the axial position reading could guarantee only the position of the stage itself rather than the position of the object point. Axial variation in the PSF is a function of the OPL from the observed object point to the back focal plane. Therefore, estimating and matching the OPL of the PSF measurement and sample imaging is an essential problem.

As the axial image position *z* is unknown, we proposed an alternating optimization solution that combined a depth-updating iterative procedure into each image-updating iterative step. The axial location *z* was updated using the Newton’s method by minimizing the log-likelihood function

*h*(

**s**;

**s**) in

_{0}*J*(

*f*) was related to its axial location

*z*, we rewrote

*J*(

*f*) as

*J*(

*f, z*). In the first image-updating iteration, the axial location

*z*was set according to the scaling axial distance with the reading of the stage depth [33, 34]. In later iterations, the initial z used the last estimation value of depth. The depth-updating iteration is expressed as

*J*(

*f, z*) with respect to depth

*z*is shown as

*h*(

**s**;

**s**) with respect to depth

_{0}*z*is written as

*k*and

_{r}*k*refer to the computed slope of the scaling parameters along the radial and axial directions. In this manner, the axial location of sample

_{z}*z*converges to an optimum to match the optical paths of the PSF and image sample. Thereafter, we can achieve a more precise image restoration than that of the conventional space-invariant method.

To solve the above ill-posed inverse problem, a regularization term must be carefully selected. The restoration result is always a trade-off between preserving details and avoiding noise amplification. The TV regularization has been proven to be highly effective in preserving the edges in the image and smoothing in homogeneous areas [38]. Therefore, we introduced a TV regularization term to the proposed SV deconvolution algorithm and the regularized likelihood function is shown as

*J*(

_{TV}*f*) refers to the TV regularization term,

*λ*is the regularization parameter, |·| refers to the L1-norm and ∇ refers to the gradient. The iteration of the proposed DV deconvolution algorithm is modified as

_{TV}## 5. Evaluation

We evaluated our proposed SV deconvolution algorithm and compared it with the conventional deconvolution algorithm on synthetic data and real data. The mean square error (MSE), peak signal-to-noise ratio (PSNR) and the normal mean I-divergence criteria were used to quantify the quality of the deblurred image. Because the ground truth is unknown in a real data experiment, real data were used to validate the algorithm without quantitative analysis. See Fig. 8 for the deblurring experiments on three synthetic test samples: three balls, a cavity and multi-objects. These synthetic object images were intentionally blurred with SV-PSFs and then degraded by Poisson noise. The restoration results from the conventional ML algorithm (SI-MLTV) and our proposed method (SV-MLTV) are compared and shown in Fig. 8 respectively.

The restored results using space-invariant algorithm show significant residual distortion related to stage drift and ring-shaped artifacts resulting from variation in the PSF. Additionally, the results restored by our proposed method show sharper edges and a more accurate distribution of intensity than the restored images using the constant PSF. For further validation, we introduced the proposed restoration method into more complicated images (Figs. 9 and 10). The cultured neuronal axons were acquired with a Leica SP8 microscope system and an Andor iXon EMCCD camera under both confocal and wide-field mode. To evaluate with the ground truth, confocal z-stack image was blurred with varying synthetic PSFs to simulate a wide-field data. Real Images was captured under the wide-filed mode was also tested. Our restoration results show sharper edges, better continuity of axonal fibers, and fewer artifacts comparing with the conventional SI methods (Figs. 9, 10 and 11). We also evaluated our proposed restoration method in a particle detection application scenario (Fig. 12). In our previous study, we reported an algorithm for detecting and tracking mitochondria, organelles critical for cellular survival and function [45]. The numerical results of particle detection in the images restored by different deconvolution algorithms is shown in Table 2. We used the F-value to measure the quality of particle detection. More particles were correctly detected from the restored image with our method, resulting much higher F-value than conventional method. We calculate the numerical criteria mentioned and our proposed method shows significant improvements in MSE, PSNR and I-divergence (Table 3).

## 6. Conclusions

In this work, we studied the 3-D SV deblurring problem for fluorescent image restoration and reconstruction. A SV-PSF model was proposed that enables generating a dense estimation of the SV-PSF from sparse measurements. In this model, the formation mechanism of axial and radial variation was considered. We proved that the similarity transformation between PSF pairs is approximately a linear function of their depth variation and is a much more accurate way to model the SV-PSF than the piecewise-invariant model and the blending model. Besides, the proposed PSF model was verified to be more precise than the PCA based model in sparse measurements. In our proposed SV-PSF estimation model, the PSF was defined as a weighted average of all of its approximations via similarity transformation. Thus, we could compute the PSF at any specific depth by effectively using all of the PSF measurements. With both synthetic and real data, our SV-PSF model was proven to be more accurate than conventional methods. We also theoretically discussed PSF estimation of the LSFM system. We defined an intrinsic PSF that characterized the system from the object point to the detector. Given the light beam intensity distribution and intrinsic PSF, we presented an efficient method for estimating the system PSF of a LSFM system. Furthermore, we proved that the stage drift in a commercial microscope system should not be ignored and presented a detailed drift correction method.

Similar to the problem in the LTV system, the SV quasi-convolution operator was defined, and the SV version of the deconvolution algorithm was given based on the depth-variant ML restoration method. In a linear SV system, the PSF is a function of its actual position, i.e., the object point should be convolved only by the PSF that optical path matched. The axial position could not be estimated accurately from the reading of the stage or the point coordinate of the sample image. The conventional method ignored this problem, introducing significant error into the restoration results. To solve this problem, we proposed an alternating optimization solution that combine a depth-updating iteration into each conventional image-updating iterative step. The depth-updating iteration computes the image location by maximizing the conditional probability of the given observed image using the Newton’s method. We evaluated our proposed SV deconvolution method with both synthetic and real data. Compared with the results of the conventional space-invariant method, our results showed a significant improvement in abrogating out-of-focus light and artifact noise suppression.

In this work, by solving the problems of SV-PSF modeling and establishing a SV restoration algorithm, we presented a practical method to significantly improve deblurring for fluorescent image. The algorithm was coded with MATLAB. We provided the source code and the sample data online at http://ese.nju.edu.cn/yogo/svdecon.zip for evaluation. This method should be highly useful for image restoration and 3-D reconstruction of images.

## Funding

Natural Science Foundation of Jiangsu, China (BK20161402); National Natural Science Foundation of China (91132710; 31671174; 31501133; 31671452); National Institutes of Health (NIH) (R01CA175360).

## References and links

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

**2. **J. Gallagher and N. Devaney, “Wavefront sensing and deconvolution of widefield images of biological samples for microscopy,” in “Frontiers in Optics 2012/Laser Science XXVIII,” (Optical Society of America, 2012), p. FTu3A.37.

**3. **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 (2013). [CrossRef]

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

**5. **C. Roider, R. Heintzmann, R. Piestun, and A. Jesacher, “Deconvolution approach for 3d scanning microscopy with helical phase engineering,” Opt. Express **24**, 15456–15467 (2016). [CrossRef] [PubMed]

**6. **M. Laasmaa, M. Vendelin, and P. Peterson, “Application of regularized richardson-lucy algorithm for deconvolution of confocal microscopy images,” J. Microsc. **243**, 124–140 (2011). [CrossRef] [PubMed]

**7. **G. Vicidomini, P. P. Mondal, and A. Diaspro, “Fuzzy logic and maximum a posteriori-based image restoration for confocal microscopy,” Opt. Lett. **31**, 3582–3584 (2006). [CrossRef] [PubMed]

**8. **Y. Wu, P. Wawrzusin, J. Senseney, R. S. Fischer, R. Christensen, A. Santella, A. G. York, P. W. Winter, C. M. Waterman, and Z. Bao, “Spatially isotropic four-dimensional imaging with dual-view plane illumination microscopy,” Nat. Biotechnol. **31**, 1032–1038 (2013). [CrossRef] [PubMed]

**9. **A. Kumar, Y. Wu, R. Christensen, P. Chandris, W. Gandler, E. Mccreedy, A. Bokinsky, D. A. Colonramos, Z. Bao, and M. Mcauliffe, “Dual-view plane illumination microscopy for rapid and spatially isotropic imaging,” Nat. Protoc. **9**, 2555 (2014). [CrossRef] [PubMed]

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

**11. **M. Castello, A. Diaspro, and G. Vicidomini, “Multi-images deconvolution improves signal-to-noise ratio on gated stimulated emission depletion microscopy,” Appl. Phys. Lett. **105**, 234106 (2014). [CrossRef]

**12. **M. Gu, “Three-dimensional space-invariant point-spread function for a single lens,” J. Opt. Soc. Am. A **12**, 1602–1604 (1995). [CrossRef]

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

**14. **E. Maalouf, “*Contribution to fluorescence microscopy, 3d thick samples deconvolution and depth-variant psf*,” Mulhouse (2010).

**15. **M. Temerinac-Ott, O. Ronneberger, P. Ochs, W. Driever, T. Brox, and H. Burkhardt, “Multiview deblurring for 3-d images from light-sheet-based fluorescence microscopy,” IEEE Transactions on Image Process . **21**, 1863–1873 (2012). [CrossRef]

**16. **S. B. E. Hadj, L. Blancferaud, E. Maalouf, B. Colicchio, and A. Dieterlen, “Depth-variant image restoration in 3d fluorescence microscopy: Two approaches under gaussian and poissonian noise conditions,” in “international symposium on biomedical imaging,” (2012), pp. 1671–1674.

**17. **C. Preza and J. A. Conchello, “Image estimation accounting for point-spread function depth cariation in three-dimensional fluorescence microscopy,” Proc. SPIE - The Int. Soc. for Opt. Eng. **4964**, 135–142 (2003).

**18. **S. Yuan and C. Preza, “Point-spread function engineering to reduce the impact of spherical aberration on 3d computational fluorescence microscopy imaging,” Opt. Express **19**, 23298–23314 (2011). [CrossRef] [PubMed]

**19. **N. Patwary, S. V. King, G. Saavedra, and C. Preza, “Reducing effects of aberration in 3d fluorescence imaging using wavefront coding with a radially symmetric phase mask,” Opt. Express **24**, 12905–12921 (2016). [CrossRef] [PubMed]

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

**21. **B. Kim and T. Naemura, “Blind deconvolution of 3d fluorescence microscopy using depth-variant asymmetric psf,” Microsc. Res. Tech. **79**, 480–494 (2016). [CrossRef] [PubMed]

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

**23. **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. The Opt. Soc. Am. A-optics Image Sci. Vis. **9**, 154–166 (1992). [CrossRef]

**24. **B. Richards and E. E. Wolf, “Electromagnetic diffraction in optical systems. ii. structure of the image field in an aplanatic system,” Proc. The Royal Soc. A: Math. Phys. Eng. Sci. **253**, 358–379 (1959). [CrossRef]

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

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

**27. **E. Maalouf, B. Colicchio, and A. Dieterlen, “Fluorescence microscopy three-dimensional depth variant point spread function interpolation using zernike moments,” J. The Opt. Soc. Am. A-optics Image Sci. Vis. **28**, 1864–1870 (2011). [CrossRef]

**28. **Z. Chen and H. Chen, “The parameter estimation method of gaussian point spread function in microscopic images,” J. Biomed. Eng. **31**, 53 (2014).

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

**30. **R. Cavicchioli, C. Chaux, L. Blancferaud, and L. Zanni, “Ml estimation of wavelet regularization hyperparameters in inverse problems,” in “international conference on acoustics, speech, and signal processing,” (2013), pp. 1553–1557.

**31. **A. Wong, X. Y. Wang, and M. Gorbet, “Bayesian-based deconvolution fluorescence microscopy using dynamically updated nonstationary expectation estimates,” Sci. Reports **5**, 10849 (2015). [CrossRef]

**32. **S. Joshi and M. I. Miller, “Maximum a posteriori estimation with good’s roughness for three-dimensional optical-sectioning microscopy,” J. Opt. Soc. Am. A **10**, 1078–1085 (1993). [CrossRef] [PubMed]

**33. **M. A. Model, J. Fang, P. Yuvaraj, Y. Chen, and B. M. Z. Newby, “3d deconvolution of spherically aberrated images using commercial software,” J. Microsc. **241**, 94–100 (2011). [CrossRef]

**34. **T. H. Besseling, J. Jose, and A. Van Blaaderen, “Methods to calibrate and scale axial distances in confocal microscopy as a function of refractive index,” J. Microsc. **257**, 142–150 (2014). [CrossRef] [PubMed]

**35. **J. A. Conchello, “Fast regularization technique for expectation maximization algorithm for optical sectioning microscopy,” Proc. SPIE - The Int. Soc. for Opt. Eng. **23**, 511–513 (1996).

**36. **G. M. Van Kempen and L. J. Van Vliet, “The influence of the regularization parameter and the first estimate on the performance of tikhonov regularized non-linear image restoration algorithms,” J. Microsc. **198** (Pt 1), 63–75 (2000). [CrossRef] [PubMed]

**37. **G. M. Van Kempen and L. J. Van Vliet, “Background estimation in nonlinear image restoration,” J. Opt. Soc. Am. A-optics Image Sci. Vis. **17**, 425 (2000). [CrossRef]

**38. **N. Dey, L. Blancferaud, C. Zimmer, P. Roux, Z. Kam, J. Olivomarin, and J. Zerubia, “Richardson-lucy algorithm with total variation regularization for 3d confocal microscope deconvolution,” Microsc. Res. Tech. **69**, 260–266 (2006). [CrossRef] [PubMed]

**39. **A. Matakos, S. Ramani, and J. A. Fessler, “Accelerated edge-preserving image restoration without boundary artifacts,” IEEE Transactions on Image Process . **22**, 2019–2029 (2013). [CrossRef]

**40. **S. Liu, J. Li, Z. Zhang, Z. Wang, Z. Tian, G. Wang, and D. Pang, “Fast and high-accuracy localization for three-dimensional single-particle tracking,” Sci. Reports **3**, 2462 (2013). [CrossRef]

**41. **K. Levenberg, “A method for the solution of certain non-linear problems in least squares,” J. Hear. & Lung Transplantation Off. Publ. Int. Soc. for Hear. Transplantation **31**, 436–438 (1944).

**42. **D. W. Marquardt, “An algorithm for least square estimation of non-linear parameters,” J. Soc. for Ind. & Appl. Math. **11**, 431–441 (1963). [CrossRef]

**43. **C. J. Engelbrecht and E. H. K. Stelzer, “Resolution enhancement in a light-sheet-based microscope (spim),” Opt. Lett. **31**, 1477–1479 (2006). [CrossRef] [PubMed]

**44. **P. A. Bello, “Characterization of randomly time-variant linear channels,” IEEE Transactions on Commun. **11**, 360–393 (1963). [CrossRef]

**45. **M. Chen, Y. Li, M. Yang, X. Chen, Y. Chen, F. Yang, S. Lu, S. Yao, T. Zhou, and J. Liu, “A new method for quantifying mitochondrial axonal transport,” Protein & Cell **7**, 804–819 (2016). [CrossRef]