Image Denoising Using Fourth Order Wiener Filter with Wavelet Quadtree Decomposition

Visual information transmitted in the form of digital images is becoming a major method of communication in the modern age, but the image obtained after transmission is often corrupted with noise. The received image needs processing before it can be used in applications. Image denoising involves the manipulation of the image data to produce a visually high quality image. This paper uses the fourth order nonlinear wiener filter with wavelet quadtree decomposition and median absolute deviation. It will be shown that this new algorithm is comparable to other algorithms like BM3D, LPG-PCA, and KSVD.


Introduction
In digital image processing, image denoising is a very important issue.Certain degradation will happen to the transmitted images [1].This degradation can be noise or blurring.Blurring is a kind of bandwidth reduction due to some errors related to methods of capturing the photos.However, noise is related to the transmission medium errors or errors of measurements [2].
Image denoising is the estimation of the original image from the noisy image.Many methods have been used.Signal to noise ratio and mean square error are used as performance measures.The denoising concept can be represented by a system where the input is the noisy image and the output is the reconstructed image.We assume that the noise is an additive Gaussian noise given by where  is the pixel of the image to which is added the noise,  is the mean of the Gaussian noise, and  2 is the variance of the noise.Usually the variance of the noise is unknown.It is estimated via a method called median absolute deviation (MAD) given by MAD  =  med         − med         . ( The MAD has the best possible breakdown point (50%, twice as much as the interquartile range), and its influence function is bounded, with the sharpest possible bound among all scale estimators [3].The constant  is needed to make the estimator consistent for the parameter of interest.In the case of the usual parameter  at Gaussian distributions, we need to set  = 1.4826 [3].
1.1.Denoising Using Wavelet Transform.Wavelet coefficients calculated by a wavelet transform represent change in the time series at a particular resolution.By considering the time series at various resolutions, it is then possible to filter out noise.Wavelet thresholding is explained as decomposition of the data or the image into wavelet coefficients, comparing the detail coefficients with a given threshold value, and shrinking these coefficients close to zero to take away the effect of noise in the data [4].The image is reconstructed from the modified coefficients.This process is also known as the inverse discrete wavelet transform.During thresholding, a wavelet coefficient is compared with a given threshold and is set to zero if its magnitude is less than the threshold; otherwise, it is retained or modified depending on the threshold rule.Thresholding distinguishes between the coefficients due to noise and the ones consisting of important signal information.The choice of a threshold is an important point of interest.It plays a major role in the removal of noise in images because denoising most frequently produces smoothed images, reducing the sharpness of the image.Care should be taken so as to preserve the edges of the denoised image.There exist various methods for wavelet thresholding, which rely on the choice of a threshold value.Some typically used methods for image noise removal include Visushrink and Sureshrink [4].
In Visushrink or Sureshrink, the image is first subjected to a discrete wavelet transform, which decomposes the image into various subbands [4] as shown in Figure 1.
The subbands HHk, HLk, LHk,  = 1, 2, . . .,  are called the details, where k is the scale and j denotes the largest or coarsest scale in decomposition.Note that LLk is the low resolution component.Thresholding is now applied to the detail components of these subbands to remove the unwanted coefficients, which contribute to noise.And as a final step in the denoising algorithm, the inverse discrete wavelet transform is applied to build back the modified image from its coefficients.
It should be noted that Visushrink uses a hard thresholding rule with a universal threshold ℎ = √2 log , where  2 is the variance of the noise and  is the signal size [4].
The Sureshrink is a combination of the universal threshold and the SURE threshold.This method specifies a threshold value  for each resolution level j in the wavelet transform which is referred to as level dependent thresholding [5].
In [6] a new threshold operator has been used.This method consists of thresholding the wavelet coefficients using fourth order polynomial and gave better results than the previous thresholding techniques.It should be noted that our fourth order Wiener filter is performed directly on the wavelet coefficients without thresholding.

Denoising Using Neighboring Wavelet
Coefficients.Wavelet denoising by thresholding tends to kill too many wavelet coefficients that might contain useful image information.As a solution for this problem, wavelet thresholding is done by incorporating neighboring co-efficients [7].The idea of considering the influence of other wavelet coefficients on the current wavelet coefficient to be thresholded is of great importance.A large wavelet coefficient will probably have large wavelet coefficients at its neighbor locations.The reason is that wavelet transforms produce correlated wavelet coefficients [7].
In the following wavelet denoising scheme, the neighboring coefficients are incorporated into the thresholding process.Suppose that  , is the set of wavelet coefficients of the noise 1D signal: If ( 3) is less than or equal to  2 , then the wavelet coefficient  , is set to zero.Otherwise,  , =  , (1 −  2 / 2 , ), where  = √2 2 log , and  is the length of the signal [7].

KSVD or K-Means Clustering
Process.Nowadays, it is important to search for sparse representations of signals using a dictionary matrix  ∈  × that contains  prototype signal atom for columns, {  }  −1 , and a signal  ∈   can be represented as a sparse linear combination of these atoms.The representation of  may either be exact  = , or approximated.The vector  ∈   contains the representation coefficients of the signal .Extraction of the sparse representation is a hard problem that has been extensively investigated in the past few years.It was assumed that the dictionary is known and fixed.Here the issue is to design the proper dictionary in order to better fit the sparsity model imposed.A full description of this algorithm is described in Figure 3 [8].
Each iteration (containing the sparse coding and the dictionary update) improves the denoising results because in this algorithm the work is to optimize || − || 2  2 , the higher the number of iterations, the best performance is achieved [9].

Image Denoising by Sparse 3D Transform-Domain Collaborative
Filtering.This is done by grouping similar 2D fragments of the image into 3D data arrays called "groups." In order to deal with 3D groups, collaborative filtering is used, and it includes 3 successive steps: (1) 3D transformation of a group, (2) shrinkage of transform spectrum, (3) inverse 3D transformation.
After doing these steps, an array of jointly filtered 2D fragments is obtained.Due to the similarity between the grouped blocks, the transform can achieve a highly sparse representation of true signal, so that the noise can be well separated by shrinkage [10].In this way, the collaborative filtering reveals even the finest details shared by grouped fragments and at the same time it preserves the essential unique features of each of individual fragments.In this case, denoising using BM3D filter is not effective.That is why in order to increase the sparsity of the true signal and improve the BM3D filter, similar adaptive-shape neighborhoods are grouping together and PCA is a part of the 3D transform used for collaborative filtering [11].This method is called BM3D-SAPCA.

Two-Stage Image Denoising by Principal Component Analysis with
Local Pixel Grouping.BM3D algorithm achieves important results in denoising, but the problem is that its implementation is complex [12].Another algorithm has appeared with high efficiency and less complexity, the LPG-PCA algorithm or PCA based denoising method with local pixel grouping.In the PCA domain, noise can be removed from an image due to the reservation of only the most significant principal components [12].
The first stage is used to remove the most of noise in the image and the second stage to improve the output.The 2

Wiener Filter in the Wavelet Domain
In this method, wavelet coefficients are considered conditionally independent Gaussian random variables.The noise is also considered as an independent stationary zero mean Gaussian variable [14].Let us assume the relation between the noisy image and the original one as follows: where  , represent the coefficients of the noisy image in the wavelet domain,  , represents the coefficient of the original image, and  , represents the coefficient of the noise [14].
, 's are normal with zero mean and variance  2  .Due to the decorrelation property of the orthogonal wavelet transform, the signal components  , are uncorrelated.
An optimal linear estimator for a signal in additive noise is formed as  , =  , + ; the noise and the signal are assumed to be independent and  and  are chosen in a way to minimize the mean squared estimation error [15]: In order to get the minimum of, we have to get the derivatives with respect to  and  and set them equal to zero: This will give the following.
As a result we have the following.
Consider  , =   , , where  is the best linear estimator of the signal component ; because the noise and the signal are independent, we can say that  (  (,  2 ) .
Without loss of generality, we can assume that the { 2 , }'s can be determined by averaging the squared values of  , in a window centered at (, ).This information can be expressed as = (2 + 1) 2 is the number of coefficients in the kernel: Restricting the values of ( 2 , ) to only positive values, the numerator of the above equation takes the form: ŝ, = max (0, )  , . (10)

Fourth Order Nonlinear Wiener
Filtering with QTD 3.1.Quadtree Decomposition.Applying QTD to an image will splits a parent block into four children blocks if the intensity gradient within the block is greater than a predefined threshold.The decomposition will stop when the final blocks are composed of pixels or coefficients that are close to each other, or the difference between the maximum value and the minimum value is less than a certain small threshold [3].An example is shown in Figure 2.
In Figure 2, the right image represents the QTD of Lena image.It is clear how the values in each blocks are close and how the block sizes differ from smooth region to region where there are edges.

Summary of the Wavelet QTD + Wiener Algorithm.
QTD is used with the wiener filtering algorithm in order to get a better performance.Instead of filtering the same block size, different block size related to the quadtree decomposition is used.The algorithm can be summarized as follows: (1) Apply the discrete wavelet transform to the noisy image.It should be noted that we have used: The dwt2(I, ' d b4') matlab functions which uses the Daubechies Db4 wavelet.

It decomposes the image I into four subbands LL, LH, HL, and HH.
(2) Apply the QTD to each of the high frequency subbands (LH, HL, and HH).It should be noted that we have used: The qtdecomp(b,Threhold) matlab function which returns a sparse matrix S containing the top left corner of each block coordinate and the size of each block.This function splits a block b if the maximum value of the block elements minus the minimum value of the block elements is greater than threshold.Threshold is specified as values between 0 and 1.We have used a threshold of 0.2.The default value of the minimum block size is used.It should be noted that the effect of the variance of added noises will be decreased using the fourth order Wiener filter.
(3) Apply the Wiener filter on each variable size block: We start with the center of the first block as the pixel to be estimated, after estimating the value, we shift the block by one.
It should be noted that we have used the fourth order filter because of its good PSNR performance with different noise variances.The first order filter for example, will give good PSNR for small values of noise variances.Its performance degrades when the noise variances increase.
(4) Apply the inverse wavelet transform on the filtered image.

Fourth Order Nonlinear Wiener Filter.
In order to get the coefficients of this filter, the steps are summarized as follows: The filter output is given by Ŝ, =  4 , +  3 , +  2 , +  , + ; the error function is given by  (, , , , ) =  {[ 4  , +  3 , +  2 , In order to get MMSE, we have to get the partial derivatives of ( 11 The solution of the system above is as follows: However, we only have the noisy image, not the noise.Knowing that the noise is Gaussian noise with zero mean leads to where (p-1)!! is the multiplication of all the odd numbers between (p-1) and 0. This will give  ( 2 ) = 2

Conclusion
In this paper, we presented a new denoising method: the fourth order nonlinear wiener filter with wavelet quadtree decomposition and median absolute deviation.It is based on It was shown that this algorithm is comparable to other algorithms like BM3D, LPG-PCA, and KSVD algorithms.

1. 5 .
BM3D-SAPCA.Square image blocks containing fine image details, or sharp, or edges are examples where a nonadaptive transform is not able to deliver a sparse representation[11].

Figure 3 :
Figure 3: The test images used.

Figures 4 and 5
Figures 4 and 5 show the denoising results of the different methods for the "house" and "Lena" images using a variance of 400.(From top to bottom, left to right: Original image, Noisy image, KSVD, BM3D, LPG-PCA, 4th order Weiner results).
(a) applying the discrete wavelet transform to the noisy image, (b) applying the QTD to each of the high frequency subbands, (c) applying the 4th order Wiener filter on each variable size block.,(d) applying the inverse wavelet transform.

Table 1 :
Comparison between the 4 h order Wiener and other methods.