1. Introduction

An image acquisition process can be modeled as a two-dimensional (2-D) convolution of a true image with a point-spread function (PSF) of the imaging system. If the PSF is known a priori, there exist several well-known methods for image restoration, such as inverse filtering, Wiener filtering, and iterative reconstruction (see the excellent review of [1] and the references therein). However, in practice, there are many situations in which PSF is difficult to measure, and only its partial information is available at best [1]. In such cases, we have to estimate both the unknown PSF and the unknown true image. Generally, the problem is termed a blind image deconvolution [2].

Blind image deconvolution methods can be categorized into parametric and non-parametric approaches, depending on their specific assumptions. The parametric approach for optical or electron microscopy exploits the parametric decomposition of the OTF (Optical transfer function, the Fourier transform pair of the PSF) [1, 2], or the location of frequency nulls [3]. The non-parametric approach is more general and does not require any prior knowledge of the PSF. However, the unknown PSF and the unknown true images are alternatively identified using a time consuming iterative update. Classical Lucy-Richardson [4] and Maximum-likelihood (ML) approaches [5] are among this class. Usually, these types of blind deconvolution approaches are computationally extensive, since they require many iterations due to slow convergence [1, 2].

Recently, a revolutionary multichannel blind deconvolution method [6, 7] has been proposed that guarantees the exact reconstruction of the unknown PSF and the true image from multiple blurred images through a number of distinct blur channels. Harikumar and Bresler [6, 7] showed that the perfect estimation of the unknown PSF is almost surely guaranteed if there are at least two or three distinct measurements of the identical image through distinct blur channels. Unlike other blind deconvolution methods, the multichannel blind deconvolution approach does not require the existence of frequency nulls in OTF or parametric decomposition of PSF. Furthermore, the algorithm is in principle non-iterative and computationally and statistically more efficient than the conventional blind deconvolution approaches [8].

However, in many practical situations, it is difficult to obtain multiple distinct blur measurements of an identical scene. For example, in electron-microscopy, the specimens are prone to damage from prolonged exposure, and often drift during the multiple defocus settings [3]. In fluorescent microscopy imaging, the fluorescent intensity could vary due to fast physiological changes between the multiple acquisitions.

The main contribution of this paper is a new non-parametric blind deconvolution theory that eliminates the need of multiple blur measurements, but still guarantees the exact estimation of PSF. The basic idea behind this novel method is to exploit the fact that PSFs in many imaging systems are radially symmetric. This simplifies the PSF estimation problem to a 1-D channel identification problem with multiple excitations, which can be solved by a subspace method [8]. We employ the existing 1-D subspace method in communication theory [8] and modify it to address our 2-D image deconvolution problem.

Since radially symmetric PSFs are quite often encountered in many practical applications, such as optical imaging systems, and electron microscopy, our theory may have great influence on many practical imaging applications. Experimental results using synthetic and real microscopy data also confirm our findings.

2. Multichannel exact blind deconvolution - a review

For a shift invariant point spread function (PSF), an observed noisy image y(r) through a blur channel is described by a 2-D convolution:

where r = (x,y) and ŕ = (x́,ý) denote the 2-D coordinates of the image and measurement domains, respectively; ∗∗ denotes the 2-D convolution; and y(r),h ^2D(r),x(r), and n(r) denote the measured image, the 2-D point-spread function (or blur kernel), the true image, and the additive noise, respectively.

The multichannel blind deconvolution problem can be formulated as recovering x from multiple blur measurements y _(i) through distinct blur kernels h _(i), 1 ≤ i ≤ p,p ≥ 2 [6, 7]:

where n _(i) denotes the noise from i-th blur channel and ∗∗ denotes the 2-D convolution. Harikumar and Bresler [6,7] showed that the multiple blur kernels {h _(i) ^2D}^p _i=1 can be estimated perfectly using subspace methods, after which the unknown image x can be easily reconstructed using the conventional deconvolution methods.

In contrast to the multichannel blind deconvolution theory, our main contribution is to eliminate the necessity of multiple distinct blur measurements but to still guarantee the exact recovery of the blur kernel by exploiting the radial symmetry of the PSF.

3. A novel single channel blind deconvolution theory

3.1. Radially symmetric PSF

Note that due to the nature of the convolution, the nonzero support of (h ^2D ∗∗ x)(r) is usually larger than that of x(r). A full data scenario denotes those cases where our field of view (FOV) includes all the measurements from the extended supports, whereas in a partial data scenario we only measure a part of them. Our theory starts from the full data scenario, and then can be extended to partial data cases.

For the full data case, we can convert the convolution Eq. (1) into the Fourier relationship:

where Y(k),H ^2D(k),X(k) and N(k) are the 2-D Fourier transforms of y(r),h ^2D(r),x(r) and n(r) at the 2-D spatial frequency k = (k_x,k_y), respectively. On the polar coordinate, Eq. (3) can be converted into

where $k=\parallel \mathbf{k}\parallel =\sqrt{k_x^2+k_y^2}\mathrm{and}\Theta {=\tan }^{-1}\frac{k_y}{k_x},$ respectively. Many optical imaging systems such as cameras, optical microscopy, and astigmatism corrected electron microscopy have radially symmetric PSFs. If a PSF is radially symmetric, the optical transfer function (OTF), H ^2D(k) is also radially symmetric. More specifically, since H ^2D(k,Θ) = H(k) for all Θ, Eq. (4) becomes

which implies that the 2-D image convolution can be simplified as a separable 1-D convolution at each Θ:

where ∗ denotes the 1-D convolution, and y ^Θ(r),x ^Θ(r) and n ^Θ(r) correspond to inverse Fourier transforms of Y(k, Θ) ,X(k, Θ) and N(k, Θ) along the radial direction, respectively.

Note the differences between Eq. (2) and Eq. (6). Even though the original 2-D multichannel blind deconvolution approach using Eq. (2) requires the estimation of multiple 2-D filter kernels {h ^2D _(i)}^p _i=1, Eq. (6) only needs identification of a 1-D filter, h, thanks to the power of radial symmetry. This implies that our problem becomes a 1-D blind deconvolution problem, which is much simpler than its 2-D counterpart is. Furthermore, as will be explained, only a single noiseless blur measurement is required to obtain the exact blur kernel.

3.2. Subspace based blind deconvolution

One-dimensional blind deconvolution or identification problems have been extensively investigated in the wireless communication community to deal with the multipath interference [8]. The conventional solution is to use the pilot signal with a training sequence, which sacrifices the communication bandwidth. Blind approaches have been studied to deal with this issue. Among various blind approaches, subspace methods have been investigated extensively due to their computational and statistical efficiency [8]. The key idea of the subspace methods is that an observed signal can be modeled as a vector in the signal subspace spanned by the impulse response function. Since the signal subspace and the noise subspace are orthogonal to each other, the unknown blur kernel can be simply estimated by exploiting the orthogonality [8].

The main contribution of the present paper is to employ the theoretical developments in [8] and apply them to the blind image deconvolution problem by simplifying the 2-D image deconvolution problem as a 1-D channel identification problem. All such developments are possible by exploiting the radial symmetry of the blur kernels. More specifically, in a discrete implementation of Eq. (6), the 1-D blur kernel is a non-causal finite impulse response (FIR) filter given by

since the image blurring is non-causal. The term iΔΘ represents the i-th discretized angle, where ΔΘ is the angular step size and i = 1, ⋯,P. Accordingly, Eq. (6) can be rewritten in a matrix form:

where y ⁽ⁱ⁾ and x ⁽ⁱ⁾ are (N + 2M)×1 measurement vector and N × 1 the true signal vector at the i-th angle, respectively:

and the convolution matrix H ∈ ℂ^(N+2M)×N of the filter kernel h is given by

If the noise is independent of the true signal and has the variance of σ ² _n, then the covariance matrix of the observed signal R _y is given by

where ${\mathbf{R}}_{\mathbf{x}}=\frac{1}{P}{\sum }_{i=1}^P{\mathbf{x}}^{\left(i\right)}{\left({\mathbf{x}}^{\left(i\right)}\right)}^H,$ and the superscript ^H denotes the Hermitian transpose. Let λ₀ ≥ λ₁ ⋯ ≥ λ_N+2M-1 be eigenvalues of R _y. The signal subspace S is then spanned by the eigenvectors corresponding to eigenvalues larger than σ² _n, and the noise subspace G is spanned by the remaining eigenvectors. Under a noiseless scenario, if R _x is full ranked, then the eigenvectors corresponding to λ₀ ≥ λ₁ ⋯ ≥ λ_N-1 span the signal subspace. More specifically, we have

where s_i and g _i denote the eigenvectors of the signal and the noise subspace, respectively. Since the noise and the signal subspace are orthogonal to each other, for a full-ranked R _x, any vectors in the noise space are orthogonal to the columns of the convolution matrix:

Therefore, it is possible to construct a quadratic cost function such that for the true convolution matrix the value of the cost function will be zero. More specifically, the cost function is given by

There is an interesting property of a convolution matrix that can be used to simplify the cost function [8]:

where h denotes the FIR filter given by Eq. (7) and g _i = [g ⁱ ₀, g ⁱ ₁, ⋯ ,g ⁱ _{N + 2M-1}]^T. The matrix ℊ_i can be calculated as follows:

Since our blur kernel h has even symmetry due to the radial symmetry of PSF, we can further impose the symmetry condition:

where ℊ^e _i = (ℊ_i + $\acute{ℊ}$ _i)/2 and $\acute{ℊ}$ _i are constructed by reordering ℊ_i from the last to the first rows. Therefore, Eq. (15) can be reformulated as follows:

where Q = (∑^2M-1 _i=0ℊ^e _iℊ_{eH i}). Under the constraint of ∥h∥₂ = 1, the minimization problem of Eq. (19) is to find the eigenvector Q, which corresponds to the minimum eigenvalue. Additionally, we can define the cost function in terms of the signal subspace [8]:

where Q̃ = (∑^N-1 _i=0 S ^e _i S _{eH i}) and S ^e _i = (S ^e _i + Ś_i)/2 are constructed using S _i, similar to Eq. (17). The optimal h of Eq. (20) can be easily obtained by calculating the eigenvector corresponding to the maximum eigenvalue of Q̃.

The main strength of the subspace method is that if the measurement is free of noise, the solution is exact in the sense that the estimated filter ĥ only differs by a constant scaling factor from the true filter [8]. Furthermore, the filter estimation problem involves finding an eigenvector corresponding to the minimum (or maximum) eigenvalues, which is computationally more efficient compared to other iterative blind identification approaches [4, 5].

4. Implementation issues

The measured signal, {y ⁽ⁱ⁾}^P _i=1, could be obtained using the inverse Fourier transform of the radial lines in Fourier space, {Y(k,iΔΘ)}^P _i=1, along the k direction. However, the direct Fourier transform of an image to obtain the radial lines is prone to error, since the discrete Fourier transform (DFT) approximation of the continuous time Fourier transform (CTFT) assumes a periodic repetition of the images, which breaks the radial symmetry of the OTFs. In addition, there exists considerable errors originating from the interpolation between the Cartesian to polar coordinate transform.

Hence, we implement all the estimation steps in Radon space without ever using the DFT. This is owing to the powerful Fourier slice theory [3], which states that a one-dimensional Fourier transform along the radial direction corresponds to the the radial line of the 2-D continuous Fourier transform. This implies that the 1-D measurement signal y⁽ⁱ⁾ corresponds to the projection data along the i-th projection angle. Hence, we propose the following blind deconvolution algorithm:

Similarly to the multichannel blind deconvolution approach by Harikumar and Bresler [6,7], our single channel exact blind deconvolution approach suffers from some technical difficulties under a partial data measurement scenario, where only a part of blurred image is available [6]. In this case, the Fourier relationship between Eq. (2) and Eq. (3) does not hold, and the 1-D convolution relationship of Eq. (6) is only an approximation of the true blurring process. However, as shown in [6, 7], the performance loss due to the partial data problem is usually limited to the area along the boundary, and we can obtain fairly good reconstruction results in other areas, especially when a penalized least squares algorithm is used for the restoration step.

In our method, we assume that the filter size is known a priori. In practice, however, this should be estimated. There exist several methods for filter size estimation, for example, eigenvalue based techniques [10] or residual based techniques [7]. A detailed discussion of filter size estimation is beyond the scope of this paper, and in this study, we simply chose the filter size by trial and error. Use of efficient algorithms to find the optimal filter size will be addressed elsewhere.

In a noiseless scenario, we can design an inverse filter for image restoration using the estimated blur kernel, ĥ. In practice, measurements are noisy and the inverse filter becomes unstable. A penalized least square algorithm ameliorates the ill-posedness of the restoration process by imposing a penalty term. The resultant cost function can be minimized by any nonlinear optimization algorithm, and we have chosen the ICD algorithm - originally proposed for image reconstruction from noisy sinogram data [9]². In our implementation, the ICD algorithm is used for deblurring the sinogram in Radon space (unlike the original application of ICD in [9]), after which the restored sinogram is filtered and backprojected to obtain the final reconstruction image.

Our algorithm has been implemented using MATLAB (Mathworks, Natick, MA). The ICD routine was implemented using C++ and called from MATLAB using an MEX interface. A Linux computer with an Intel XEON 3.2GHz CPU was used for simulation. The computation time is dependent on the image size. For example, the total computation time for estimating the PSF using an image size of 755 × 755 pixels with 360 views is about 15 seconds, which includes the time to generate the sinogram data (about 10 seconds) and to estimate the filter using the subspace method (about 5 seconds), respectively. One iteration of ICD takes about 2 seconds for 755 × 755 images, and the algorithm converges at around 100 iterations.

5. Experimental results

5.1. Computer simulation

In order to validate the performance of the proposed algorithm, we have performed extensive computer simulations. The original Lena image for simulation is shown in Fig. 1.

We consider two radially symmetric PSFs. First, a circular symmetric Gaussian PSF (see its frequency response in Fig. 2(a)) is calculated as follows:

In Eq. (21), as σ grows large, the image becomes more blurred. We use this model to simulate out-of-focus blur in optical imaging systems. Second, the contrast transfer function (CTF) is considered (see its frequency response in Fig. 2(b)). This is a model to describe the phase contrast in a defocused cryo-electron microscope [3]. Specifically, the radial frequency response H(k) is given by H(k) = CTF(k)E(k), where the envelop function is defined by E(k) = e ^{-Bk 2}

 

Fig. 1. Original Lena image.

Download Full Size | PDF

for some constant B, and the contrast transfer function CTF(k) is given by

where C_s,λ,Δz and Q denote the spherical aberration constant, electron wavelength, defocus, and amplitude constant, respectively [3].

 

Fig. 2. OTFs for (a) a Gaussian blur, and (b) an electron microscope.

Download Full Size | PDF

In order to simulate the blurring process using radially symmetric blurs, we have tested two different approaches. In the first approach, the blurring process is implemented in Radon space. In this case, the blur kernel is a one-dimensional kernel, and the convolution is applied along the radial direction using Eq. (6). This process is summarized as follows:

The second approach is using Fourier domain implementation of a 2-D convolution:

 

Fig. 3. Gaussian blurred image (left column), and proposed blind deconvolution results (right column). Reconstruction at SNR=50dB (1st row), 30dB (2nd row), and 20dB(3rd row), respectively.

Download Full Size | PDF

Ideally, both simulations should provide identical results. However, due to the periodic extension and the interpolation error, the latter approach is more prone to noise. Hence, if the inverse filter is used for restoration from noiseless blur measurements, the latter approach provides inferior reconstruction results due to the approximation error. However, if the penalized least squares is used, we found that the restoration results from both radon and Fourier based blurring provide comparable results thanks to the power of regularization. Hence, we only consider the radon domain approach to simulate the radial blurring.

To evaluate robustness of the algorithm against noise, additive Gaussian noise is added to the blurred image. Noisy images from Gaussian blur and the restored images using the proposed algorithm are illustrated in Fig. 3 at SNR= 50dB, 30dB, and 20dB, respectively. High quality reconstruction was obtained at 50dB, and the reconstruction quality gradually decreases for decreasing SNR. The matrix size for this simulation was 383 × 383.

Similar simulations have been performed for blind estimation of CTF and restoration, as shown in Fig. 4, at SNR = 30dB, 10dB, and 0 dB. Better reconstruction results were obtained even at SNR=30dB compared to that of Gaussian blur. This is because the CTF used in our simulation has non-zero spectral contents for the entire frequency range (see Fig. 2(b)); hence, recovery of the original image is more stable than that of the Gaussian blur.

In order to demonstrate the performance improvement of the proposed algorithm compared to the classical blind-deconvolution method, we have also conducted reconstruction experiments using deconvblind routine in MATLAB, which corresponds to Lucy-Richardson blind deconvolution algorithm [4]. The root-mean square error (RMSE) plots in Fig. 5 demonstrate that our algorithm is better in all SNR ranges. Especially, for the CTF blurred images, the Lucy-Richardson algorithm has failed to provide any reasonable reconstruction result, and our algorithm significantly outperforms the conventional one.

5.2. Optical microscopy experiments

In order to demonstrate the performance of our algorithm in a real measurement scenario, we have applied our algorithm to deconvolution experiments of an optical microscope. More specifically, a slice of mouse brain with a thickness of 60μm is imaged using an Olympus BX51 optical microscope with a magnifying power of 200. During the experiment, we obtained focused and out-of-focus images, as shown in Fig. 6(a) and 6(b), respectively, by changing the focus of the microscope. The in-focus image was used as a ground-truth for evaluating the performance of our blind deconvolution algorithm. Figure 6(c) illustrates the restoration results using the proposed method from the blurred image in Fig. 6(b). Even though the reconstruction is not as perfect as the ground-truth, we can clearly see that the resolution is improved. Since the specimen is 60μm thick, the true blurring process in Fig. 6(b) is three-dimensional, which makes our 2-D assumption invalid. However, the deconvolution result in Fig. 6(c) indicates that our algorithm is still useful in restoring transversal resolution. Another important observation from Fig. 6(c) is the artifact of partial data. As discussed before, artifacts due to the partial data appeared near the boundaries, while the central region of the image was not affected by the artifacts.

5.3. Cryo-electron microscopy experiments

Cryo-electron microscopy (cryo-EM) is a form of electron microscopy (EM) where the sample is studied at cryogenic temperatures using cryofixation [3]. Cryofixation has an important advantage in that the native structure of the sample is preserved by vitrification of the water content. However, due to the low dose data collection using defocus, the collected micrographs are very noisy and are often distorted by the complicated contrast transfer function.

In order to demonstrate that our algorithm is still useful under this extreme measurement condition, we applied our algorithm to the micrographs of 1,900Å diameter Paramecium bursaria chlorella virus, type 1 (PBCV-1) [11] provided by Prof. Timothy S. Baker at the University of California at San Diego. All data was collected using a microscope voltage of 200K and a pixel size of 1.78Å. Figures 7(a) and 7(b) show the defocused micrograph data and the restoration results obtained using the proposed blind deconvolution method, respectively. The image matrix size was 1073 × 1073.

 

Fig. 4. CTF blurred image (left column), and proposed blind deconvolution results (right column). Reconstruction at SNR=30dB (1st row), 10dB (2nd row), and 0dB(3rd row), respectively.

Download Full Size | PDF

Since the original focused image of the specimen is unknown, it is difficult to verify the result. However, we can observe that the outer capsid and inner structures of the virus are more emphasized in the restored image. Considering that the SNR for the cryo-EM is very low, our result indicates that the proposed blind deconvolution algorithm is a still useful tool at low SNR.

 

Fig. 5. Root mean square error (RMSE) of the reconstruction using the classic blind de-convolution (Lucy-Richardson) and the proposed method: (left) Gaussian blur, and (right) CTF blur experiments.

Download Full Size | PDF

 

Fig. 6. Optical microscope images from Olympus BX51 microscope. (a) In-focus image, (b) the defocused image, and (c) blind deconvolution result by the proposed method.

Download Full Size | PDF

 

Fig. 7. Micrograph of PBCV-1 virus [11]: (a) Original de-focused measurement, and (b) the refined image. White arrows are added to indicate some examples of refined areas.

Download Full Size | PDF

6. Conclusion

In this paper, we have proposed a novel single channel blind deconvolution algorithm from a radially symmetric point spread function. Our approach is radically different from the existing blind deconvolution methods, such as that of Lucy-Richardson [4] or the parametric approach [3], since our algorithm is non-iterative and does not require a priori information about the PSF other than radially symmetry. Furthermore, the new method guarantees the exact recovery of unknown blurring kernel in a noiseless measurement scenario. Additionally, compared to the existing multichannel exact blind deconvolution methods, our approach is simpler: since only a single blurred image is required and the estimation process is one-dimensional. Based on the Fourier slice theorem, the blur kernel estimation as well as the restoration procedure were implemented in Radon space, after which the final image was obtained using the filtered backprojection of the restored sinogram. In order to ameliorate the ill-posedness of the deconvolution, a regularized least squares method was employed. Simulation and real experimental results using an optical and an electron microscope demonstrated that reasonably good quality reconstructions were obtained even from noisy measurements.