2-D Impulse Noise Suppression by Recursive Gaussian Maximum Likelihood Estimation

An effective approach termed Recursive Gaussian Maximum Likelihood Estimation (RGMLE) is developed in this paper to suppress 2-D impulse noise. And two algorithms termed RGMLE-C and RGMLE-CS are derived by using spatially-adaptive variances, which are respectively estimated based on certainty and joint certainty & similarity information. To give reliable implementation of RGMLE-C and RGMLE-CS algorithms, a novel recursion stopping strategy is proposed by evaluating the estimation error of uncorrupted pixels. Numerical experiments on different noise densities show that the proposed two algorithms can lead to significantly better results than some typical median type filters. Efficient implementation is also realized via GPU (Graphic Processing Unit)-based parallelization techniques.


Introduction
Median type filters are known as the effective methods to suppress impulse noise which often arises from sensor damage, malfunctioning or timing errors in signal acquisition [1]. Suppressing impulse noise routinely includes two stages -noise detection and noise removal [2][3][4]. The Standard Median Filter (SMF) is ineffective in cases with high noise densities, and the SMF processed images are often seriously degraded by detail loss and new introduced artifacts when the noise density is over 50% [5].
Several methods have been developed to overcome the drawback of SMF. Adaptive Median Filter (AMF) was proposed in [6] to improve the performance of median filter by only replacing corrupted pixel values by median values while leaving uncorrupted pixels unchanged. Additionally, the AMF recursively enlarges the filter window until one uncorrupted median is identified. The Decision Based Algorithm (DBA) proposed in [7] is similar to the AMF method except that DBA replaces the corrupted pixel value by a neighborhood pixel value when the median is tagged as a corrupted pixel. In [8], a method named Modified Decision Based Unsymmetric Trimmed Median Filter (MDBUTMF) was developed to improve DBA method by removing uncorrupted pixels in median calculation and setting the corrupted pixel value to mean of the selected window when all the window elements are corrupted. To achieve edge-preserving noise suppression, Chan and Nikolova put forward in [9] a twophase algorithm which includes a SMF based noise detection and an edge-preserving term regularized optimization. In [10], a multi-scale adaptive median-based filtering was developed to suppress salt-and-pepper noise.
In suppressing impulse noise, median-type methods suffer from the image degradation of detail blurring and staircase artifacts. This problem might be more serious for high noise density cases [11], or the case combined with structure blurring [12]. The reason can be ascribed to the fact that only the median information from one single neighboring pixel is selected for median filtering and no structure continuity information is utilized. To overcome this, we propose in this paper an approach termed Recursive Gaussian Maximum Likelihood Estimation (RGMLE) for impulse noise removal. We focus this study on noise removal stage, so only salt-and-pepper type impulse noise with pre-known intensities is considered here. The proposed RGMLE approach solves the maximum likelihood estimation (MLE) of the corrupted points through recursive variance-weighted averaging. A robust strategy of controlling the recursion in RGMLE is also developed by analyzing the estimation error of the uncorrupted pixels. Under the framework of RGMLE, we derived two noise removal algorithms RGMLE-C and RGMLE-CS by respectively using certainty and joint certainty & similarity to estimate the variance. GPU (Graphic Processing Unit)-based technique is also applied to accelerate the proposed processing. Extensive experiments are performed to show that the RGMLE-C and RGMLE-CS algorithms can produce significantly improved results than the comparative median type filters. Among all the methods, the RGMLE-CS algorithm demonstrates the best restoration performance.
The structure of this paper is organized as follows: In section II, we explain the RGMLE approach together with the two derived algorithms RGMLE-C and RGMLE-CS. Experimental results are given and discussed in section III. Section IV concludes this paper with a brief description of its contributions and some open questions for future work.

Reursive Gaussian Maximum Likelihood Estimation
Based on the theory of statistical estimation, the sample median corresponds to the maximum likelihood estimators of locations for independent and identically distributed (i.i.d.) observations obeying the Laplacian distributions with probability density function (pdf) shown in Fig.1 [13][14]. Let x 1 ,::::,x N be N i.i.d. Laplacian distributed observations with unknown location parameter m and the same variance 2g 2 , the MLE of m is denoted bym m, which maximizes the likelihood function: L x 1 ,::: This is equivalent to minimizing (2): Minimizing (2) with respect to m leads to the simple median m m = MEDIAN x 1 ,::::,x N ð Þ [13]. Here for an ordered list, the MEDIAN operator is defined as producing the middle value for odd number of values or the arithmetic mean of the middle two observations for even number of values. Assuming the Laplacian distribution of the observations in the 3|3 window, the 2D median filtering of the center pixel value m is in fact the maximizer of the likelihood function with respect to each location parameter m [15][16].
Likewise, with the Gaussian distribution shown in Fig. 1, we can build the likelihood function as follows [17][18]: L x 1 ,::: The value m minimizing (3) leads to the average operator: Then we can easily obtain the variance s 2 m m of m m: Here, equation (5) shows that the variance estimation for Gaussian distribution decreases linearly with the available observation numbers.
However, in 2-D impulse noise suppression, assuming an identical variance for all samples in (3), the averaging filtering often leads to a non-discriminative averaging operation and might tend to produce oversmoothing in the restored images(we would illustrate this in the following section on experiment). Improved performance can be expected if a discriminative estimation of variance is applied based on the different reliabilities of different observations. Furthermore, considering the fact that the estimation error (or variance) decreases as the available uncorrupted observations increase, the information from corrected intensities should be exploited in a recursive way to refine the final restoration. In this work, we develop a Recursive Gaussian Maximum Likelihood Estimation (RGMLE) for impulse noise suppression. In this 2-D noise suppression case, with y the noisy image and m the target restored image, to estimate each corrupted pixel m j we can rewrite the likelihood function with variances s 2 ij for each observation x ij in the filter window N j : L x 1j ,::: Here DN j D denotes the point number in filter window N j . In (6), spatially-adaptive variance is used based on joint certainty and similarity information of observations, which is to be elaborated in detail below. Minimizing (6),with respect to m j leads to a normalized variance-weighted average: Here, for each observation i, the weight w ij is a positive value calculated as the direct reciprocal of the variance s 2 ij : The recursive update of image m in (7) should be stopped at certain recursion number to avoid the deterioration caused by the successive involvement of those irrelevant uncorrupted pixels. But the original true image is not available to guide the recursion to stop at the point when the optimal restoration (the closest to the original true image) is reached. Here, we use the estimation error of the ground truth uncorrupted pixels to evaluate the restoration deterioration in recursion. After each recursion, we recalculate the intensities of the uncorrupted pixels using the restored image via (7), and analyze their peak signal-to-noise ratio with respect to the original ground truth uncorrupted pixels (PSNR_U, calculated via (9) in the below algorithm outline). The recursion is stopped once the PSNR_U stops to increase, which represents the time when the contribution from corrected pixels starts to be outweighted by the deterioration from the irrelevant uncorrupted pixels. This strategy works well for the cases with noise densities ranging from 10% to 70%. For the cases with very high noise densities (.70%) this strategy does not work well because the re-estimation of the uncorrupted pixels mainly come from irrelevant pixel intensities, and can not reflect the recursion deterioration. In this case a large number of recursion is required to get a stable restoration. An illustration will be given in the below experiment to validate this. We also assume intensity value 255 and 0 for the corrupted pixels in the original image y, and let x denote the image to be restored. Two binary matrixes M o and M u represent the originally uncorrupted pixel mask and the uncorrupted pixel mask in recursion. We should note that the corrupted pixels will receive some estimation and become uncorrupted pixels when they receive some estimation in the processing. So M u will be updated to an all-ones matrix when all corrupted pixels get estimated as the recursion proceeds. DMD, DM u D and DM o D denote the total pixel number in image, and the total non-zero pixel numbers in M u and M o , respectively. The flowchart of the RGMLE strategy is given in Fig.2. Here, a large recursion number T is set when the noise density is larger than 70%. We judge the PSNR increase for the estimation of the originally uncorrupted pixels by evaluating the error improvement C. The two algorithms RGMLE-C and RGMLE-CS are detailed as follows: Algorithm RGMLE-C In RGMLE-C algorithm, we consider the certainty measure as whether the pixels are corrupted in the original image y. In estimating each m j via (7), the variance s 2 ij is estimated based on the certainty measure of each observation x ij , and can be calculated as follows: where the pixel x ij is considered having no certainty if it is originally corrupted and has not yet been corrected (M u ij~0 ), and in this case the variance s 2 ij is set to z?. We set s 2 ij to 1 if the originally corrupted pixel x ij has obtained some certainty via the recursive estimation ((M u ij~1 )and(M o ij~0 )). As to those originally uncorrupted pixels (M o ij~1 ), they should be assigned a high certainty, then a low variance is assigned to them with certainty parameter l (l,1). The corresponding weight w ij can then be calculated as the reciprocal of variance s 2 ij : In re-estimating each uncorrupted pixel m _ j via (8), to accurately reflecting the restoration deterioration in recursion, we exclude the original uncorrupted pixels and only use the restored ones, and the following weight w _ ij is used: Algorithm RGMLE-CS Within a filter window for one specific point, two pixels belonging to the same structure tend to come from the same class, and should be assigned low variances (or large weights) in estimating each other. We borrow the patch similarity idea in [19], [20], and quantify the variance by the exponential function of the ' 1 norm distance between the patches centered at the two pixels.
Here we use ' 1 norm for its good edge-preserving property [21], [22]. By introducing this patch similarity information into the variance estimation in the above RGMLE-C algorithm, we calculate each variance s 2 ij : where the variance s 2 ij to estimate m j for pixel x ij is set to z? if the corrupted pixel has still not been corrected (M u ij~0 ). Also if the originally corrupted pixel has obtained some certainty through the correction in recursion ((M u ij~1 )and(M o ij~0 )), we set the variance to the exponential value of the ' 1 norm of the distance d ij between the uncorrupted pixels in the two 2D patches P j and P ij centered at point j and its neighboring point ij. In (14), two binary patch masks P u ij and P u j denote the positions of the uncorrupted points and are used to exclude the uncorrupted elements in P j and P ij in distance calculation. Practically, parameter h in (13) should be set to modulate the attenuation of the exponent function with respect to distance d ij . As to the originally uncorrupted pixel (M o ij~1 ), we also lower the value s 2 ij via multiplying it by a certainty parameter l (l,1). Then, the weight w ij is given as Also, in re-estimating each uncorrupted pixel m _ j , we exclude the original uncorrupted pixels and only use the restored ones. In this case, we use the weight w where, two binary patch masks P P o ij ( P P o ij~1 ÀP o ij ) and P P o j ( P P o j~1 ÀP o j ) are used to exclude the original uncorrupted pixels in distance calculation.

A. Experiment configuration
Three 512|512grayscale images (Lena, Pepper and MRI in the below Fig.3-Fig.8)  format and was collected from a head scanning of a 57 years old patient in the first hospital of Nanjing. The intensity range of the MRI image is [0, 1332]. The patient in this manuscript has given the written informed consent of using the MRI image data. A wide range of salt-and-pepper noise with density from 10% to 90% (10% increments) is considered. The same simulation as the analysis in above Section II is used. The AMF, DBA and MDBUTMF methods are implemented based on [6][7][8] for comparison. The MDBUTMF method is actually the standard median filter when any pixel in window is uncorrupted, and switches to averaging filtering otherwise [8]. Replacing the median operator in MDBUTMF by the average operator leads to an average based method and we termed it MDBUTAF. We implement this MDBUTAF method to provide a fair comparison with median and average filtering. We also compare our method with the method based on steering kernel regression (SKR method) in [23], which has demonstrated good performance in denoising, interpolation and super-resolution. The parameters involved in different methods are all set to give stable results with the highest peak signal-to-noise ratio (PSNR, calculated via (18)). We use a 3|3 filtering window N in the AMF, DBA, MDBUTMF and MDBUTAF methods. For the SKR method, the global smoothing parameter and the kernel size were manually tuned to give the results with the highest PSNR. For the RGMLE-C and RGMLE-CS algorithms, several parameters need to be set, namely the size of filtering window W, the certainty parameter l and the recursion limit T. Here 3|3 and 5|5 respectively used in the RGMLE-C and RGMLE-CS algorithms. For RGMLE-CS algorithm, the 5|5 of more distant pixels via the above similarity-based weighting strategy. Also two additional parameters (the smoothing parameter h and the size of patch P) need to be set for the RGMLE-CS method. In RGMLE-CS algorithm, a |7 patch P is set to give an effective characterization of local structures. And the smoothing parameter h is set to 3.5 for Lena and Peppers images, and to 10 for the MRI image with a larger intensity range. Specially, for RGMLE-CS algorithm, we set P and N to 1|1 and 3|3 in the first recursion to overcome the lack of structure information in the original corrupted noisy image. We set the certainty parameter l in (11) and (15) to 0.2 to increase the weights of the uncorrupted pixels. The recursion number limit T is set to 50 to include a long recursion in restoration using RGMLE-C and RGMLE-CS algorithms.
The pixel-wise operation (7) in RGMLE-C and RGMLE-CS algorithms allows for easy parallelization using GPU technique [24], [25]. Under Compute Unified Device Architecture (CUDA) framework, we set the total number of blocks in grid to the row size of the image, and the total number of threads in each block to the column size of the image. In implementing the RGMLE-C and RGMLE-CS algorithms, all threads in the block-grid structure execute simultaneously to perform all the pixel-wise operations, and only half of the distance d is calculated because of the symmetry property of distance calculation (d ij = d ji ). All the images were processed in a PC workstation (Intel Core 2 Quad CPU and 4096 Mb RAM, GPU (NVIDIA GTX465)) with Visual C++ as the developing language (Visual Studio 2008 software; Microsoft).
Restoration performances are quantitatively measured by PSNR and the structural similarity index comparisons (SSIM) proposed in [26]: where, m I denote the restored image and the original true image. m m and m I denote the mean intensities of images mand I, s m and s I are the standard deviation of images m I, s mI is the covariance of images m and I, c 1~( K 1 L) 2 and c 2~( K 2 L) 2 with L the dynamic range of the pixel values (255 for 8-bit grayscale images), K 1 and K 2 are set to 0.01 and 0.03 as suggested in [25]. Fig. 3 and Fig. 4 illustrate the restoration results of Lena image with 80% and 90% noise densities for different methods. Likewise, Fig. 5 and Fig. 6 present the results for Pepper image, and Fig. 7 and Fig. 8 present the results for MRI image. The calculated PSNR and SSIM with respect to the original true images are also given in the captions below the figures. The original true images and the corrupted images are illustrated in (a) and (b) for reference. In Fig.3-Fig.8, the denoised images from different methods are given in (c)-(i), and (az)-(iz) denote the zoomed local parts for the corresponding images in (a)-(i).We can see that the traditional median-type filters (AMF, DBA and MDBUTMF methods with results (c)-(e) in Fig.3-Fig.8) can not effectively estimate corrupted points, and tend to introduce many artifacts in restoration. Comparing the images (e) and (f) in Fig.3-Fig.8, we find that under high noise densities the MDBUTAF method works much better than MDBUTMF method (over 2.5dB PSNR enhancement). But we also find the MDBUTMF method with non-discriminative averaging operation can not effectively restore image details and overcome oversmoothing. Comparing the images (f) and (g) in Fig.3-Fig.8, we can see that, by recursively incorporating the certainty measure into estimation, the RGMLE-C algorithm produces notable improvement with more than 2.5dB PSNR enhancement over the MDBUTAF method. Obviously, the proposed RGMLE-CS algorithm gives good performance in terms of both noise suppression and detail preservation(see images (i) in Fig.3-Fig.8). With patch similarity information used, some tiny features (e. g. the hair structures in Lena image and the brain structures in MRI image) can be successfully restored by the RGMLE-CS algorithm. Among all the methods, the RGMLE-CS algorithm leads to the results with significantly higher PSNR values than other methods except the SKR method. We note that SKR method can give higher PSNR values than the RGMLE-CS algorithm in Lena image (80% noise density), Peppers image (90% noise density) and MRI image (80% and 90% noise density). But in terms of SSIM, which is deemed a more robust quality metric than PSNR, the proposed RGMLE-CS algorithm provides the best results. Compared to the results from SKR method, the processed results from RGMLE-CS algorithm present structures with higher contrast in the case of 80% noise density (see the arrows in images (iz) in Fig.3, Fig.5 and Fig.7), and perceptively comparable results under 90% noise density.

B. Restoration Result
In Fig. 9, we summarize the PSNR and SSIM values of the restored images for noise densities ranging from 10% to 90%. We can observe in the plots that the obtained PSNR and SSIM values decrease as noise level rises for all the methods, and our proposed RGMLE-C and RGMLE-CS algorithms obtain significantly higher PSNR and SSIM than the other methods except the SKR method in almost all noise densities. This advantage becomes more prominent as noise density increases. The SKR method can sometimes provide higher PSNR values than RGMLE-CS algorithm under high noise densities (80% and 90%), but has poor performance in the cases of low noise densities. We can see the proposed RGMLE-CS algorithm can give the best results in terms of SSIM for all the noise densities.  (18)) with respect to the recursion numbers (T = 50) for the whole restored image (denoted by PSNR_I, blue lines) and the PSNR_U (via above (9)) of the reestimated uncorrupted pixels (green lines) under different noise densities. Lena image with the same parameter setting as above is used in this validation. Fig. 10 shows that the PSNR-I decreases when the recursion proceeds across one specific point, which shows the deterioration will appear in recursion for the proposed RGMLE-C and RGMLE-CS methods. In Fig. 10, we can note a good match between the highest obtained PSNR of the restored images and that of the re-estimated uncorrupted pixels, but such match degree tends to obviate as the noise density increases to over 70%. The reason of the match deviation for noise density higher than 70% is that under very high noise densities the re-estimation of the uncorrupted pixels mainly relies on irrelevant pixel   intensities, which can not accurately reflect the overall deterioration in restoration. This observation validates the above proposed stopping strategy for the RGMLE algorithm. Table 1 lists the CPU time costs in implementing all the methods in the above experiments with different noise densities. The CPU time costs for implementing both the parallelized RGMLE methods (parallelized RGMLE-C and RGMLE-CS) and the original serial versions (serial RGMLE-C and RGMLE-CS) are included. The AMF and SKR methods require much more intensive computation than other median-type filters because the AMF method recursively enlarges the filter window to obtain an uncorrupted median and the SKR method requires a computationally costly operation of matrix inversion. Table 1 also shows that the proposed RGMLE-C and RGMLE-CS methods require more computation than other methods, and about 150 times acceleration is achieved after GPU based parallelization. We can also notice that the required computation cost of RGMLE-C and RGMLE-CS methods increases as the noise density increases because more recursions need to be run based on the recursion stopping rule.

Conclusions
This paper proposed a Recursive Gaussian Maximum Likelihood Estimation (RGMLE) for 2-D impulse noise removal. Assuming Gaussian distribution for neighboring intensities, the proposed RGMLE approach solves the maximum likelihood estimation (MLE) of corrupted points through recursive varianceweighted averaging. Simulation in Section II reveals that Gaussian distribution assumption is theoretically preferable to Laplacian distribution assumption in estimating corrupted points under high noise densities. In validation we use salt-and-pepper type impulse noise with pre-known intensities. Two noise removal algorithms RGMLE-C and RGMLE-CS are respectively derived by using certainty based or certainty&similarity based variance estimations.
An effective stopping strategy for recursion is also developed based on the estimation error of uncorrupted pixels using restored image information. Extensive experiments have been performed to show the RGMLE-C and RGMLE-CS algorithm can obtain significantly better results than the comparative median type filters. Among all the methods, the RGMLE-CS algorithm achieves good restoration performance through the incorporation of patch similarity information. The application of the proposed RGMLE methods can be easily extended to suppress general impulse noise with random values when combined with the detection methods proposed in [3], [4], [11], [27]. An extension to Gaussian noise suppression would also be performed by considering a spatially varying noise corruption degree for each pixel.
With GPU parallelization, the RGMLE-CS algorithm can be efficiently applied and should be suggested for impulse noise suppression. A constant smoothing parameter h is currently used in the RGMLE-CS algorithm, and we can improve the algorithm by optimizing the estimation of parameter h [28], [29]. The denoising performance of the RGMLE-CS algorithm can also be further enhanced via the combination with multiscale or compressed sensing techniques [30], [31]. Current recursion stopping method is based on the knowledge of estimation error of uncorrupted pixels, and from above we know that this strategy does not apply to the case of very high noise density. More effective method is to be studied to control recursion to obtain the exactly optimal restoration for different noise levels. Computation efficiency still holds a future issue to be considered. Though significant acceleration has been obtained by GPU parallelization technique, further computation optimization is still required to realize realtime processing of the RGMLE-CS algorithm.