OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editor: Gregory W. Faris
  • Vol. 2, Iss. 5 — May. 17, 2007
« Show journal navigation

Single channel blind image deconvolution from radially symmetric blur kernels

Kwang Eun Jang and Jong Chul Ye  »View Author Affiliations


Optics Express, Vol. 15, Issue 7, pp. 3791-3803 (2007)
http://dx.doi.org/10.1364/OE.15.003791


View Full Text Article

Acrobat PDF (2296 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

The multichannel exact blind image deconvolution theory tells us that exact recovery of unknown blur kernels is possible from multiple measurements of an identical scene through distinct blur channels. However, in many biological applications, there often exist technical difficulties in obtaining multiple distinct blur measurements, since the image content may vary for various reasons, including specimen drift between snapshots, specimen damage due to prolonged exposure, or physiological changes in live cell imaging. The main contribution of this paper is a new non-iterative single channel blind deconvolution method that eliminates the need of multiple blur measurements, but still guarantees an accurate 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 with multiple excitations, which can be solved using a standard subspace method. 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. Simulation results as well as real experimental results using optical and electron microscopy confirm our theory.

© 2007 Optical Society of America

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

01. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Processing Mag. 23,32–45 May (2006). [CrossRef]

] 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

01. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Processing Mag. 23,32–45 May (2006). [CrossRef]

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

02. D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Mag. 13,43–64 May (1996). [CrossRef]

].

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

01. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Processing Mag. 23,32–45 May (2006). [CrossRef]

, 2

02. D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Mag. 13,43–64 May (1996). [CrossRef]

], or the location of frequency nulls [3

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

04. D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A 12,58–65 (1995). [CrossRef]

] and Maximum-likelihood (ML) approaches [5

05. J.A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A 15,2609–2620 (1998). [CrossRef]

] are among this class. Usually, these types of blind deconvolution approaches are computationally extensive, since they require many iterations due to slow convergence [1

01. P. Sarder and A. Nehorai, “Deconvolution methods for 3-D fluorescence microscopy images,” IEEE Signal Processing Mag. 23,32–45 May (2006). [CrossRef]

, 2

02. D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Mag. 13,43–64 May (1996). [CrossRef]

].

Recently, a revolutionary multichannel blind deconvolution method [6

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

, 7

07. 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). [CrossRef]

] 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

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

, 7

07. 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). [CrossRef]

] 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

08. 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). [CrossRef]

].

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

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

08. 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). [CrossRef]

]. We employ the existing 1-D subspace method in communication theory [8

08. 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). [CrossRef]

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

y(r)=h2D(rr´)x(r´)+n(r)=(h2Dx)(r)+n(r)
(1)

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 ≤ ip,p ≥ 2 [6

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

, 7

07. 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). [CrossRef]

]:

y(i)=h(i)2Dx+n(i),
(2)

where n (i) denotes the noise from i-th blur channel and ∗∗ denotes the 2-D convolution. Harikumar and Bresler [6

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

,7

07. 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). [CrossRef]

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

Y(k)=H2D(k)X(k)+N(k),
(3)

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 = (kx,ky), respectively. On the polar coordinate, Eq. (3) can be converted into

Y(k,Θ)=H2D(k,Θ)X(k,Θ)+N(k,Θ)
(4)

where k=k=kx2+ky2andΘ=tan1kykx, 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

Y(k,Θ)=H(k)X(k,Θ)+N(k,Θ),
(5)

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

yΘ=hxΘ+nΘ,
(6)

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

08. 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). [CrossRef]

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

08. 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). [CrossRef]

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

08. 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). [CrossRef]

].

The main contribution of the present paper is to employ the theoretical developments in [8

08. 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). [CrossRef]

] 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

h=[hMh0hM]T,
(7)

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:

y(i)=Hx(i)+n(i),
(8)

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:

y(i)=[yM(i)y0(i)yN+M1(i)]T,x(i)=[x0(i)xN1(i)]T
(9)

and the convolution matrix H ∈ ℂ(N+2MN of the filter kernel h is given by

H=[hM00hMhM0hM](N+2M)×N.
(10)

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

Ry=1Pi=1Py(i)(y(i))H=HRxHH+σn2I
(11)

where Rx=1Pi=1Px(i)(x(i))H, and the superscript H denotes the Hermitian transpose. Let λ0 ≥ λ1 ⋯ ≥ λN+2M-1 be eigenvalues of R y. 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

S=[s0s1sN1]
(12)
G=[g0g1g2M1]
(13)

where si 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:

giHH=0,i=0,,2M1.
(14)

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

c(H)=i=02M1giHH2=i=02M1giHHHHgi.
(15)

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

08. 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). [CrossRef]

]:

giHH=hHi
(16)

where h denotes the FIR filter given by Eq. (7) and g i = [g i 0, g i 1, ⋯ ,g i N + 2M-1]T. The matrix ℊi can be calculated as follows:

i=[g0ig1igN1ig1ig2igNig2Mig2M+1igN+2M1i]
(17)

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

giHH=hHie
(18)

where ℊe i = (ℊi + ´ i)/2 and ´ i are constructed by reordering ℊi from the last to the first rows. Therefore, Eq. (15) can be reformulated as follows:

c(h)=hTQh
(19)

where Q = (∑2M-1 i=0e ieH i). Under the constraint of ∥h2 = 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

08. 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). [CrossRef]

]:

c(h)=Nh2hHQ˜h
(20)

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

08. 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). [CrossRef]

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

04. D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A 12,58–65 (1995). [CrossRef]

, 5

05. J.A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A 15,2609–2620 (1998). [CrossRef]

].

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

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

    09. C. Bouman and K. Sauer, “A unified approach to statistical tomography using coordinate descent optimization,” IEEE Trans. on Image Processing ,5,480–492 (1996). [CrossRef]

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

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

,7

07. 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). [CrossRef]

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

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

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

06. G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

, 7

07. 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). [CrossRef]

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

10. M. Gurelli and C. Nikias, “EVAM: An eigenvector-based algorithm for multichannel blind deconvolution of input colored signals,” IEEE Trans. on Signal Processing ,43,134–149 (1995). [CrossRef]

] or residual based techniques [7

07. 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). [CrossRef]

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

09. C. Bouman and K. Sauer, “A unified approach to statistical tomography using coordinate descent optimization,” IEEE Trans. on Image Processing ,5,480–492 (1996). [CrossRef]

]2. In our implementation, the ICD algorithm is used for deblurring the sinogram in Radon space (unlike the original application of ICD in [9

09. C. Bouman and K. Sauer, “A unified approach to statistical tomography using coordinate descent optimization,” IEEE Trans. on Image Processing ,5,480–492 (1996). [CrossRef]

]), 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:

h(x,y)=12πσ2exp(x2+y22σ2).
(21)

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

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

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

CTF(k)=sin(γ(k)+ψ)
(22)
γ(k)=π2(Csλ3k42Δzλk2)
(23)
ψ=tan1(Q/1Q2)
(24)

where Cs,λz and Q denote the spherical aberration constant, electron wavelength, defocus, and amplitude constant, respectively [3

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

].

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

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:

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

04. D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A 12,58–65 (1995). [CrossRef]

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

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

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). [CrossRef] [PubMed]

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

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

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

04. D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A 12,58–65 (1995). [CrossRef]

] or the parametric approach [3

03. J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

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

Acknowledgments

This work was supported by grant No. R01-2005-000-10035-0 from the Basic Research Program of the Korea Science and Engineering Foundation. The author would like to thank Prof. Timothy S. Baker and Dr. Bob Sinkovits at UCSD for providing PBCV-1 virus micrographs.

Footnotes

2Note that the other computationally more efficient optimization methods rather than the ICD could be also used for deblurring the sinogram data.

References and links

01.

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

02.

D. Kundur and D. Hatzinakos, “Blind image deconvolution,” IEEE Signal Processing Mag. 13,43–64 May (1996). [CrossRef]

03.

J. Frank, Three Dimensional Electron Microscopy of Macromolecular Assemblies, (Academic Press, San Diego, CA, USA, 1996).

04.

D. A. Fish, A. M. Brinicombe, E. R. Pike, and J. G. Walker, “Blind deconvolution by means of the Richardson-Lucy algorithm,” J. Opt. Soc. Am. A 12,58–65 (1995). [CrossRef]

05.

J.A. Conchello, “Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images,” J. Opt. Soc. Am. A 15,2609–2620 (1998). [CrossRef]

06.

G. Harikumar and Y. Bresler, “Exact image deconvolution from multiple FIR blurs,” IEEE Trans. on Image Processing 8,846–862, (1999). [CrossRef]

07.

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). [CrossRef]

08.

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). [CrossRef]

09.

C. Bouman and K. Sauer, “A unified approach to statistical tomography using coordinate descent optimization,” IEEE Trans. on Image Processing ,5,480–492 (1996). [CrossRef]

10.

M. Gurelli and C. Nikias, “EVAM: An eigenvector-based algorithm for multichannel blind deconvolution of input colored signals,” IEEE Trans. on Signal Processing ,43,134–149 (1995). [CrossRef]

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). [CrossRef] [PubMed]

OCIS Codes
(100.1830) Image processing : Deconvolution
(100.2000) Image processing : Digital image processing
(100.3020) Image processing : Image reconstruction-restoration

ToC Category:
Image Processing

History
Original Manuscript: January 16, 2007
Revised Manuscript: March 21, 2007
Manuscript Accepted: March 22, 2007
Published: April 2, 2007

Virtual Issues
Vol. 2, Iss. 5 Virtual Journal for Biomedical Optics

Citation
Kwang Eun Jang and Jong Chul Ye, "Single channel blind image deconvolution from radially symmetric blur kernels," Opt. Express 15, 3791-3803 (2007)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-15-7-3791


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  1. P. Sarder and A. Nehorai, "Deconvolution methods for 3-D fluorescence microscopy images," IEEE Signal Processing Mag.23, 32-45 May (2006). [CrossRef]
  2. D. Kundur and D. Hatzinakos, "Blind image deconvolution," IEEE Signal Processing Mag.13, 43-64 May (1996). [CrossRef]
  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, "Blind deconvolution by means of the Richardson- Lucy algorithm," J. Opt. Soc. Am. A 12, 58-65 (1995). [CrossRef]
  5. J. A. Conchello, "Superresolution and convergence properties of the expectation-maximization algorithm for maximum-likelihood deconvolution of incoherent images," J. Opt. Soc. Am. A 15, 2609-2620 (1998). [CrossRef]
  6. G. Harikumar and Y. Bresler, "Exact image deconvolution from multiple FIR blurs," IEEE Trans. on Image Processing 8, 846-862 (1999). [CrossRef]
  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). [CrossRef]
  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). [CrossRef]
  9. C. Bouman, and K. Sauer, "A unified approach to statistical tomography using coordinate descent optimization," IEEE Trans. on Image Processing 5, 480-492 (1996). [CrossRef]
  10. M. Gurelli and C. Nikias, "EVAM: An eigenvector-based algorithm for multichannel blind deconvolution of input colored signals," IEEE Trans. on Signal Processing  43, 134-149 (1995). [CrossRef]
  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). [CrossRef] [PubMed]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited