Single channel blind image deconvolution from radially symmetric blur kernels

The multichannel exact blind image deconvolution theory te lls us thatexact recovery of unknown blur kernels is possible from multiple measurements of an identical scene through distinct blur ch annels. However, in many biological applications, there often exist te chnical difficulties in obtaining multiple distinct blur measurements, since th image content may vary for various reasons, including specimen drift betw e n snapshots, specimen damage due to prolonged exposure, or physiologica changes in live cell imaging. The main contribution of this paper is a ne w non-iterative single channel blind deconvolution method that eliminates the need of multiple blur measurements, but still guarantees an accura te estimation of the blurring kernel. The basic idea behind this novel method is to exploit the radial symmetry of a certain class of PSFs. This approach simplifies the PSF estimation to a 1-D channel identification problem wi th multiple excitations, which can be solved using a standard subspace m ethod. 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 app lications. Simulation results as well as real experimental results usi ng optical and electron microscopy confirm our theory. © 2007 Optical Society of America OCIS codes: (100.1830) Deconvolution; (100.2000) Digital image proces sing; (100.3020) Image reconstruction-restoration. References and links 1. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D flu orescence microscopy images,” IEEE Signal Processing Mag. 23, 32–45 May (2006). 2. D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Mag. 13, 43–64 May (1996). 3. J. Frank,Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996). 4. D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Bl ind deconvolution by means of the RichardsonLucy algorithm,” J. Opt. Soc. Am. A 12, 58–65 (1995). 5. J.A. Conchello, “Superresolution and convergence prope rties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A 15, 2609–2620 (1998). #78978 $15.00 USD Received 16 January 2007; revised 21 March 2007; accepted 22 March 2007 (C) 2007 OSA 2 April 2007 / Vol. 15, No. 7 / OPTICS EXPRESS 3791 6. G. Harikumar and Y. Bresler, “Exact image deconvolution fro m multiple FIR blurs,” IEEE Trans. on Image Processing 8, 846–862, (1999). 7. G. Harikumar and Y. Bresler, “Perfect blind restoration of images blurred by multiple filters: Theory and efficient algorithms,” IEEE Trans. on Image Processing 8, 202–219 (1999). 8. E. Moulines, P. Duhamel, J.-F. Cardoso, and S. Mayrargue, “ Subspace methods for the blind identification of multichannel FIR filters,”IEEE Trans. on Signal Processing 43, 516–525 (1995). 9. C. Bouman, and K. Sauer, “A unified approach to statistical t omography using coordinate descent optimization,” IEEE Trans. on Image Processing, 5, 480–492 (1996). 10. M. Gurelli, and C. Nikias, “EVAM: An eigenvector-based a lgorithm for multichannel blind deconvolution of input colored signals,”IEEE Trans. on Signal Processing, 43, 134–149 (1995). 11. X. Yan, N. Olson, J. Van Etten, M. Bergoin, M. Rossmann, and T. Baker, “Structure and assembly of large lipid-containing dsDNA viruses,”Nat. Struct. Biol., 7, 101–103 (2000).


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.

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 r ′ = (x ′ , y ′ ) 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 2D (i) } 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.

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 = ||k|| = k 2 x + k 2 y and Θ = tan −1 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.

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 (i) and x (i) are (N + 2M) × 1 measurement vector and N × 1 the true signal vector at the i-th angle, respectively: If the noise is independent of the true signal and has the variance of σ 2 n , then the covariance matrix of the observed signal R y is given by where

H
, and the superscript H denotes the Hermitian transpose.Let The signal subspace S is then spanned by the eigenvectors corresponding to eigenvalues larger than σ 2 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 λ 0 ≥ λ 1 • • • ≥ λ 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 The matrix G i can be calculated as follows: where G e i = (G i + G ′ i )/2 and G ′ i are constructed by reordering G i from the last to the first rows.Therefore, Eq. ( 15) can be reformulated as follows: where . Under the constraint of ||h|| 2 = 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 and S e i = (S i + S ′ i )/2 and S i 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].

Implementation issues
The measured signal, {y (i) } 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 (i) corresponds to the projection data along the i-th projection angle.Hence, we propose the following blind deconvolution algorithm: 1. Apply a Radon transform to obtain the measured sinogram data {y (i) } P i=1 .
2. Estimate the filter kernel h using the sinogram by minimizing Eq. (20).
3. Obtain a restored sinogram by solving the inverse problem of Eq. (8) using a penalized least squares algorithm.In this paper, we employ the iterative coordinate descent (ICD) algorithm by Bouman and Sauer [9].
4. Calculate a restored image using a filtered backprojection of the restored sinogram.
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]2 .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.

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 where the envelop function is defined by for some constant B, and the contrast transfer function CT F(k) is given by where C s , λ , ∆z and Q denote the spherical aberration constant, electron wavelength, defocus, and amplitude constant, respectively [3].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: 1. Apply Radon transform of a true image at every ∆Φ angle to obtain the sinogram of a true image.
2. Apply 1-D convolution along the radial direction of the sinogram to obtain a blurred sinogram.
3. Apply an inverse Radon transform to obtain a blurred image.
The second approach is using Fourier domain implementation of a 2-D convolution: 1. Zero-pad the true image.
2. Multiply the frequency domain response of the zero padded image with 2-D OTF.
3. Take the inverse Fourier transform to obtain the blurred image.
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.

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.

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

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.

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