BLUR IDENTIFICATION BASED ON KURTOSIS MINIMIZATION Dalong Li, Russell M. Mersereau
Steven Simske
School of Electrical and Computer Engineering Georgia Institute of Technology Atlanta, GA 30332 {dalong,rmm}@ece.gatech.edu
Imaging Systems Laboratory Hewlett-Packard Laboratories Fort Collins, CO 80528
[email protected] ABSTRACT In this paper, we describe an algorithm for identifying a parametrically described blur based on kurtosis minimization. Using different choices for the parameters of the blur, the noisy blurred image is restored using Wiener filter. We use the kurtosis as a measurement of the quality of the restored image. From the set of the candidate deblurred images, the one with the minimum kurtosis is selected. The proposed technique is tested in a simulated experiment on a variety of blurs including atmospheric turbulence blurs, Gaussian blurs, and out-of-focus blurs. The proposed approach is also tested on real blurred images. Moreover, we test the performance when a wrong blur model is given. Our experiments show that the kurtosis minimization measurements match well with methods that maximize PSNR. 1. INTRODUCTION The primary objective of image restoration is to recover the visual information from a degraded image. It has important applications in photographic deblurring, remote sensing, and medical imaging. The blur may be caused by degradations in the imaging process, such as defocusing, atmospheric turbulence, object motion, or diffraction. In many imaging applications, the observed image g(x, y) can be approximated as the sum of a two-dimensional convolution of the true image,f (x, y), with a linear shift-invariant blur, (also known as the point-spread-function (PSF)), h(x, y), and additive noise. Mathematically, ~ +n g = f ∗ h(θ)
(1)
θ~ ∈ Ω in which * denotes the two-dimensional linear convolution operation, n denotes the additive noise, and θ~ is a parametrization of the blur. The problem of recovering the true image f (x, y) from the degraded image g(x, y) is called image deconvolution or image restoration. Classical restoration methods require complete knowledge of the blur PSF
prior to restoration. However, it is often impossible, or, in some cases, impractical to determine the blur a priori. These could be because various practical constraints such as the difficulty of characterizing atmospheric turbulence in aerial imaging or the potential health hazard of using a stronger incident beam to improve the image quality in X-ray imaging. In these situations, blind image deconvolutions are essential for recovering visual information from the degraded images. Many blind restoration algorithms have been proposed in the past [1]. All these algorithms make certain assumptions to make the problem tractable. Principle component analysis (PCA) [2] assumes that the pixel values in the true image are spatially uncorrelated. Training images are required to use support vector regression [3], which is then able to perform deconvolution of similar, but unknown blurs. Variational approximation [4] is proposed for Bayesian blind image deconvolution which assumes that the PSF is partially known. Then the probabilistic law relating the observations and the quantities to be estimated can be formulated. A commonly used model for the deblurring problem is an auto-regressive moving average (ARMA) model, in which the image is modeled as an autoregressive process and the blur is modeled as moving average process. Many algorithms use this ARMA formulation, including maximum likelihood [5] and generalized cross validations [6]. Let us use fˆ(θ~h ) to denote the deblurred image with the hypothetical parameter vector θ~h . Given the set {Ψ : fˆ(θ~h |θ~h ∈ Ω)} the question that we are interested in is: which restoration should we choose. The problem is easy if f is available and you can use PSNR as a measurement. The one with maximum PSNR is the best one, since, by this criterion it is the closet fit. However, f is not available. We have discovered that the kurtosis, which is a measure of how outlier-prone a distribution is, is a useful criterion for selecting the best deblurred image. The reminder of the paper is organized as follows. In Section 2, the kurtosis is reviewed and related to our blur identification problem . In Section 3, the implementations
and the experiments results of our approach are reported. Some brief conclusions is contained in Section 4.
The Wiener filter is given by: G(u, v) =
2. KURTOSIS BASED BLUR IDENTIFICATION
H ∗ (u, v) |H(u, v)|2 + Pn (u, v)/Pf (u, v)
(3)
Where Pf and Pn are the power spectral of the signal and noise, respectively. The kurtosis minimization rule can be written as
8000 histogram of the original image (k=2.0915) histogram of the 3 by 3 averaging filtered image (k=2.1607) histogram of the 7 by 7 averaging filtered image (k=2.1991)
7000
~ θ~k = min k(fˆ(θ)) ~ θ∈Ω
number of the pixel
6000
The benchmark rule is the PSNR maximization rule,
5000
~ R = max P SN R(fˆ(θ)) ~ θP SN θ∈Ω
4000
(5)
PSNR is defined as:
3000
PM PN 2 i=1 j=1 255 ˆ P SN R(f ) = 10 log10 PM PN 2 ˆ i=1 j=1 (f (i, j) − f (i, j)) (6) where M and N are the height and width of the image.
2000 1000 0
(4)
0
10
20 30 bin of the histogram
40
50
Fig. 1. Histogram of the original and the blurred cameraman image.
3. EXPERIMENTS 3.1. Evaluation of the blur identification
The kurtosis of a random variable is defined as the normalized fourth central moment k=
E((x − µ)4 ) σ4
(2)
where µ is the mean of x, σ is its standard deviation , and E(x) represents the expected value operation. The kurtosis measures the peakedness of a distribution. The normal distribution (k = 3) has a moderate tail and is called mesokurtic. A distribution with a small tail has a small kurtosis (k < 3) and is called platykurtic, and one with a long tail (k > 3) is called leptokurtic. For a platykurtic distribution, the larger the kurtosis, the smoother the data. Fig.1. shows the histograms of the original cameraman image and those of two smoothed versions. As the smoothing increases, the corresponding kurtosis increases as the data becomes more Gaussian. Given a noisy blurred image g with a known functional form for the PSF, we seek the best estimate of the blur pa~ The parameter is searched within a reasonable rameter θ. space Ω. At each step in the search loop, the image is deblurred using a Wiener filter G(u, v) and the kurtosis of the ~ is computed and stored. Then the dedeblurred image fˆ(θ) blurred image with the smallest kurtosis is chosen as the restored image and the corresponding parameter is regarded as the identified blurring parameter θ~k .
From the degradation model (1), it is clear that both h~(θ) and n contribute to the degradation. The restoration depends on both of them. Even if θ~ is perfectly identified, it is not guaranteed to obtain the best deblurred image if the noise is not well estimated. Therefore, it makes sense to compare the PSNR of the deblurred image, rather than to ~ compare how close estimated θ~ is to the real θ. To measure the accuracy of the blur identification, we measure the PSNR of the corresponding deblurred image ~ Since PSNR maximization rule is the benchmark, we fˆ(θ). ~ R )). compare P SN R(fˆ(θ~k )) with P SN R(fˆ(θP SN ~ R )) − P SN R(fˆ(θ~k )))| (7) δP SN R = |(P SN R(fˆ(θP SN 3.2. Atmospheric turbulence blur identification The optical transfer function (OTF) of the blur is H(u, v) = e−λ(u
2
+v 2 )5/6
(8)
λ controls the severity of the blur. When λ is zero, there is no blur; as λ increases, the blur becomes stronger. λ is the only parameter to be determined. The cameraman image is blurred with the turbulence OTF model and Gaussian random noise (σ = 0.0025) is added to the blurred image. λ varies from 0.0012 to 0.004 in steps of 0.0002. 15 different blured images are made. The parameters for the turbulence
blur in these 15 blurred images are to be estimated. The parameter search space Ωλ = {λ : 0.001 + i ∗ 0.0002|i = 0, 1, 2, . . . , 15}. The accuracy of the 15 identifications are shown in Fig.2. The deblurred image selected based on the kurtosis minimization rule are very close to those selected by the PSNR maximization rule. Fig.3 shows the result on a real turbulence blurred image for which the model parameters are unknown. 0.8 0.7 0.6
0.4
δ
PSNR
(db)
0.5
0.3
hsize 2 3 4 5 6 7 σ hsize 2 3 4 5 6 7 σ
Table 1. PSNR and Kurtosis of fˆ(hsize , σ) P SN R 16.974 16.974 16.974 16.974 16.974 16.192 16.105 16.078 16.064 16.06 15.105 15.842 16.27 16.501 16.671 17.163 22.964 25.09 24.273 23.248 15.6 16.204 14.851 13.787 13.224 13.071 14.858 12.107 10.891 10.29 1.5 2 2.5 3 3.5 Kurtosis 2.2452 2.2452 2.2452 2.2452 2.2452 2.2476 2.2486 2.2497 2.2525 2.2545 2.2323 2.2235 2.2163 2.2139 2.2128 2.1852 2.1518 2.1417 2.1487 2.1554 2.2062 2.1642 2.194 2.2115 2.2055 2.1948 2.1768 2.1929 2.1897 2.1701 1.5 2 2.5 3 3.5
0.2 0.1 0
1
1.5 2 2.5 3 3.5 λ: the real atmospheric turbulence blur parameter
4 −3
x 10
Fig. 2. Identification of atmospheric turbulence blur with different λ
The search space for hsize is Ωhsize = {2, 3, 4, 5, 6, 7} and that for σ is Ωσ = {1.5, 2, 2.5, 3, 3.5}. Table.1 shows that the minimum value of the kurtosis is achieved for (5, 2.5), where the PSNR is maximum. 3.4. Out-of-focus blur identification Out-of-focus blur is usually modeled as a circular averaging filter (pillbox) within the square matrix of side 2 ∗ γ + 1, where γ is the radius of the blur. Two examples are shown here. The true radii are 2 and 3 respectively. The cameraman image is blurred with the pillbox filter. Gaussian noise (σ = 0.0025) is added to the blurred image. The search space for γ is Ωγ = {γ : 1 + i ∗ 0.1|i = 0, 1, 2, . . . , 30}. As seen in Fig.4, although the estimated radius (3.1) is different from the real radius (3) the PSNR is nonetheless higher. The reason is noise estimation error. In example (b), the kurtosis correctly identifies that the radius is 2 but it does not correspond to the maximum PSNR.
(a)
(b) 3.5. Blur identification with incorrect blur model
Fig. 3. (a) Real atmospheric turbulence blurred image. (b) Restored image, estimated λ = 0.0024
3.3. Gaussian Blur identification We use a rotationally symmetric Gaussian lowpass filter of size hsize with standard deviation σ (positive) to simulate a Gaussian blur. There are 2 parameters to be determined. hsize and σ. The cameraman image is blurred by the Gaussian filter (hsize = 5, σ = 2.5). Gaussian noise (σ = 0.0025) is added to the blurred image.
In this experiment, we blur the cameraman image with turbulence (λ = 0.0025) and add white noise (σ = 0.002). Then we use the kurtosis to identify the blur, assuming that the blur is an averaging filter. The search space of hsize Ωhsize = {3, 4, 5, 6}. The estimated hsize of the blur is 5 and the PSNR of the restored image is 20.78. The PSNR of the restored image using the correct model and parameter is 22.86. This mismatch is not surprising since the model is not correct, but the image is still improved as shown in Fig.5. This illustrates that the technique is robust to model errors.
4. CONCLUSION 1.2
We propose a higher order statistical based criterion (kurtosis) to select the best deblurred image from a set of candidate deblurred images. Each candidate deblurred image is made using a Wiener filter with a hypothesized parameter of the blur taken from a parameter search space. The proposed method is tested on a variety of different blurs and real blurred images. The kurtosis minimization produces very similar results to those selected by the PSNR maximization criterion. We also observe that when evaluating blur identification, it makes more sense to compare PSNR of the deblurred image to compare the closeness of the estimated parameters to the real parameters, since the effect of noise must be considered.
1 0.8
0.6 0.4 0.2 0 −0.2 1
[1] D. Kundur and D. Hatzinakos, “Blind Image Deconvolution Revisited,” IEEE Signal Processing Magazine, vol. 13, no. 6, pp. 61–63, Nov 1996.
1.2
[2] Dalong Li, Steven Simske and Russell M. Mersereau, “Blind Image Deconvolution using Constrained Variance Maximization,” in the Proc. of the 38th Asilomar Conf. on Signals, Systems and Computers, Pacific Grove, CA, USA, 2004.
0.8
4
Normalized Kurtosis Normalized PSNR
1
0.6 0.4 0.2 0 −0.2 1
[4] A.C. Likas and N.P. Galatsanos, “A Variational Approach for Bayesian Blind Image Deconvolution,” IEEE Transactions on Signal Processing, vol. 52, pp. 2222–2233, Aug. 2004.
1.5 2 2.5 3 3.5 r: hypothetical radius of the circular averaging filter
(γ = 3, estimated radius: 3.1)
5. REFERENCES
[3] Dalong Li, Russell M. Mersereau and Steven Simske, “Blind Image Deconvolution Using Support Vector Regression,” in the Proc. of the IEEE ICASSP, Philadelphia, PA, USA, 2005.
Normalized Kurtosis Normalized PSNR
1.5 2 2.5 3 3.5 r: hypothetical radius of the circular averaging filter
4
(γ = 2, estimated radius: 2) Fig. 4. Out-of-focus blur identification
[5] R.L. Lagendijk, A.M. Tekalp and J. Biemond, “Maximum Likelihood Image and Blur Identification: a Unifying Approach,” Optical Engineering, vol. 29, no. 5, pp. 422–435, May 1990. [6] S.J. Reeves and R.M. Mersereau, “Blur Identification by the Method of Generalized Cross-validation,” IEEE Trans Image Processing, vol. 1, no. 3, pp. 301–311, July 1992.
(a)
(b)
Fig. 5. (a) Atmospheric turbulence blurred image. (b) Restored image using averaging filter blur model