Measure and model a 3-D space-variant PSF for fluorescence microscopy image deblurring

Yemeng Chen; Mengmeng Chen; Li Zhu; Jane Y. Wu; Sidan Du; Yang Li

doi:10.1364/OE.26.014375

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.

Fig. 1 The optical path difference (OPD) for the point source located at different depths under the same condition in WFM.

Download Full Size | PDF

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.

Fig. 2 A framework of proposed SV deblurring for fluorescence microscopy images. Sparsely distributed fluorescent beads were recorded and used to generate a dense estimation of the SV-PSF. The PSF model considered the axial and radial variation. The proposed space-variant restoration algorithm considered the matching object point with its corresponding PSF by using an alternating optimization scheme.

Download Full Size | PDF

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} : = {\tilde{h_{i}} (s) | i = 1, \dots, N}$ , where $\tilde{h_{i}} (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} : = {\tilde{p_{i}} = (\tilde{x_{i}}, \tilde{y_{i}}, \tilde{z_{i}}) | i = 1, \dots, N}$ where ${\tilde{p}}_{n}$ refers to the fine center position of $\tilde{h_{i}} (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_i refer to the center position of the center-adjusted PSF h_i (s), to prepare for the PSF modeling. See Fig. 3 for the diagram of the PSF estimation process.

Fig. 3 A diagram illustrating the process of SV-PSF estimation. The PSF at any specific depth is defined as a weighted average of all its approximations from measurements using the similarity transformation. Stage drift was also considered in most cases.

Download Full Size | PDF

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 $T (\cdot)$ that allows a measured PSF h_i (s) to express the PSF h_j (s) at any specific position.

h_{j} (s) = T (h_{i} (s))

where h_i (s) refers to the intensity of light at s = (u, v, w)^T of the measured PSF.

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 $T (\cdot)$ that allows h_i (s) to approximate h_j (s) was defined by a similarity transformation,

h_{j} (s) \approx T (h_{i} (s); θ_{i j}) = h_{i} (θ_{i j} s)

where θ_ij is the similarity transformation matrix,

θ_{i j} = d i a g (α_{i j}, α_{i j}, β_{i j})

where (α_ij, β_ij) refer to the scaling parameters which transform h_i (s) to approximate h_j (s) in the radial and axial direction, respectively. Then, the residual error is

C (θ_{i j}) = \frac{1}{2} {‖ h_{i} (θ_{i j} s) - h_{j} (s) ‖}^{2}

where ‖·‖ refers to the L2-norm and the parameter θ_ij can be computed as

{\hat{θ}}_{i j} = \underset{θ_{i j}}{\arg \min} C (θ_{i j})

To solve the nonlinear fitting problem, we minimize the least-square criterion using the Levenberg-Marquardt algorithm [41,42]. Table 1 gives the derivative expression for each parameter in θ_ij and

\frac{\partial h_{i} (θ_{i j} s)}{\partial (α u)}

,

\frac{\partial h_{i} (θ_{i j} s)}{\partial (α v)}

,

\frac{\partial h_{i} (θ_{i j} s)}{\partial (β w)}

refer to the gradient of h_i (θ_ij s) along the x-axis, y-axis and z-axis, respectively.

Table 1. Derivative expressions for similarity-transformation based PSF model.

View Table | View all tables in this article

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_r and k_z 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.

Fig. 4 PSF model evaluation with synthetic data (a)–(c), (g)–(i) and experimental data (d)–(f). (a)(d) Coefficients of similarity transformation with different relative depths. (b)(e) MSE between PSFs with different relative depths. (c)(f) Comparison of our proposed similarity transformation based PSF model (STM) with the piecewise-invariant model and the blending model. PSFs in sub-region 1 are approximated with the PSF of z_p = 1000nm and PSFs in sub-region 2 are approximated with the PSF of z_p = 3500nm (z_p refers to the axial location of the object point). The synthetic data was generated by PSF Generator, the ImageJ plugin with acquisition parameters: the emission wavelength λ=530nm, the numerical aperture of the microscope NA = 1.4, the refractive index of the oil-immersion layer n_i = 1.515, the refractive index of the specimen layer n_s = 1.3, the working distance t_i = 150µm and the voxel size 25 × 25 × 100nm³. Besides, the particle position is set from 1000 nm to 3500 nm whose step is 50 nm. S.T. refers to the acronym of similarity transformation. The real data was acquired by Leica SP8 microscope system and a NA 1.4 oil-immersion objective under wide-field mode with the same λ, n_i and voxel size as the synthetic data. (g) Slopes for coefficients of similarity transformation along with different depths. (h-i) Comparison between our proposed similarity transformation based PSF model (STM) with the PCA based model of approximating PSFs at 1 − 3.5 µm with step equal to 100 nm using six synthetic PSFs at 1, 1.5, 2, 2.5, 3, 3.5 µm for parameter fitting.

Download Full Size | PDF

To reduce error, the proposed method uses all of the PSF measurements efficiently. The PSF estimation $\hat{h} (s; s_{0})$ whose point source is located at the coordinate s₀ = (u₀, v₀, w₀)^T in the image domain, is defined as a weighted average of all its approximations $T (h_{i} (s); θ_{i})$ from measurements using the similarity transform.

\hat{h} (s; s_{0}) = A \sum_{i = 1}^{N} ω_{i} T (h_{i} (s); θ_{i})

where θ_i refers to the similarity transformation matrix that allows h_i (s) to approximates

\hat{h} (s; s_{0})

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

ω_{i} = {\begin{matrix} 1 / | z - z_{i} |, & | z - z_{i} | \leq Δ z_{t h r e s h} \\ 0, & otherwise \end{matrix}

where z refers to the axial position of the object point s₀, z_i refers to the axial center position of h_i (s), and Δz_thresh refers to the threshold of the relative depth Δ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.

Fig. 5 One example of PSF estimation using different methods on synthetic data (a) and real data (b). The ground truth of PSF located at the specific depth (first column), the estimation result with the piecewise-invariant PSF model (second column), the estimation result with the blending PSF model (third column) and the estimation result with the proposed SV similarity transformation based PSF model (fourth column) are shown in the corresponding panels.

Download Full Size | PDF

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.

Fig. 6 Influence of stage-drift. (a) & (c) Stage drift paths along the y-axis and x-axis. Blue lines represent the drift paths of the PSFs extracted from one bead image, and red lines are the drift paths of PSFs extracted from another bead image. For observability, the stage drift is amplified 100 times. (b) & (d) Meridional cross-section of the PSF before (upper) and after (lower) stage drift correction along the y-axis (b) and x-axis (d). Red lines refer to the projection of the drift path along the y-axis (b) and x-axis (d). The PSF size is 83 × 83 × 61, and the voxel size is 25 × 25 × 100nm³.

Download Full Size | PDF

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

S H (s; x_{0}, z_{0}) = I H (s; z_{0}) * I D (s; x_{0})

where SH (s; x₀, z₀) is the system PSF located at x = x₀, z = z₀, IH (s; z₀) refers to the intrinsic PSF at depth z = z₀, and ID (s; x₀) refers to the intensity distribution of its excitation beam located at the cross-section x = x₀.

Fig. 7 Spatial variation for the PSF in LSFM, wherein the system PSF is estimated using the pointwise product of the intrinsic PSF and the intensity distribution of its excitation beam.

Download Full Size | PDF

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.

(f \otimes h) (s) = \sum_{m = 1}^{M} f_{m} (s) * h_{m} (s)

where M is the number of strata and f_m (s) contains the mth stratum of the object.

f_{m} (s) = {\begin{matrix} f (s) & for z = m \\ 0 & otherwise \end{matrix}

* denotes the 3D convolution operation and h_m (s) is the system response to a point source placed at the mth 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

\begin{matrix} (f \otimes h) (s) = A \sum_{i = 1}^{N} \sum_{m = 1}^{M} ω_{i} (f_{m} (s) * h_{i} (θ_{i}^{m} s)) \\ = A \sum_{i = 1}^{N} \sum_{m = 1}^{M} ω_{i} F^{- 1} (F (f_{m} (s)) \cdot F (h_{i} (θ_{i}^{m} s))) \\ A \sum_{i = 1}^{N} \sum_{m = 1}^{M} \frac{ω_{i}}{| θ_{i}^{m} |} F^{- 1} (F (f_{m} (s)) \cdot F (h_{i} (s))) \end{matrix}

Thus, the Fourier transformation of 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

f^{(k + 1)} (s) = \sum_{m = 1}^{M} {[\frac{g (s)}{(f^{(k)} \otimes h) (s)}] * h_{m} (- s)} f_{m}^{(k)} (s)

where

f_{m}^{(k)} (s)

refers to the mth 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

J (f) = \sum_{s \in S} {- g (s) \ln [(f \otimes h) (s)] + (f \otimes h) (s)}

As h(s; s₀) in 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

z^{(j + 1)} = z^{(j)} - {[\frac{\partial J (f, z)}{\partial z} / \frac{\partial^{2} J (f, z)}{\partial z^{2}}] |}_{z = z^{(j)}}

The partial derivative and the second partial derivative of J (f, z) with respect to depth z is shown as

\frac{\partial J (f, z)}{\partial z} = \sum_{s \in S} [f (- s) * (1 - \frac{g (s)}{(f \otimes h) (s)})] \frac{\partial h (s; s_{0})}{\partial z}

\frac{\partial^{2} J (f, z)}{\partial z^{2}} = \sum_{s \in S} {[f (s) * (1 - \frac{g (s)}{(f \otimes h) (s)})] \frac{g (s) f (- s)}{{[(f \otimes h) (s)]}^{2}} {[\frac{\partial h (s; s_{0})}{\partial z}]}^{2} + [f (- s) * (1 - \frac{g (s)}{(f \otimes h) (s)})] \frac{\partial^{2} h (s; s_{0})}{\partial z^{2}}}

where the partial derivative and the second partial derivative of the PSF h (s; s₀) with respect to depth z is written as

\frac{\partial h (s; s_{0})}{\partial z} = A \sum_{i = 1}^{N} ω_{i} [(\frac{k_{r} \partial}{\partial (α u)} u + \frac{k_{r} \partial}{\partial (α v)} v + \frac{k_{z} \partial}{\partial (β w)} w) h_{i} (θ_{i} s)]

\frac{\partial^{2} h (s; s_{0})}{\partial z^{2}} = A \sum_{i = 1}^{N} ω_{i} [(\frac{k_{r} \partial}{\partial (α u)} u + \frac{k_{r} \partial}{\partial (α v)} v + \frac{k_{z} \partial}{\partial (β w)} w) h_{i} (θ_{i} s)]

where k_r and k_z refer to the computed slope of the scaling parameters along the radial and axial directions. In this manner, the axial location of sample 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 (f) + J_{T V} (f) = \sum_{s \in S} {- g (s) \ln [(f \otimes h) (s)] + (f \otimes h) (s)} + λ_{T V} \sum_{s \in S} | \nabla f (s) |

where J_TV (f) refers to the TV regularization term, λ_TV 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

f^{(k + 1)} (s) = \sum_{m = 1}^{M} \frac{f_{m}^{(k)} (s_{0})}{1 - λ_{T V} d i v (\frac{\nabla f^{(k)} (s)}{| \nabla f^{(k)} (s) |})} {[\frac{g (s)}{(f^{(k)} \otimes h) (s)}] * h_{m} (- s)}

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.

Fig. 8 Simulations for deconvolution of synthetic data using the space-invariant and space-variant methods in widefield system. (a)–(c) The raw image (upper-left), the degraded image blurred by a series of space-variant PSFs and Poisson noise (lower-left), the restored image using SI-MLTV (upper-right) and the restored image using SV-MLTV (lower-right). Dotted lines in the XZ section (upper) show where the XY sections are taken. The lower-left XY section is taken where the upper dotted line is located, and the lower-right XY section is taken where the lower dotted line is located. The image size is 199 × 199 × 99, and the voxel size is 25 × 25 100nm³. (a) Deconvolution of three balls, where the intensity of the raw image is 255 and the background is 1. The PSNR is 21.2872. (b) Deconvolution of a cavity, where the intensities of a light spot inside, at the central portion, shell and background are 255, 240, 200 and 1, respectively. The PSNR is 23.2761. (c) Deconvolution of the multi-objects with different intensities: 255 for the ball, 211 for the cube, 180 for the cylinder, 225 for the triangle, 150 for cross and 1 for the background. The PSNR is 24.8370.

Download Full Size | PDF

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

Fig. 9 Simulation for deconvolution of complicated images of neuronal axons in primary neuronal cultures acquired using a Leica SP8 confocal system: (a) the z-stack of the degraded image blurred by space-variant PSFs and Poisson noise with PSNR of 21.5493, (b) the z-stack of the restored image using SV-MLTV, (c) the z-projection (upper) and its partial views (lower) of the restored image using SI-MLTV, (d) the z-projection (upper) and its partial views (lower) of the restored image using SV-MLTV. The raw image size is 601 × 601 × 8, and the voxel size is 50 × 50 × 200nm³. Before image degradation, the raw image is scaled into 601 × 601 × 40 using the software ImageJ by setting the option named “Z Scale” equal to 5. PSF size is 33 × 33 × 25 and the voxel size is 50 × 50 × 40nm³.

Download Full Size | PDF

Fig. 10 Deconvolution of the complicated image where the raw sample is the measured confocal image acquired with the Leica SP8. (a) Lateral cross-sections of the raw sample image. (b) Lateral cross-sections of the degraded image blurred by space-variant PSFs and Poisson noise with PSNR of 21.5493. (c) Lateral cross-sections of the restored image using SI-MLTV. (d) Lateral cross-sections of the restored image using SV-MLTV. The raw image size is 601 × 601 × 8, and the voxel size is 50 × 50 × 200nm³.

Download Full Size | PDF

Fig. 11 Deconvolution of real data acquired using a Leica SP8 system in wide-field mode: (a) the z-stack of the degraded image, (b) the z-stack of the restored image using SV-MLTV, (c) the z-projection (upper) and its partial views (lower) of the restored image using SI-MLTV, (d) the z-projection (upper) and its partial views (lower) of the restored image using SV-MLTV.

Download Full Size | PDF

Fig. 12 Deconvolution of the synthetic image of mitochondria: (a) the raw images, (b) the degraded images blurred by space-variant PSFs and Poisson noise with PSNR of 29.2674, (c) the restored images using SI-MLTV, and (d) the restored images using SV-MLTV. The meridional cross-sections of the PSF along the y-axis (upper panels) and x-axis (lower panels). The image size is 199 × 199 × 29, and the voxel size is 50 × 50 200nm³. SI-MLTV: space-invariant ML algorithm with TV regularization; SV-MLTV: space-variant ML algorithm with TV regularization.

Download Full Size | PDF

Table 2. Quantitative analyses of particle detection.

View Table | View all tables in this article

Table 3. Quantitative analyses of deconvolution with different methods.

View Table | View all tables in this article

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]

Data	Real no.	Detected no.	Matched no.	Precision	Recall	F-value
SI-MLTV	1500	1495	456	0.3050	0.3040	0.3045
SV-MLTV	1500	1500	1497	0.9980	0.9980	0.9980

Samples		Tri-balls	Cavity	Multi-obj	Confocal	Mito
MSE	SI-MLTV	128.0413	142.3446	67.2999	504.1297	374.5518
	SV-MLTV	54.9510	28.1196	10.8789	377.5328	183.5430
	%	57.1	80.3	83.8	25.1	51.0
PSNR / dB	SI-MLTV	27.0573	26.5974	29.8507	21.1054	22.3957
	SV-MLTV	30.7310	33.6407	37.7650	22.3613	25.4934
	%	13.6	26.5	26.5	6.0	13.8
I-divergence	SI-MLTV	2.7013	2.8292	2.4247	37.8429	62.4136
	SV-MLTV	1.6458	1.5885	1.1278	28.2951	26.0499
	%	39.1	43.9	53.5	25.2	58.3

Data	Real no.	Detected no.	Matched no.	Precision	Recall	F-value
SI-MLTV	1500	1495	456	0.3050	0.3040	0.3045
SV-MLTV	1500	1500	1497	0.9980	0.9980	0.9980

Samples		Tri-balls	Cavity	Multi-obj	Confocal	Mito
MSE	SI-MLTV	128.0413	142.3446	67.2999	504.1297	374.5518
	SV-MLTV	54.9510	28.1196	10.8789	377.5328	183.5430
	%	57.1	80.3	83.8	25.1	51.0
PSNR / dB	SI-MLTV	27.0573	26.5974	29.8507	21.1054	22.3957
	SV-MLTV	30.7310	33.6407	37.7650	22.3613	25.4934
	%	13.6	26.5	26.5	6.0	13.8
I-divergence	SI-MLTV	2.7013	2.8292	2.4247	37.8429	62.4136
	SV-MLTV	1.6458	1.5885	1.1278	28.2951	26.0499
	%	39.1	43.9	53.5	25.2	58.3

Measure and model a 3-D space-variant PSF for fluorescence microscopy image deblurring

Abstract

1. Introduction

2. General framework of space-variant deblurring

3. Estimating the SV-PSF from sparse measurements

3.1. Modeling and estimating the axial variation of PSF

3.2. Stage drift estimation and correction

3.3. The radial variation of the PSF in LSFM

4. 3-D space-variant restoration

5. Evaluation

6. Conclusions

Funding

References and links

Cited By

Figures (12)

Tables (3)

Equations (20)

Optics Express

Parameter	Derivative expressions
α	$\frac{\partial C (θ_{i j})}{\partial α} = ‖ h_{i} (θ_{i j} S) - h_{j} (S) ‖ (\frac{\partial h_{i} (θ_{i j} S)}{\partial (α u)} u + \frac{\partial h_{i} (θ_{i j} S)}{\partial (α v)} v)$
β	$\frac{\partial C (θ_{i j})}{\partial β} = ‖ h_{i} (θ_{i j} S) - h_{j} (S) ‖ \frac{\partial h_{i} (θ_{i j} S)}{\partial (β w)} w$