IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. X, NO. XX, XXXX 2006

1

Atmospheric Turbulence Degraded Image Restoration using Principal Components Analysis Dalong Li, Russell M. Mersereau, Fellow, IEEE and Steven Simske, Member, IEEE

Abstract— Our earlier work revealed a connection between blind image deconvolution and Principal Components Analysis (PCA). In this letter, we explicitly formulate multichannel and single-channel blind image deconvolution as a PCA problem. Although PCA is derived from blur models that do not contain additive noise, it can be justified both on theoretical and experimental grounds that the PCA-based restoration algorithm is actually robust to the presence of white noise. The algorithm is applied to the restoration of atmospheric turbulence degraded imagery and compared to an adaptive Lucy-Richardson maximum likelihood (LR) algorithm on both real and simulated atmospheric turbulence blurred images. It is shown that the PCAbased blind image deconvolution runs faster and is more robust to noise. Index Terms— Principal Components Analysis, blind image deconvolution, Lucy-Richardson algorithm, atmospheric turbulence.

I. INTRODUCTION It is well known that atmospheric turbulence degrades the quality of long-distance surveillance imagery [1]. Atmospheric turbulence blur is caused primarily by the random fluctuations of the refraction index. These fluctuations, which can be modeled as a dynamic random process, perturb the phase of the incoming light. The restoration of atmospheric turbulence degraded images has been actively studied [2-7]. Based on the refraction index structure functions, Hufnagel and Stanley [2] derived a long-exposure optical transfer function (OTF) 2 2 5/6 H(u, v) = e−λ(u +v ) (1) to model the long-term effect of turbulence in optical imaging. Here u and v are the horizontal and vertical frequency variables and λ parameterizes the severity of the turbulence blur. Since the refraction index fluctuation is a random process, the blurring parameter λ is often unknown. In such situations, blind image deconvolution algorithms [3] are applicable. Lucy [4] and Richardson [5] independently developed a non-blind deconvolution method based on maximum likelihood estimation. The algorithm maximizes the likelihood that the resulting image, when convolved with the blurring point spread function (PSF), is an instance of the blurred image. This method can be effective when the PSF is known, but little is known about any additive noise. Ayers and Dainty Manuscript received April 19, 2006; revised July 11, 2006. Dalong Li and Russell M. Mersereau are with Center for Signal and Image Processing, School of Electrical and Computer Engineering at the Georgia Institute of Technology, Atlanta, GA 30332 (e-mail: [email protected]). Steven Simske is with the Digital Printing and Imaging Laboratory, HewlettPackard Laboratories, Fort Collins, CO 80528

[6] proposed an iterative blind deconvolution (IBD) method, which they applied to the restoration of turbulence-degraded images. The method alternates between estimating the OTF and estimating the image. Image-domain constraints of nonnegativity of the OTF and the image are used during the iterations. Although promising results have been obtained in simulations, the uniqueness and convergence properties of the IBD algorithm remain unclear. Moreover, noise amplification is often observed. To make IBD more robust to noise, the Lucy-Richardson algorithm can be inserted into the IBD framework. We refer to this combination as the adaptive Lucy-Richardson maximum likelihood (LR) algorithm. In this work, this IBD algorithm is compared with our Principal Components Analysis (PCA) method for restoring atmospheric turbulence degraded images. Deconvolution of an unknown finite impulse response (FIR) blur can be viewed as a blind source separation/extraction problem. With this point of view, the source (the original image) is repeatedly shifted, then linearly mixed. The mixing matrix corresponds to the blurring PSF, which is unknown. From the mixture, an optimal estimation of the true image needs to be computed. The first principal component can be used as such an estimate since it has maximum variance and it contributes most to the variance of the observed dataset. Intuitively, blurring is a smoothing process in which the high-frequency components are removed and the variance is reduced. Thus, deblurring should be able to boost highfrequency components to some extent to restore the image. Variance can be viewed as a measurement of high-frequency components in an image. Previously, we developed a constrained variance maximization method for blind image deconvolution [7]. That variance maximization is equivalent to a PCA. In this letter, for the first time, we explicitly formulate blind deconvolution as a PCA problem in both multichannel and single-channel cases. In [7], a spatially varying blur was considered. Although it might appear to be helpful if the deblurring were performed locally so that each local region (row, column or block) could be processed independently of others, practically, blurs are not separable. Moreover, artifacts are created at block boundaries when the image is processed by a local deblurring process. These artifacts are sometimes misinterpreted as noise amplification. In this letter, we do not attempt to address spatially varying blur and only consider shift-invariant blurs and apply PCA to the entire image. This letter is organized as follows. In Section II, the PCA-based blind image deconvolution algorithms for both multichannel and single-channel cases are presented. Noise

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. X, NO. XX, XXXX 2006

robustness and computational complexity are also analyzed in this section. In Section III, the details of the comparative experiments of PCA and the adaptive Lucy-Richardson algorithm on both real and simulated degraded images are presented. A concluding remark is presented in Section IV. II. B LIND I MAGE D ECONVOLUTION WITH PCA A. Multichannel Blind Image Deconvolution using PCA To simplify notation, we use one-dimensional discrete convolution in the blurring model g = f ∗ h,

(2)

where f is the original image and h is the blurring filter. Additive noise is also generally present. For the moment, it is ignored, but the effect of additive noise is analyzed later on after the PCA algorithm is derived. The length of the blurring filter (M ) is assumed to be 3 to make notation simpler, although other values are used in the simulation. h can be written as a sum of delayed unit impulses that are scaled by the sample values: h(n) =

∞ X

h(k)δ(n − k)

of the data (the observed blurred imagery). Principal Components Analysis (PCA) [8], [9] is an orthogonal transformation technique that can be used for dimensionality reduction of a dataset. In the context of blind image deconvolution, the dataset is a set of observed blurred images and the goal is to reduce the dimensionality of the image dataset to 1. PCA is a linear transformation that transforms the data to a new orthogonal coordinate system. In the new coordinate system, the first coordinate is the first principal component, the eigenvector associated with the maximum eigenvalue. The eigenvalue indicates the importance of the eigenvector to the variance of the data. Given multiple blurred versions of the same true image, the PCA-based multichannel blind image deconvolution takes the following steps: 1) Form an ensemble from the blurred images. Reshape each N × N image into a column vector gi . The ensemble is denoted as {g1 , g2 , g3 , . . . , gM }, a set of M observations of g. g is a column vector with N 2 scalar random variable components. 2) Compute the average image of the ensemble.

(3) Ψ=

k=−∞

=

2

h(−1)δ(n + 1) + h(0)δ(n) + h(1)δ(n − 1)

Applying the shift property of convolution, we can rewrite Eq. 2 as: g(n) = f (n) ∗ (h(−1)δ(n + 1) + h(0)δ(n) + h(1)δ(n − 1)) = h(−1)f (n + 1) + h(0)f (n) + h(1)f (n − 1) = [h(1), h(0), h(−1)][f (n − 1), f (n), f (n + 1)]T Assuming that there are 3 observations of the blurred images, we then have g1 (n) = [h1 (1), h1 (0), h1 (−1)][f (n − 1), f (n), f (n + 1)]T g2 (n) = [h2 (1), h2 (0), h2 (−1)][f (n − 1), f (n), f (n + 1)]T (4) g3 (n) = [h3 (1), h3 (0), h2 (−1)][f (n − 1), f (n), f (n + 1)]T where T denotes the vector transpose, and g1 [n], g2 [n] g3 [n] are three observed blurred versions of f [n]. The blurring filters for the three blurred images are h1 [n], h2 [n] and h3 [n], respectively. The equations in Eq. 4 can be put together in matrix-vector form as:      g1 h1 (1) h1 (0) h1 (−1) f (n − 1)  g2  =  h2 (1) h2 (0) h2 (−1)   f (n)  (5) g3 h3 (1) h3 (0) h3 (−1) f (n + 1) or G = HF

(6)

Note that two-dimensional blurs can also be cast in the form of Eq. 6. The matrix-vector form suggests that convolution can be viewed as a mixing process where the original image is shifted and linearly combined to generate the observed blurred images. Thus, the deblurring problem is to extract or estimate the original image from the mixtures. The optimal linear estimate using second-order statistics is the first principal component, which has the largest contribution to the variance

M 1 X gi M i=1

(7)

This sample mean, Ψ, is an estimate for E(g), the mathematical expectation of g. 3) “Centering” the ensemble. We subtract the mean Ψ from each gi to get Φi . This step prepares for the computation of the covariance matrix as in the next step. Put {Φi } together to form a matrix A = [Φ1 , Φ2 , . . . , ΦM ]. A is N 2 × M . 4) Construct the empirical covariance matrix. The covariance matrix of g is defined as E((g − E(g))(g−E(g))T ), which is estimated by the empirical covariance matrix: C=

1 AAT . M

(8)

C is N 2 × N 2 . If the image is 256 × 256, then C is 65536 × 65536. It is computationally intractable to determine the eigenvectors and eigenvalues of a matrix of such a size. However, it has been found [9] that the computation can be dramatically reduced by first computing the eigendecomposition of AT A, which is M × M. 5) Calculate the eigenvector. Let vi and µi denote the ith eigenvalue and eigenvector of the matrix AT A, AT Aµi = µi vi

(9)

Premultiplying both sides by A, we have AAT Aµi = Aµi vi

(10)

from which we can see that Aµi is an eigenvector of 1 AAT , and thus of C. The scalar M does not change the eigenvector.

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. X, NO. XX, XXXX 2006

6) Compute the estimated image. The eigenvectors are often sorted by their eigenvalues in decreasing order. Among all the eigenvectors, the first one (Aµ1 ) associated with the maximum eigenvalue contributes most to the variance of the ensemble. Therefore, we use it as the estimate of the image f [i, j]. Since the mean image Ψ was removed before the covariance matrix was computed, it needs to be added back: ˆ f = Aµ1 + Ψ.

B. Single-channel Blind Image Deconvolution using PCA Ideally, multiple observations of the blurred images should be available when PCA is applied in blind image deconvolution. However, when there is only one blurred image, an ensemble may be created by shifting the blurred image. The length of the 1-D blurring filter is assumed to be 3, as before, for simpler notation. The PCA model for single-channel blind image deconvolution is

(11)

Aµ1 can be interpreted as the projection of A onto the first eigenvector of AT A, and ˆ f provides an estimate of the ideal image. Finally, the vector ˆ f is reshaped into the 2D matrix to get the restored image. In Eq. 11, there are two components in the restored image. Aµ1 can be interpreted as the high-frequency component and the low-frequency component is Ψ. µ1 can be viewed as a filter to boost high-frequency content. Now we are ready to analyze the robustness to noise, which is assumed to be white. Since Ψ is a sample mean, it will have reduced noise because of averaging. The white noise added to the blurred images is assumed to be uncorrelated with the image. The covariance matrix of the noisy blurred images is made up of 2 parts: C = UΛUT + σ 2 I = UΛUT + σ 2 UIUT = U(Λ + σ 2 I)UT

3

(12)

in which U is the transformation matrix of the noise-free part of the data and Λ is the diagonalized matrix. It can be seen from Eq. 12 that the noise variance σ 2 is added to eigenvalues in Λ and the eigenvectors in U remain the same. Since σ does not change the order of the eigenvalues, the first principal eigenvector remains the same. This explains the noise robustness of the PCA algorithm. In the experiment section, noise robustness experiments results are reported. For the multichannel PCA restoration algorithm, the number of images in the ensemble (M ) is the only parameter. Generally, it is preferable for M to be large since the covariance matrix can be better estimated with more samples. However, a larger M can lead to artifacts in the PCA single-channel deconvolution as shown in subsection B in the experiments section. Only the first eigenvector of a rather small matrix (M × M ) needs to be computed. The singular value decomposition (SVD) of the matrix can be computed by the standard LAPACK [10] subroutine dgesvd with complexity O(M 3 ) [11]. By way of comparison, the computational complexity of the adaptive Lucy-Richardson algorithm is dominated by the fast Fourier transform (FFT), which is O(N 2 log N ) for a N × N image. Since the number of images (M ) is usually much smaller than the number of pixels in the image (M ¿ N 2 ), the computation time of PCA is generally much lower than that of the adaptive Lucy-Richardson algorithm. For convenient notation, the size of the image was assumed to be square (N × N ). However, there is no restriction regarding the size of the image in the PCA restoration algorithm.

G = HF where G

 =  

H

=  

F

  =   

 g(n − 1) g(n)  g(n + 1) h(1) h(0) h(−1) 0 h(1) h(0) 0 0 h(1)  f (n − 2) f (n − 1)   f (n)   f (n + 1)  f (n + 2)

(13)

(14) 0 h(−1) h(0)

 0  (15) 0 h(−1)

(16)

The following steps describe the PCA algorithm to restore an image blurred with an unknown spatially invariant PSF. When applying the PCA algorithm to a single image, most of the steps remain the same as in the multichannel algorithm except for the first few steps. 1) Create an ensemble. An ensemble is first created by shifting the image. Assuming that the PCA restoration filter PSF support is (2R+1)×(2S+1), there will be (2R+1)(2S+1) images denoted as {g[i−k, j −l]|(k, l) ∈ {−R, . . . , 0, . . . , R}× {−S, . . . , 0, . . . , S}}. 2) Compute the average of the ensemble. Ψ=

R S X X 1 g[i − k, j − l] (17) (2R + 1)(2S + 1) k=−R l=−S

3) The remaining steps are the same as in the multichannel case. III. E XPERIMENTS The PCA blind deconvolution algorithm was implemented and tested on a variety of real and simulated turbulence blurred images. For the real blurred images and video clips, there was no information about the atmospheric turbulence condition (λ of the turbulence etc) or the imaging system (wavelength, angle, exposure etc) used when they were acquired. For computer simulation of the turbulence blur, the long-exposure optical transfer function in Eq. 1 is used. The number of images in the ensemble (M ) is the only parameter for the PCA restoration algorithm. The adaptive Lucy-Richardson maximum likelihood algorithm (LR) is used to restore the same images for comparison purposes. The LR algorithm requires an initial estimate of the blur PSF. We have observed that the size of the support of the initial PSF is more important to the

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. X, NO. XX, XXXX 2006

4

are used. As the PSF support increases, the edges in the PCA restored image become stronger; this effect is less obvious for the adaptive Lucy-Richardson algorithm. The Lucy-Richardson algorithm uses the discrete Fourier transform (DFT), which assumes that the image is periodic. This assumption results in boundary-related ringing in the deblurred images. This ringing effect is observed in Fig. 2(b) and Fig. 2(d). PCA based restoration does not introduce DFT-related ringing effects. C. Noise Robustness Experiments on Simulated Atmospheric Turbulence Degraded Images

(a)

The moon surface image in the Miscellaneous volume from the USC-SIPI database [12] was used in this test. The image was blurred by the OTF in Eq. 1 (λ = 0.003). Gaussian random noise was added at a level of 10 dB BSNR. The blurred signal-to-noise ratio (BSNR) is defined as BSN R = 10 log

(b) Fig. 1. (a) The 8th frame of the atmospheric turbulence degraded video. (b) The multichannel PCA restored image.

ultimate success of the restoration algorithm than the actual values selected. Because of this fact, we used an averaging filter with the appropriate support as the initial estimate. The performance of the LR algorithm is also very sensitive to the number of iterations. A. Multichannel Deconvolution on Real Turbulence Degraded Images A number of consecutive frames are taken from an atmospheric turbulence degraded video clip. The information about the turbulence and the imaging device is unknown. We used the first 9 frames in one of the video clips. Fig. 1(a) shows the 8th frame in the video clip. Fig. 1(b) shows the image restored by the multichannel PCA blind deconvolution method. The image looks sharper and fine details can be observed. B. Single-channel Deconvolution on Real Turbulence Degraded Images An example of the restoration on a real turbulently blurred image is shown in Fig. 2. Both methods boost the highfrequency contents of the image. Two different PSF supports

σb2 , σn2

(18)

where σb2 and σn2 are the variances of the blurred image and noise respectively. Fig. 3 shows the comparative result on this simulated degraded image. Noise amplification is obvious in the adaptive Lucy-Richardson algorithm as shown in Fig. 3(c). Much of the noise was removed in the PCA result as seen in Fig. 3(d). On the other hand, some artifacts are observed in the PCA restored image where the edges appear unnatural. As analyzed previously, the PCA method has lower computational complexity when compared with the adaptive LucyRichardson algorithm. The programs were run in Matlab on a notebook computer. The CPU is a 1400MHz Pentium M processor, with 768MB of RAM. On average, it takes 0.22 seconds to deblur a 240 by 240 pixel image; while on average the Lucy-Richardson algorithm takes 4.4 seconds (a factor of 20 longer) using 10 iterations. IV. C ONCLUSION We have introduced a blind image deconvolution method based on PCA and applied it to the task of restoring atmospheric turbulence degraded images. Some of the main features of the PCA based blind deconvolution algorithm are the following: 1) PCA is fast. 2) PCA is robust to noise. 3) PCA is primarily a multichannel blind image deconvolution method, but it can also be applied in the singlechannel case. So far, no knowledge about atmospheric turbulence is assumed and used when PCA is applied to restore the image. However, information about the turbulence can be used in the PCA method to select a proper deblurring filter PSF support to achieve better result as the experiments show that an improper PSF support can introduce artifacts in the restoration result. One of well-known limitations of variance is that it is phase blind. That is, variance does not reflect information about phase; it is only determined by the magnitude of an image in the frequency domain. It is expected that higher

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. X, NO. XX, XXXX 2006

5

(a)

(b)

(c)

(d)

(a)

(b)

(c)

(d)

(e)

Fig. 2. (a) An atmospheric turbulence degraded image. (b) The adaptive Lucy-Richardson algorithm restored image, 10 iterations, 3 × 3 PSF support. (c) The image restored with PCA, 3 × 3 PSF support. (d) The adaptive LucyRichardson algorithm restored image, 8 iterations, 5 × 5 PSF support. (e) The image restored with PCA, 5 × 5 PSF support.

order statistics might yield better restoration result since higher order statistics such as kurtosis can reflect information about phase. However, variance, or second order central moment, has obvious advantage since it leads to an analytic closedform solution while higher order statistics does not. V. ACKNOWLEDGMENTS We greatly appreciate the helpful comments from the anonymous reviewers. R EFERENCES [1] M. C. Roggemann and B. Welsh, Imaging Through Turbulence, CRC Press, 1996.

Fig. 3. (a) The original moon surface image from the USC-SIPI database. (b) The degraded image (λ = 0.003), noise is added at the level of BSNR = 10 dB. (c) The adaptive Lucy-Richardson algorithm restored image, 10 iterations, 5 × 5 PSF support. (d) The image restored with PCA, 5 × 5 PSF support.

[2] R. E. Hufnagel and N. R. Stanley, “Modulation transfer function associated with image transmission through turbulence media,” J. Opt. Soc. Amer. A., vol. 54, pp. 52–61, 1964. [3] D. Kundur and D. Hatzinakos, “Blind image deconvolution revisited,” IEEE Signal Process. Mag., vol. 13, no. 6, pp. 61–63, Nov 1996. [4] L. Lucy, “An iterative technique for the rectification of observed distribution,” Astron. J., vol. 79, no. 6, pp. 745–754, 1974. [5] W. H. Richardson, “Bayesian-based iterative method of image restoration,” J. Opt. Soc. Amer. A., vol. 62, pp. 55–59, 1972. [6] G. R. Ayers and J. C. Dainty, “Iterative blind deconvolution method and its applications,” Opt. Lett., vol. 13, pp. 547–549, July 1988. [7] D. Li, S. Simske, and R. M. Mersereau, “Blind image deconvolution using constrained variance maximization,” in Proc. Asilomar Conf. on Signals, Syst. and Computers, 2004. [8] H. Anton, Elementary Linear Algebra, John Wiley & Sons, 1987. [9] M. Turk and A. Pentland, “Eigenfaces for recognition,” J. Cognitive Neuroscience, vol. 3, no. 1, pp. 71–86, 1991. [10] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, and D. Sorenson, Lapack User’s Guide, SIAM, 1995. [11] G. H. Golub and C. F. Van Loan, Matrix Computation, John Hopkins University Press, 1996, Third Edition. [12] University of Southern California - Signal & Image Processing Institute The USC-SIPI Image Database, http://sipi.usc.edu/services/database/Database.html.