# Deconvolution of Defocused Image with Multivariate Local Polynomial Regression and Iterative Wiener Filtering in DWT Domain

Authors: Liyun Su; Fenglan Li
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

A novel semiblind defocused image deconvolution technique is proposed, which is based on multivariate local polynomial regression (MLPR) and iterative Wiener filtering (IWF). In this technique, firstly a multivariate local polynomial regression model is trained in wavelet domain to estimate defocus parameter. After
obtaining the point spread function (PSF) parameter, iterative wiener filter is adopted to complete the restoration. We experimentally illustrate its performance on simulated data and real blurred image. Results show that the proposed PSF parameter estimation technique and the image restoration method are effective.


---

## Body

## 1. Introduction

Shooting a real world image with a camera through an optical device gives a 2D image where at least some parts are affected by a blur and noise. Images can be blurred by atmospheric turbulence, relative motion between sensors and objects, longer exposures, and so on, but the exact cause of blurring may be unknown. Restoration of blurred noisy images [1–3] is one of the main topics in many processing. The literatures [4–6] have given good methods to improve image qualities. The purpose of image restoration is to reconstruct an unobservable true image from a degraded observation. An observed image can be written, ignoring additive noise, as the two-dimensional (2D) convolution of the true image with a linear space-invariant (LSI) blur, known as the PSF. Restoration in the case of known blur, assuming the linear degradation model, is called linear image restoration and it has been presented extensively in the last three decades giving rise to a variety of solutions [7–10]. In many practical situations, however, the blur is unknown. Hence, both blur identification and image restoration must be performed from the degraded image. Restoration in the case of unknown blur is called blind image restoration [11–13]. Existing blind restoration methods can be categorized into two main groups: (i) those which estimate the PSF a priori independent of the true image so as to use it later with one of the linear image restoration methods, such as zero sheet separation, generalized cross validation, and maximum likelihood and expectation maximization based on the ARMA image model [14–16], and (ii) those which estimate the PSF and the true image simultaneously, such as nonnegative sand support constraints recursive inverse filtering, maximum likelihood and conjugate gradient minimization, and simulated annealing [17, 18]. Algorithms belonging to the first class are computationally simple, but they are limited to situations in which the PSF has a special form, and the true image has certain features. Algorithms belonging to the second class, which are computationally more complex, must be used for more general situations. In this paper, a kind of semi-blind image restoration algorithm is proposed in case of known the blur type (defocused blurring).In general, discrete model for a linear degradation caused by blurring can be given by the following equation:(1.1)y(i,j)=h(i,j)*f(i,j)+n(i,j),
where * indicates two-dimensional convolution, f(i,j) represents in original image, y(i,j) is the degraded image, h(i,j) represents the two-dimensional PSF, and n(i,j) is the additive noise. In this paper, we deal only with additive Gaussian noise, as it effectively models the noise in many different imaging scenarios. The difficulty in solving the restoration problem with a spatially varying blur commonly motivates the use of a stationary model for the blur. This leads to the following expression for the degradation system:
(1.2)y(i,j)=h(i,j)*f(i,j)+n(i,j)=∑k=1M∑l=1Nh(i-k,j-l)f(k,l)+n(i,j).The use of linear techniques for solving the restoration problem is facilitated by using space-invariant model. Models that utilize space-variant degradations are also common, but lead to more complex solutions. As for defocus blur,PSF is modeled as a uniform intensity distribution within a circular disk,(1.3)h(i,j)={1πR2,ifi2+j2≤R,0,otherwise,
where disk radius R is the only unknown parameter for this type of blur.Many existing image restoration algorithms assume that thePSF is known, but in practical it is not always the case. The restoration without knowing of the PSF is called blind image restoration. Fourier methods can be used to estimate the defocus parameter R through calculating a ratio of power of high frequencies portion to that of low frequencies portion. However, a main drawback of the method is its bad noise immunity. To solve this problem, [19] proposed a novel algorithm to overcome this shortcoming, the RBF neural network is applied to fit R. This scheme has good fitting, but bad prediction. To avoid the weak generalization ability, an more efficient method for estimating parameter R is proposed in this paper. The prediction ability of these two methods is compared with the trained five images. The steps of the presented algorithm in this paper is as follows: Firstly we construct feature vectors of several blurred images with known defocus radius R in wavelet domain, then a multivariate local polynomial estimation model [20–22] is trained using the vectors as inputs and defocus parameters as outputs. After the model is trained, the new defocused images are applied to the trained model for predicting the parameter R. For a semi-blind defocused image, R can be estimated through calculating the feature vectors and using it as input of the trained model. With known radius R, many traditional algorithm could be applied to restore the degraded image. In this paper, iterative Wiener filtering [23] is adopted to image restoration.


## 2. Relationship between Wavelet Coefficients andR

The wavelet transform provides a powerful and versatile framework for image processing. It is widely used in the fields of image denoising, compression, fusion, image restoration [24–26], and so forth.The two-dimensional discrete wavelet transform (DWT) [27, 28] hierarchically decompose an input image into a series of successively lower resolution images and their associated detail images. DWT is implemented by a set of filters, which are convolved with the image rows and columns. An image is convolved with low-pass and high-pass filters and the odd samples of the filtered outputs are discarded resulting in downsampling the image by a factor of 2. The l level wavelet decomposition of an image I results in an approximation image Xl and three detail images Hl,Vl, and Dl in horizontal, vertical, and diagonal directions respectively. Decomposition into l levels of an original image results in a down sampled image of resolution 2l with respect to the image as well as detail images.When an image is defocused, edged in it are smoothed and widened. The amount of high frequency band decreased, and that corresponding to low frequency band increases.In order to denote the relationship between wavelet coefficients and defocused radiusR, we define five variables named v1,v2,v3,v4, and v5 as follows:
(2.1)v1=|V2|s|H2|s,v2=|H2|s|X2|s,v3=|H1|snum{H1},v4=|H2|snum{H2},v5=|D1|snum{D1},
where |·|s represents the summation of all coefficients' absolute value, num{·} is total number of coefficients.An original image is blurred artificially by a uniform defocusPSF with R whose value rang from 1 to 20. The relationship between v1,v2,v3,v4,v5 and R is shown in Figure 1, where the curves are normalized in [0,1] interval. When R increases, v2,v3,v4 and v5 decrease monotonously.Relationship Betweenv1-5 and R.
(a)(b)In order to estimate defocus parameterR, only known the roughly similar relationship is not enough. As shown in Figure 2, every image has monotonous curve between v2,v5 and R, but they are not superposition. For a degraded unknown PSF image, R cannot be calculated because the curve of the given image is not known. For example, if v2 of image “rice” has been calculated, and then we estimate R according curve if “ic” in Figure 2, wrong results are obtained obviously. To solve this problem, one of the methods is to choose neural networks. Computational artificial neural networks are known to have the capability for performing complex mappings between input and output data, but neural network method has bad generalization ability. Here we propose a multivariate local polynomial regression model to estimate R. The variables v1-5 are chosen to train the multivariate local polynomial estimation model. Prediction Comparison is made to verify the advantages of multivariate local polynomial fitting.Curvesv2 and v5 of Different Images.
(a)(b)


## 3. Training Multivariate Local Polynomial Estimation Model

Multivariate local polynomial fitting is an attractive method both from theoretical and practical point of view.  Multivariate local polynomial method has a small mean squared error compared with the Nadaraya-Watson estimator which leads to an undesirable form of the bias and the Gasser-Muller estimator which has to pay a price in variance when dealing with a random design model. Multivariate local polynomial fitting also has other advantages. The method adapts to various types of designs such as random and fixed designs, highly clustered and nearly uniform designs. Furthermore, there is an absence of boundary effects: the bias at the boundary stays automatically of the same order as the interior, without use of specific boundary kernels. The local polynomial approximation approach is appealing on general scientific grounds: the least squares principle to be applied opens the way to a wealth of statistical knowledge and thus easy generalizations. In this Section, we briefly outline and review the idea of the extension of multivariate local polynomial fitting [20–22] to the parameter R of defocused PSF.


### 3.1. Multivariate Kernel Function

To localize data in them-dimension, we need a multikernel function. Generally speaking, a multivariate kernel function refers to an  m-variate function satisfying
(3.1)∫-∞+∞⋯∫-∞+∞K(x̲)dx̲=1.Here and hereafter, we use∫ to indicate multivariate integration over the m-dimensional Euclidean space.There are two common methods for constructing multivariate kernel functions. For a univariate kernelk(x), the product kernel is given by
(3.2)K(x̲)=∏i=1mk(xi),
and the spherically symmetric kernel is defined as
(3.3)K(x̲)=cK,mK(∥x̲∥),
where cK,m={∫K(∥x̲∥)dx̲}-1 is a normalization constant and ∥x̲∥=(x12+x22+⋯+xm2)-1/2. Popular choices of K(x̲) include the standard d-variate normal density
(3.4)K(x̲)=(2π)-m/2exp(-∥x̲∥22),
and the spherical Epanechnikov kernel(3.5)K(x̲)={d(d+2)Γ(m/2)4πm/2}(1-∥x̲∥2)+.
The latter is the optimal kernel, according to Fan and Gijbels [21] and Su [22].The localization in multivariate nonparametric regression is frequently carried out by the kernel weighting. LetH̲ be a symmetric positive-definite matrix called a bandwidth matrix. The localization scheme at a point x̲ assigns the weight
(3.6)KH̲(Xi̲-x̲),withKH̲(x̲)=|H̲|-1K(H̲-1x̲),
where |H̲| is the determinant of the matrix H̲. The bandwidth matrix is introduced to accommodate the dependent structure in the independent variables. For practical problems, the bandwidth matrix H̲ is taken to be a diagonal matrix. The different independent variables will be accommodated into different scales. For simplification, the bandwidth matrix is designed into H̲=hIm̲ (Im̲ denoting the identity matrix of order m).


### 3.2. Multivariate Predictor with Local Polynomial Fitting

Suppose that the input vector isV̲=(v1,v2,v3,v4,v5), the model is fitted by the function(3.7)R=f(V̲).Our purpose is to obtain the estimationR̂=f̂(V̲) of function f, this paper, we use the dth order multivariate local polynomial f(V̲) to predict the defocused parameter RT value based on the point VT̲ of the test image. The polynomial function can be described as(3.8)f(V̲)≈∑0≤|j̲|≤d1j̲!D(j̲)fi(VT̲)(V̲-VT̲)j̲=∑0≤|j̲|≤db̲j̲(VT̲)(V̲-VT̲)j̲,
where(3.9)m=5,j̲=(j1,j2,…,jm),j̲!=j1!j2!⋯jm!,|j̲|=∑l=1m‍jl,∑0≤|j̲|≤d=∑|j̲|=0d(∑j1=0|j̲|∑j2=0|j̲|⋯∑jm=0|j̲|)|j̲|=j1+j2+⋯+jm,V̲j̲=v1j1v2j2⋯vmjm,D(j̲)fi(VT̲)=∂|j̲|fi(V̲)∂v1j1∂v2j2⋯∂vmjm|V̲=VT̲,b̲j̲(VT̲)=1j̲!D(j̲)fi(VT̲).In the multivariate prediction method,VTa̲(a=1,2,…,A) denotes the trained image feature vectors. Using A pairs of (VTa̲,Ra), for which the values are already known, the coefficients of fi is determined by minimizing
(3.10)∑a=1A[Ra-∑0≤|j̲|≤db̲j̲(VT̲)(VTa̲-VT̲)j̲]2·KH̲(VTa̲-VT̲).For the weighted least squared problem, a matrix form can be described by(3.11)W̲1/2·Y̲=W̲1/2·X̲·B̲+ɛ̲,
where(3.12)Y̲=(y1,y2,…,yA)T,ya=Ra,B̲=(b̲0(VT̲),b̲1(VT̲),…,b̲d(VT̲))T,W̲=diag{KH̲(VT1̲-VT̲),KH̲(VT2̲-VT̲),…,KH̲(VTA̲-VT̲)},
and X̲ is the A×S(S=∑0≤|j̲|≤d|j̲|/j̲!),(3.13)X̲=(1(VT1̲-VT̲)1⋯(VT1̲-VT̲)d1(VT2̲-VT̲)1⋯(VT2̲-VT̲)d⋮⋮⋱⋮1(VTA̲-VT̲)1⋯(VTA̲-VT̲)d).We then have the least squared solution with multivariate local polynomial fitting as follows:(3.14)B̲̂=(W̲1/2X̲)†Y̲
or, when X̲TW̲X̲ is inverse, the estimation can be written by(3.15)B̲̂=(X̲TW̲X̲)-1X̲TW̲Y̲,
then we can get the estimation R̂T=f̂(VT̲)(3.16)R̂T=f̂(VT̲)=E1̲(X̲TW̲X̲)-1X̲TW̲Y̲,
where E1̲=(1,0,0,…,0)1×S.Computing theB̲̂ will suffer from large computational cost. we can use the recursive least squared method to reduce the computation complexity, and it is very powerful especially in the real time prediction problems. There are several important issues about the bandwidth, the order of multivariate local polynomial function and the kernel function which have to be discussed. The three problems will be presented in Section 3.3.


### 3.3. Parameters Selections

For the multivariate local polynomial predictor, there are three important problems which have significant influence to the prediction accuracy and computational complexity. First of all, there is the choice of the bandwidth matrix, which plays a rather crucial role. The bandwidth matrixH̲ is taken to be a diagonal matrix. For simplification, the bandwidth matrix is designed into H̲=hIm̲. So the most important thing is to find the bandwidth h. A too big bandwidth underparameterizes the regression function, causing a large modeling bias, while a too small bandwidth overparameterizes the unknown function and results in noisy estimates. In theory, there exists an optimal bandwidth hopt in the meaning of mean squared error, such that
(3.17)hopt=argminh∫(f(x̲)-f̂(x̲))2dx̲.But the optimal bandwidth cannot be solved directly. So we discuss how to get the asymptotically optimal bandwidth. There are quite a few important techniques for selecting the bandwidth. such as cross-validation and plug-in bandwidth selectors. a conceptually simple technique, with theoretical justification and good empirical performance, is the plug-in technique.Another issue in multivariate local polynomial fitting is the choice of the order of the polynomial. Since the modeling bias is primarily controlled by the bandwidth, this issue is less crucial however. For a given bandwidthh, a large value of d would expectedly reduce the modeling bias, but would cause a large variance and a considerable computational cost. Since the bandwidth is used to control the modeling complexity, and due to the sparsity of local data in multidimensional space, a higher-order polynomial is rarely used. We use the local quadratic regression to indicate the flavor of the multivariate local polynomial fitting, that is to say, d=2.The third issue is the selection of the kernel function. In this paper, of course, we choose the optimal sphericalEpanechnikov kernel function, which minimizes the asymptotic MSE of the resulting multivariate local polynomial estimators, as our kernel function.


### 3.4. Estimating the Defocused Parameter

Twenty original images are chosen to train the model. The images are defocused artificially withR whose value ranging from 2 to 7. So the total number of training samples are 120. Then feature vectors are constructed using variables v1-5 of each image:(3.18)V̲=(v1,v2,v3,v4,v5).The defocused parameterR is the model output.When training samples{VTa̲,Ra}a=1120 are given, obtaining weights matrix B, according to the relationship between the V and R, then the defocused parameter R can be calculated using the trained model.


## 3.1. Multivariate Kernel Function

To localize data in them-dimension, we need a multikernel function. Generally speaking, a multivariate kernel function refers to an  m-variate function satisfying
(3.1)∫-∞+∞⋯∫-∞+∞K(x̲)dx̲=1.Here and hereafter, we use∫ to indicate multivariate integration over the m-dimensional Euclidean space.There are two common methods for constructing multivariate kernel functions. For a univariate kernelk(x), the product kernel is given by
(3.2)K(x̲)=∏i=1mk(xi),
and the spherically symmetric kernel is defined as
(3.3)K(x̲)=cK,mK(∥x̲∥),
where cK,m={∫K(∥x̲∥)dx̲}-1 is a normalization constant and ∥x̲∥=(x12+x22+⋯+xm2)-1/2. Popular choices of K(x̲) include the standard d-variate normal density
(3.4)K(x̲)=(2π)-m/2exp(-∥x̲∥22),
and the spherical Epanechnikov kernel(3.5)K(x̲)={d(d+2)Γ(m/2)4πm/2}(1-∥x̲∥2)+.
The latter is the optimal kernel, according to Fan and Gijbels [21] and Su [22].The localization in multivariate nonparametric regression is frequently carried out by the kernel weighting. LetH̲ be a symmetric positive-definite matrix called a bandwidth matrix. The localization scheme at a point x̲ assigns the weight
(3.6)KH̲(Xi̲-x̲),withKH̲(x̲)=|H̲|-1K(H̲-1x̲),
where |H̲| is the determinant of the matrix H̲. The bandwidth matrix is introduced to accommodate the dependent structure in the independent variables. For practical problems, the bandwidth matrix H̲ is taken to be a diagonal matrix. The different independent variables will be accommodated into different scales. For simplification, the bandwidth matrix is designed into H̲=hIm̲ (Im̲ denoting the identity matrix of order m).


## 3.2. Multivariate Predictor with Local Polynomial Fitting

Suppose that the input vector isV̲=(v1,v2,v3,v4,v5), the model is fitted by the function(3.7)R=f(V̲).Our purpose is to obtain the estimationR̂=f̂(V̲) of function f, this paper, we use the dth order multivariate local polynomial f(V̲) to predict the defocused parameter RT value based on the point VT̲ of the test image. The polynomial function can be described as(3.8)f(V̲)≈∑0≤|j̲|≤d1j̲!D(j̲)fi(VT̲)(V̲-VT̲)j̲=∑0≤|j̲|≤db̲j̲(VT̲)(V̲-VT̲)j̲,
where(3.9)m=5,j̲=(j1,j2,…,jm),j̲!=j1!j2!⋯jm!,|j̲|=∑l=1m‍jl,∑0≤|j̲|≤d=∑|j̲|=0d(∑j1=0|j̲|∑j2=0|j̲|⋯∑jm=0|j̲|)|j̲|=j1+j2+⋯+jm,V̲j̲=v1j1v2j2⋯vmjm,D(j̲)fi(VT̲)=∂|j̲|fi(V̲)∂v1j1∂v2j2⋯∂vmjm|V̲=VT̲,b̲j̲(VT̲)=1j̲!D(j̲)fi(VT̲).In the multivariate prediction method,VTa̲(a=1,2,…,A) denotes the trained image feature vectors. Using A pairs of (VTa̲,Ra), for which the values are already known, the coefficients of fi is determined by minimizing
(3.10)∑a=1A[Ra-∑0≤|j̲|≤db̲j̲(VT̲)(VTa̲-VT̲)j̲]2·KH̲(VTa̲-VT̲).For the weighted least squared problem, a matrix form can be described by(3.11)W̲1/2·Y̲=W̲1/2·X̲·B̲+ɛ̲,
where(3.12)Y̲=(y1,y2,…,yA)T,ya=Ra,B̲=(b̲0(VT̲),b̲1(VT̲),…,b̲d(VT̲))T,W̲=diag{KH̲(VT1̲-VT̲),KH̲(VT2̲-VT̲),…,KH̲(VTA̲-VT̲)},
and X̲ is the A×S(S=∑0≤|j̲|≤d|j̲|/j̲!),(3.13)X̲=(1(VT1̲-VT̲)1⋯(VT1̲-VT̲)d1(VT2̲-VT̲)1⋯(VT2̲-VT̲)d⋮⋮⋱⋮1(VTA̲-VT̲)1⋯(VTA̲-VT̲)d).We then have the least squared solution with multivariate local polynomial fitting as follows:(3.14)B̲̂=(W̲1/2X̲)†Y̲
or, when X̲TW̲X̲ is inverse, the estimation can be written by(3.15)B̲̂=(X̲TW̲X̲)-1X̲TW̲Y̲,
then we can get the estimation R̂T=f̂(VT̲)(3.16)R̂T=f̂(VT̲)=E1̲(X̲TW̲X̲)-1X̲TW̲Y̲,
where E1̲=(1,0,0,…,0)1×S.Computing theB̲̂ will suffer from large computational cost. we can use the recursive least squared method to reduce the computation complexity, and it is very powerful especially in the real time prediction problems. There are several important issues about the bandwidth, the order of multivariate local polynomial function and the kernel function which have to be discussed. The three problems will be presented in Section 3.3.


## 3.3. Parameters Selections

For the multivariate local polynomial predictor, there are three important problems which have significant influence to the prediction accuracy and computational complexity. First of all, there is the choice of the bandwidth matrix, which plays a rather crucial role. The bandwidth matrixH̲ is taken to be a diagonal matrix. For simplification, the bandwidth matrix is designed into H̲=hIm̲. So the most important thing is to find the bandwidth h. A too big bandwidth underparameterizes the regression function, causing a large modeling bias, while a too small bandwidth overparameterizes the unknown function and results in noisy estimates. In theory, there exists an optimal bandwidth hopt in the meaning of mean squared error, such that
(3.17)hopt=argminh∫(f(x̲)-f̂(x̲))2dx̲.But the optimal bandwidth cannot be solved directly. So we discuss how to get the asymptotically optimal bandwidth. There are quite a few important techniques for selecting the bandwidth. such as cross-validation and plug-in bandwidth selectors. a conceptually simple technique, with theoretical justification and good empirical performance, is the plug-in technique.Another issue in multivariate local polynomial fitting is the choice of the order of the polynomial. Since the modeling bias is primarily controlled by the bandwidth, this issue is less crucial however. For a given bandwidthh, a large value of d would expectedly reduce the modeling bias, but would cause a large variance and a considerable computational cost. Since the bandwidth is used to control the modeling complexity, and due to the sparsity of local data in multidimensional space, a higher-order polynomial is rarely used. We use the local quadratic regression to indicate the flavor of the multivariate local polynomial fitting, that is to say, d=2.The third issue is the selection of the kernel function. In this paper, of course, we choose the optimal sphericalEpanechnikov kernel function, which minimizes the asymptotic MSE of the resulting multivariate local polynomial estimators, as our kernel function.


## 3.4. Estimating the Defocused Parameter

Twenty original images are chosen to train the model. The images are defocused artificially withR whose value ranging from 2 to 7. So the total number of training samples are 120. Then feature vectors are constructed using variables v1-5 of each image:(3.18)V̲=(v1,v2,v3,v4,v5).The defocused parameterR is the model output.When training samples{VTa̲,Ra}a=1120 are given, obtaining weights matrix B, according to the relationship between the V and R, then the defocused parameter R can be calculated using the trained model.


## 4. Iterative Wiener Filter

Wiener filtering (minimizing mean square error) is commonly used to restore linearly-degraded images. To obtain optimal results, there must be accurate knowledge of the covariance of the ideal image. In this section, the so-called iterative Wiener filter [19, 23] is used to restore the original image.The imaging systemH is assumed to be linear shift invariant with additive, independent, white noise processes of known variance. the model for the observed image y is given in matrix notation by
(4.1)y=Hf+n,
where f is the ideal image. The optimal linear minimum mean-squared error, or Wiener restoration filter given by
(4.2)f̂=By,
where B=RffHT[HRffHT+Rnn]-1, requires accurate knowledge of Rff, the autocorrelation of ideal image f. However, in practical situations f is usually not available and only a single copy of the blurred image to be restored, y, is provided. In the absence of a more accurate knowledge of the ideal image f, the blurred image y is often used in its place simply because there is no other information about f readily available. The signal y is subsequently used to compute an estimate of Rff and this estimate is used in place of Rff in (4.2).The following summarizes the iterative Wiener filtering procedure.Step 1 (Initialization).
Usey to compute an initial (i=0) estimate of Rff by
(4.3)Rff(0)=Ryy=E{yyT}.Step 2 (Filter Construction).
UseRff(i), the ith estimate of Rff to construct the (i+1)th restoration filter B(i+1) given by
(4.4)Bi+1=RffHT[HRffHT+Rnn]-1.Step 3 (Restoration).
Restorey by the B(i+1) filter to obtain f̂(i+1), the (i+1)th estimate of f(4.5)f̂(i+1)=B(i+1)y.Step 4 (Update).
Usef̂(i+1) to compute an improved estimate of Rff, given by
(4.6)Rff(i+1)=E{f̂(i+1)f̂T(i+1)}.Step 5 (Iteration).
Incrementi and repeat Steps 2, 3, 4, and 5.


## 5. Experimental Results and Analysis

The experiments are carried out by using the Matlab image processing toolbox. The performance of the proposed image restoration algorithm has been evaluated using the classical gray-scaleMoon image, Coins image, Saturn image, and Tire image in Matlab toolbox. To verify the good ability of restoration of the proposed algorithm, one real blurred image is used for the deconvolution procedure. The results show our method is very successful for this kind of blurred image.In image restoration studies, the degradation modelled by blurring and additive noise is referred to in terms of the metric blurred signal-to-noise ratio (BSNR). This metric for a zero-mean M×N image is given by(5.1)BSNR=10log10{(1/MN)∑m=1M∑n=1Nz2(m,n)σv2},
where z(m,n) is the noise-free blurred image and σv2 is the additive noise variance.For the purpose of objectively testing the performance of linear image restoration algorithms, the improvement in signal-to-noise ratio(ISNR) is often used. ISNR is defined as(5.2)ISNR=10log10{∑m=1M∑n=1N[f(m,n)-y(m,n)]2∑m=1M∑n=1N[f(m,n)-f̂(m,n)]2},
where f(m,n) and y(m,n) are the original and degraded image pixel intensity values and f̂(m,n) is the restored true image pixel intensity value. ISNR cannot be used when the true image is unknown, but it can be used to compare different methods in simulations when the true image is known.In order to find the good performance of the proposed multivariate local polynomial Regression method (MLPR) compared with the RBF neural network algorithm (RBFNN) which is investigated in the literature [19], the same defocused blurred images are used for the experiments. Mean squared prediction errors are shown in Table 1. From Table 1, we can conclude that the prediction results of MLPR predictor are significantly better than the RBF neural network method in the same simulated data.Table 1
MSE using both methods.
Training  imageDifferent methodseMSEMoonRBFNN4.81×10-6MoonMLPR4.13×10-8CoinsRBFNN5.06×10-6CoinsMLPR3.97×10-8SaturnRBFNN6.62×10-6SaturnMLPR5.65×10-9TireRBFNN8.04×10-6TireMLPR7.19×10-8Figures3, 4, 5 and 6, in which the true images, blurred images and estimated true images are depicted in the left, middle and right column, respectively, illustrate how the method behaves in Moon, Coins, Saturn and Tire images. It is clear from Figures 3–6 that performance of the new method is effective in different images. Figure 7 also shows that the presented algorithm is good for real blurred image.Result ofMoon. True image (a), blurred image (b), estimated image(c), BSNR = 12.35, ISNR = 22.56.
(a)(b)(c)Result forCoins. True image (a), blurred image (b), restored image (c), BSNR = 11.22, ISNR = 23.14.
(a)(b)(c)Result forSaturn. True image (a), blurred image (b), restored image (c), BSNR = 13.17, ISNR = 24.31.
(a)(b)(c)Result forTire. True image (a), blurred image (b), restored image (c), BSNR = 11.56, ISNR = 22.09.
(a)(b)(c)Result for real blurred image. blurred image (a) restored image (b).
(a)(b)


## 6. Conclusions

A new method that is based on multivariate local polynomial regression model and iterative Wiener filtering for semi-blind restoration of blurred noisy images was proposed in this paper. Defocused parameter was estimated by a multivariate local polynomial regression model trained in wavelet domain. The main advantages of the proposed technique are that it is not only robust to noise because wavelet transform has an excellent denoising ability, but also effective to artificially and practically defocused blurred image. Restoration is successfully realized by the iterative Wiener filter, resulting in improved the image quality. The algorithm was justified via simulation and real image. Defocused image parameter can be successfully estimated by using trained model. Experimental results show the proposed algorithm is reliable and robust for defocused blurred image restoration.



---

*Source: 605241-2010-09-06.xml*