Abstract
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.
1. 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 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 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 [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 HH, HL, LH, are called the details, where is the scale and denotes the largest or coarsest scale in decomposition. Note that LL 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 , where 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.
1.2. 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 coefficients [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 , then the wavelet coefficient is set to zero. Otherwise, , where , and is the length of the signal [7].
1.3. 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, , 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 , the higher the number of iterations, the best performance is achieved [9].
1.4. 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.
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].
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.
1.6. 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 stages have the same principle, only the noise parameter will be changing [13].
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 . 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. Consider . Assume that and this implies that . Consider . And: Consider . 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
Without loss of generality, we can assume that the ’s can be determined by averaging the squared values of in a window centered at . This information can be expressed as
is the number of coefficients in the kernel:
Restricting the values of to only positive values, the numerator of the above equation takes the form:
3. 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.
3.2. 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:
(3)Apply the Wiener filter on each variable size block: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.
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.
(4)Apply the inverse wavelet transform on the filtered image.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.
3.3. 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 ; the error function is given by In order to get MMSE, we have to get the partial derivatives of (11) as follows: Now the system becomes The solution of the system above is as follows: where However, we only have the noisy image, not the noise. Knowing that the noise is Gaussian noise with zero mean leads to where is the multiplication of all the odd numbers between and 0.
This will give
4. Implementation and Results
We have done a comparison between the fourth order nonlinear wiener filtering with QTD and the BM3D, LPG-PCA, and KSVD algorithms. Four different images (house, Cameraman, Pepper, and Lena) were used with 4 different variances. The PSNR’s (18) between the original image and the reconstructed image of these methods are shown in Table 1: 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).
5. 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(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.
It was shown that this algorithm is comparable to other algorithms like BM3D, LPG-PCA, and KSVD algorithms.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.