A Laplacian based image filtering using switching noise detector

This paper presents a Laplacian-based image filtering method. Using a local noise estimator function in an energy functional minimizing scheme we show that Laplacian that has been known as an edge detection function can be used for noise removal applications. The algorithm can be implemented on a 3x3 window and easily tuned by number of iterations. Image denoising is simplified to the reduction of the pixels value with their related Laplacian value weighted by local noise estimator. The only parameter which controls smoothness is the number of iterations. Noise reduction quality of the introduced method is evaluated and compared with some classic algorithms like Wiener and Total Variation based filters for Gaussian noise. And also the method compared with the state-of-the-art method BM3D for some images. The algorithm appears to be easy, fast and comparable with many classic denoising algorithms for Gaussian noise.


Introduction
Denoising is one of the most important issues in image processing. The most popular noise removal methods are Adaptive Median Filtering (AMF), Total Variation (TV) based algorithms, Kernel based methods, Bilateral and Guided filtering and recently BM3D state-of-the-art in natural image denoising. In this work, we introduce a noise removal approach using a local noise estimator. We use the noise estimator in a minimization energy functional scheme. We obtain an iterative image denoising process using Laplacian. Denoising can be seen as adding values of pixels with their relative Laplacian weighted by local noise estimator.
Section 2 is related work. Section 3 represents our idea for denoising. We show how we can define a noise estimator using the sign of change in intensity of pixels in a 3x3 window. Using local noise estimator modified by Gaussian weight, we define an energy functional, drive the final equation and use it in an iterative denoising process. We show that although Laplacian is known as edge detector, it can be used for noise removal purposes. The algorithm is implemented in section 4. In section 5, results are shown and described. Moreover, figures of the denoising method are shown in comparison with Wiener, ROF and BM3D filters based on TV value and visual performance.

Related work
A Total Variation based noise removal method (ROF) (Rudin et al. 1992) defines an energy functional that preserves edges of the image and smoothes Gaussian noisy area, based on the total variation norm minimizer. The TV regularization technique is a suitable method that can be extended to different noisy conditions such as Laplace and Poisson (Chan and Esedo Lu 2005;Li et al. 1994). The Split Bregman method (Goldstein and Osher 2008) is fast, reliable and extendable to different models of noise distribution. Split Bregman is a basic and effective tool in solving many functional-based problems such as Compressed Sensing (CS) (Candès et al. 2006). In recent years, the usage of kernel-based techniques in image denoising has developed the quality of noise removal results. The image used in kernel functioning is called the guidance image. One of the most popular approaches using the guidance image is bilateral filter (Petschnigg et al. 2004). Other important kernel-based methods are Dataadaptive kernel regression (Takeda et al. 2007), Non-Local Means (Buades et al. 2005) and Optimal Spatial Adaptation (Kervrann et al. 2006). Another novel state of the art method was recently introduced as guided filter (He et al. 2010). The US patent 6229578 "Edge Detection Based Noise Removal Algorithm" (Acharya et al. 2001) is a denoising method based on using edge detector. This method removes noise by distinguishing between edge and non-edge pixels and applying a first noise removal technique to pixels classified as non-edge and a second noise removal technique to pixels classified as edge pixels. BM3D is a well-engineered algorithm which represents the current state-of-the-art method for denoising images corrupted by Additive White Gaussian Noise (AWGN) (Dabov et al. 2007;Dabov et al. 2006;Dabov et al. 2009;Chen and Wu 2010). In another strategy, denoised image is considered as a linear combination of the original image and its average when the coefficients are determined by an edge detector (Ranjbaran et al. 2013).

Methodology
Our methodology is based on using a local noise estimator in an energy functional minimizing scheme. Here we start with explaining our idea to define a local noise estimator. Consider Figure 1 where u(x, y) is pixel intensity. A noise in x direction is estimated when the sign of the change of the image intensity for two adjacent pixels is in opposite direction. Taking two x direction gradient components g x and g x + Δx we write: A noise is detected when g x < 0, g x + Δx > 0 or g x > 0, g x + Δx < 0. Generally we have: −g x g xþΔx > 0 Using Heaviside function H the noise detector in x direction can be defined as: Noise appears in two conditions in an image, first as ideal noise, shown in Figure 1 and second as correlated noise ( Figure 2). To distinguish these two cases we need to assign a weight to the detected location. Because SWN x acts as a switch its value is between zero and one. For ideal noisy pixels in the image plane, independent of the noise intensity, the value of weight should be one.
For the cases where noise is added to an edge as shown in Figure 2 the weight should be lower than one. To find a measure of the weight we consider that for ideal case |u(x + Δx, y) − u(x − Δx, y)| = 0 independent of noise intensity, but for non-ideal noise this difference is not zero. Then we can use a Gaussian weighting for the detected noise as the following equation: Where: As noise is added to the image in two directions we should drive the similar equations for y components. By using similar notations we find the final switching noise estimator as: To find a denoising way using SWN we define a measure of noise intensity in the image plane. Since SWN ≥ 0 the intensity of the noise in the noisy image can be defined as: To reduce noise, we define the following energy functional and try to minimize it: where u and u 0 are denoised and noisy images respectively and the first part is regularization term. The energy Figure 1 Two ideal noisy pixels.
functional is minimized in Appendix A. Although there are some techniques for computing u by Equ. 47, this way is a bit sophisticated and difficult to implement. To find a more simple equation we add λ f(u 0 ) to the two sides of the equation and write: Based on SWN operation and iteratively computing f(u 0 ) and f(u) in the locations classified as noise, we approximate f(u 0 ) = f(u). Then the final equation we used in our implementation is: This relation can be interpreted as follows: Because u is disturbed by f(u) and creates noisy image u 0 (Equ. 47), a similar process can restore u from u 0 (Equ. 9). By decreasing noise after some cycles of iteration f(u 0 ) goes to a small value. Equ. 47 presents a noise cancellation method based on using Laplacian value. The algorithm decreases the noise by adding the pixels value with Laplacian that weighted by SWN. Laplacian has been known as a common second-order edge detector but it has considerable value in noisy condition. Block diagram of the method is demonstrated in Figure 3.
Adding the intensity of pixels with the relative Laplacian is an averaging process. In an iterative process we can generally consider the evolution equation as: where Δt is the evolution timing step. For the cases that SWN = 1 the evolution equation is:  Then we have the following first order differential equation: When is the average value of u. The timing response is: where λ is the time constant. For λ ≥ 0 the denoised image u exponentially approaches to it steady state value ū in the locations where SWN is high. To find the constraints on λ we note that because 0 ≤ u ≤ 1 and − 1 ≤ f (u 0 ) ≤ 1, we have two constrains on λ as: As 0 ≤ λ ≤ 1 − u 0 and we choose ≥ 0, maximum value for λ is: It is predictable that in implementing the algorithm by large number of iterations λ must be a positive small value. In real condition SWN is not constant and reduced during the evolution.

Implementation
We have implemented our method in Matlab for Gaussian noise and evaluated it on different images. For implementation the Heaviside function is approximated by inverse tangent function: According to Equ.12 λ controls timing response and the value is between 0 and 1. So choosing the middle value can be suitable. Based on experimental results C = 3 is found as an appropriate value for noisy cases. We use a 3x3 window for simplicity and fast computing. The method is implemented by 10 numbers of iteration. The algorithm of implementation can be shown as the following steps: 1. Setting parameters : window size 3×3, λ = 0.5 2. Reading image u 0 3. Adding zero mean Gaussian Noise (imnoise code) 4. Computing SWN for current pixel

Computing Laplacian for current pixel
6. Updating u = u 0 + λ SWN ∇ 2 u 0 for the current pixel 7. Going to step 4 and continuing 8. Finishing when the whole of the image is scanned for ten times.

Results and discussions
We implemented our method in two noisy conditions (0.005 and 0.1 noise variance using imnoise matlab code for Gaussian noise) and compared it with ROF model using evolve2D code with 10 iterations and Wiener filter using wiener2 ('image', [3 3]) matlab code. f(u 0 ) can be interpreted as a noise intensity pattern in the image plane. Similar to ROF model in which div ∇u 0 ∇u 0 j j can be used as a measure of noise variance (Dabov et al. 2007), ‖f(u 0 )‖ is related to noise intensity. An example is shown in Figure 4. The results including noisy and denoised images are demonstrated in Figure 5 for Boat, Figure 6 for Man and Figure 7 for House. Figures 8, 9 and 10 show the results for Cameraman, Lena and Barbara respectively. TV for noisy and denoised images is shown in Tables 1 and 2. TV of the denoised image is totally comparable with ROF model. The computation time of our model is equal to ROF method with 10 iterations. f(u 0 ) has a decreasing behavior during the filtring time.
At start, f(u 0 ) is high in noisy pixels identified by SWN. Because image intensity is reduced by Laplacian value, f(u 0 ) tends to zero after some iteration. An example of such response is shown if Figure 11. Number of iterations is the only parameter that controls smoothness and running time. Large number of iteration makes the image blurry. Since f(u 0 ) goes to zero after some iteration, we can implement the algorithm adaptively until f(u 0 ) reaches a small value. This value should be estimated experimentally to make the algorithm applicable for every noisy condition. We also compared our method with BM3D technique for Boat, House and Cameraman in two SNR 25 and 75. The results are shown in Figures 12, 13 , 14, 15, 16, 17. The denoising technique used in this work can be generalized to other applications for reducing any feature of the image (F(u)) in other image processing tasks such as deblurring. A proposing framework can be written as the following energy functional: Such that: and λ is a control parameter determined experimentally. Minimizing the energy functional results the final relation for implementation.
As an example, for deblurring, Detector(F) and Wiegth (F) can be suggested as: where the blurring locations identified when the sign of change in intensity of two adjacent pixels is equal (g x g x + Δx > 0 , g y g y + Δy > 0), and weighted by difference in their gradient components for x and y directions. Computing ∂J ∂u ¼ 0 the final relation is: This relation shows that in the blurred locations, identified by the blurring detector, u is restored by the forth order of differentiation in the blurred areas.

Conclusion
We presented a noise removal method using a local switching noise estimator in an energy functional minimizing process. We showed that in addition to using Laplacian in edge detection tasks, it can be used for noise removal applications. Smoothness can be easily controlled just by the number of iterations. We compared the  method with some classic methods like ROF and Wiener filters based on TV value and also with the state-of-the-art BM3D. The result of denoising is totally comparable with ROF model. The computation time is equal to ROF model with 10 iterations. Based on experimental results we concluded that in relation to classic filters like ROF the method appears to be easy, fast and applicable for many noisy images. We analyzed that the technique can be applied for other image processing applications like deblurring by defining the appropriate detector and weight functions. The main disadvantage of the method is its filtering action on the area of the image including low intensity edges. The method can be improved to represent better response by defining better noise detectors.

Appendix A: driving denoising equation
Here we drive the denoising equation by minimizing the energy functional. First we use Taylor series and relate u (x − Δx, y) to x component of the gradient of the image u x . We have: Then g x can be computed as: And for g x + Δx we can write: Using Taylor series similarly for y components, g y and g y + Δy are: Δy u yy ð27Þ Based on the above relations, we write: The minimizer u is found by ∂J ∂u ¼ 0 . Considering the following relations for simplicity: We can write:  Differentiating J with respect to u we have: ∂u H 2 Z 1 Z 2 þλ ∂H −g y g yþΔy ∂u H 1 Z 1 Z 2 −2λ g x þ g xþΔx À Á ∂ g x þ g xþΔx À Á ∂u H 1 H 2 Z 1 Z 2 −2λ g y þ g yþΔy ∂ g y þ g yþΔy ∂u First we compute ∂g x ∂u : Similarly for g x + Δx we have: Then we can write: And for y direction: ∂ g y þ g yþΔy Second we have: ∂ −g x g xþΔx À Á ∂u ¼ −2u xx −2u xy u x u y þ Δx 2 2 u xxx u x þ u xxy u y u xx ð38Þ ∂ −g y g yþΔy ∂u ¼ −2u yy −2u yx u y u x þ Δy 2 2 u yyy u y þ u yyx u x u yy Derivative of Heaviside Function appears the impulse function: where δ(.) is the impulse function. We Ignore its effect in the denoising process and assume λ is constant during the implementation. Therefore we have: ∂H −g x g xþΔx À Á ∂u ≈0 ; ∂H −g y g yþΔy ∂u ≈0 ; ∂λ ∂u ≈0 ð41Þ Figure 11 Decreasing behavior of f(u 0 ) during the denoising process.