Edge-Preserving Image Denoising Based on Lipschitz Estimation

: The information transmitted in the form of signals or images is often corrupted with noise. These noise elements can occur due to the relative motion, noisy channels, error in measurements, and environmental conditions (rain, fog, change in illumination, etc.) and result in the degradation of images acquired by a camera. In this paper, we address these issues, focusing mainly on the edges that correspond to the abrupt changes in the signal or images. Preserving these important structures, such as edges or transitions and textures, has signiﬁcant theoretical importance. These image structures are important, more speciﬁcally, for visual perception. The most signiﬁcant information about the structure of the image or type of the signal is often hidden inside these transitions. Therefore it is necessary to preserve them. This paper introduces a method to reduce noise and to preserve edges while performing Non-Destructive Testing (NDT). The method computes Lipschitz exponents of transitions to identify the level of discontinuity. Continuous wavelet transform-based multi-scale analysis highlights the modulus maxima of the respective transitions. Lipschitz values estimated from these maxima are used as a measure to preserve edges in the presence of noise. Experimental results show that the noisy data sample and smoothness-based heuristic approach in the spatial domain restored noise-free images while preserving edges.


Introduction
Several industrial imaging systems are widely employed in different applications of Non-Destructive Testing (NDT). The presence of noise elements or degradation added during the testing process can directly affect the performance of the evaluation. In general, denoising such data (in the form of a signal or image) involves removing these degradations while capturing and transmitting signals or images. These degradations are commonly categorized as blur or noise. Blurring, the most common type of degradation, involves bandwidth reduction due to a problem in the image formation. There can be several reasons for this imperfection. Most commonly, these problems occur due to the relative motion of the camera and the original scene or due to an out of focus optical system.
In addition to these blurring effects, the recorded image is corrupted by additional noise elements. These additional noise elements are introduced during the transmission via a noisy channel or errors during the measurement process and directly affect the performance of NDT. Therefore, it is of extreme importance to reduce noise and preserve the sharpness of edges without introducing blurring. Goyal et al. evaluated the most state-of-the-art image denoising methods in [1].
Methods proposed in the past are mainly categorized into spatial domain and transformbased techniques. The selection is dependent on the application and the nature of noise present in the image. Denoising in the spatial domain is performed using linear methods based on an averaging filter, e.g., Gabor proposed a model based on Gaussian smoothing [2], and edge detection-based denoising using anisotropic filtering [3][4][5]. Yaroslavsky and Manduchi used a neighborhood filtering method [6,7]. Image denoising in the spatial domain is complex and more time-consuming especially for real-time problems.
Some researchers worked in the frequency domain to denoise images, e.g., Wiener filters [6]. Chatterjee and Milanfar improved the performance of denoising by taking advantage of patch redundancy from the Weiner filter-based method [8]. The scope of these linear techniques is limited to stationary data, as the Fourier transform is not suitable to perform analysis on non-stationary and nonlinear data. As a result, denoising methods based on multi-scale denoising employing nonlinear operations in the transform domain came up as an alternative. These multiple scales provide sparse representation of the signal in the transform domain.
Several methods based on wavelet transform have been proposed in the past [9][10][11][12][13][14][15]. A brief survey of some of these methods is given by Buades et al. in [16,17]. Seo et al. presented a comparison and discussion on the kernel based methods [18]. The method proposed by Blu and Luisier minimized the mean square error estimate termed as -Stein's unbiased risk estimate (SURE) [19]. In another method, dual tree complex wavelets were transformed into ridgelet transforms by Chen and Kegl [20]. Starck et al. proposed image denoising based on the family of transforms-the ridgelet and curvelet transforms-as alternatives to the wavelet representation of image data [12].
Similarly, Tessens et al. used curvelet coefficients to distinguish two classes of coefficients, e.g., noise-free components, which they termed as "signal of interest," and other [10]. More recently, Pesquet et al. adopted a hybrid approach that combined frequency and multi-scale analysis [13]. They formulated the restoration problem as a nonlinear estimation problem leading to the minimization of a criterion derived from Stein's unbiased quadratic risk estimate. Alkinani et al. investigated patch-based denoising methods for additive noise reduction [21]. Recently, Naveed et al. proposed DWT-and DT-CWT-based image denoising methods and employed statistical goodness of fit tests on wavelet coefficients at multi-scales to identify noise at different levels [22].
Bnou proposed a method based on an unsupervised learning model. They presented adaptive dictionary learning-based denoising of approximation. The wavelet coefficients were denoised by using an adaptive dictionary learned over the set of extracted patches from the wavelet representation of the corrupted image [23]. Several other methods have also been proposed in the past, but there is limited study regarding the preservation of edges.
Recently, Paras and Vipin presented a brief evolution of the research on edge-preserving denoising methods in [24], and Pizurica presented an overview of image denoising algorithms ranging from wavelet shrinkage to patch-based non-local processing with the main focus on the suppression of additive Gaussian noise [25]. Pal et al. reviewed the benchmark edge-preserving smoothing algorithms and keeping the main focus on anisotropic diffusion and bilateral filtering [26].
Many other filtering algorithms exists but, as mentioned before, very few of them discussed the edges and the spurious oscillations or artifacts introduced after denoising. Gibbs phenomena explain these oscillations in the neighborhood of discontinuities. An edge separates two different areas with limited or no connections between them e.g., no common characteristics: texture, noise, etc. Figure 1a, presents an edge part of the House image. Figure 1b,c represents the same image with the addition of white Gaussian noise and denoised image [19]. It can be observed from the figure that, while denoising, some additional artifacts have been produced. This is because most linear or nonlinear methods use neighborhood pixels while denoising. It is, therefore, logical to separate edges and transitions from the rest of the background.
Most of the wavelet-based methods use Discrete Wavelet Transform (DWT) for denoising; however, DWT suffers from three main issues, lack of shift-invariant, poor directional, and lack of phase information. Stationary wavelet transform reduced the problem of partial translation in variance and considerably improved the denoising results; however, it suffers from the cost of very high redundancy and is ultimately computationally expensive. Many algorithms have been proposed to solve the shortcomings of DWT by using different forms of Complex Wavelet Transforms (CWT) [27].
CWT is nearly translation invariant with directional selectivity but at the cost of high redundancy, [27]. Similarly, through continuous wavelet transform analysis, a set of wavelet coefficients is obtained, indicating how close the signal is to a particular basis function. Since the continuous wavelet transform behaves like orthonormal basis decomposition, it can be shown that it is isometric, i.e., it preserves energy [28]. Hence, the signal can be recovered from its transform.
The wavelet transform in its continuous form can accurately represent minor variations/edges present in signal. In most of the denoising methods, these minor variations that correspond to the edges of different levels in images are not considered while denoising. In this work, we classify this separation as classes for edges and background using continuous wavelet transform. This classification is based on the information obtained from the Lipschitz estimation.
The rest of the paper is organized as follows. Section 2 introduces the method and the principle of method. Section 3 gives a detailed explanation of the multi-scale analysis based edge detection. In Section 4, we present the Lipschitz regularity in the context of images, followed by the explanation of the reconstruction method. Section 5 explains the obtained results and evaluates the performance of the proposed method, and Section 6 finally concludes the paper.

Method
We identified transitions and restored noise-free image using heuristic approach. Initially, all significant edges where the intensity changes abruptly were identified by using the modulus maxima approach using multiscale analysis. Based on the detection of edges from modulus maximas, we classified the information from the image into the edge and background. Lipschitz exponents based reconstruction resulted in noise-free images. The edges where the intensity changes abruptly preserved their Lipschitz estimation even with the addition of noise.

Principle of the Method
We assume that the image is corrupted by different types of noise. In this work, we focused on white Gaussian noise with a zero mean. The given image is specified in the model: f is the noise-free signal, uniformly sampled (e.g., an image) and ij is the white Gaussian noise with 0 mean and σ 2 variance. In the present work, the given data will always be an n × n matrix with n = 2 N . The aim is to estimate the function F : ( f (x ij )) n i,j=1 with respect to an estimatorF such that: In order to estimate the function F, the method is proposed to preserve the energy of the strong edges by separating them from the background. The restoration based on data samples and smoothing results in noise-free images. The block diagram of the method is given in Figure 2.

Multi Scale Analysis Based Edge Detection
Edges correspond to the intensity changes in an image and separate two different regions by defining their boundary. In most denoising methods, while suppressing noise elements and the smoothing process, significant changes are introduced in the highfrequency elements which correspond to the edges. The main purpose of this paper is to preserve the sharpness of these high-intensity elements and, to do so, we employed wavelet transform-based multi-scale analysis. We decomposed the image into different scales and detected the local minima and maxima at each level. The multi scale was performed by smoothing the surface with the convolution kernel θ that is dilated, θ correspond to the Gaussian function [27]. The wavelet transform components are proportional to the coordinates of the gradient vector of f smoothed by θ 2 j : If I(x, y) corresponds to the intensity image after wavelet transformation, then the discrete version can be written as [27]: G i (i, j) and G j (i, j) correspond to the gradients in the horizontal and vertical directions, respectively. The gradient image in both directions (horizontal and vertical) of the House image is given in Figure 3a,b respectively.
Let A f (u, 2 j ) be the angle of the wavelet transform vector in the plane (x 1 , x 2 ) [27].

Local Maxima Extraction
The maxima across each scale were extracted from the gradient of the image and by following the direction of the phase angle. The comparison of gradients from the neighborhood pixels classified the current pixel as the local maxima or otherwise. The current pixel can only be the local maxima if it satisfies the following condition: where − −− → |G n+ | represents the gradient in the positive phase angle of the neighborhood pixel, − −− → |G n− | represents the gradient in the negative phase angle of the neighborhood pixel, and − − → |G c | corresponds to the gradient of the current pixel/point. Figure 4b illustrates the detected maxima of the image representing the edges of the House, the background, and noise. We employed Median Absolute Deviation (MAD)-based thresholding to eliminate noise and background maxima from the main structure of the image. Figure 4 illustrates the corresponding maxima before and after applying thresholding.
In signal feature extraction, the selection of the mother wavelet function is directly affected by the characteristics of the signal/image. In this context, to extract the most useful information, the selection of wavelet function is of extreme importance. At the same time, while computing the modulus maxima, it is not guaranteed that the modulus maxima is located at a point belonging to a maxima line that propagates toward a finer line. Mallat highlighted that, when the scale decreases, W f (u, s) may have no more maxima in the neighborhood of u = u 0 . He proved that by using the diffusion equation that this would never be the case if the wavelet function is a derivative of Gaussian. Hence, the derivative of a Gaussian function has an important property of continuous differentiability, which makes this function suitable for the analysis of most type of signals [27].

Lipschitz Based Image Restoration
In order to perform the restoration process, we computed Lipschitz exponents to identify strong edges. These exponents were estimated using modulus maxima lines. The modulus maxima curves are defined as the lines along which all points correspond to the modulus maxima. In the case of images, these maxima lines can be traced by following the neighborhood of local maxima in successive scales as shown in Figure 5. The boundary of some important structures can be traced by following the series of maxima that are normally distributed along the edge curve. Wavelet maxima that correspond to the common edge at different scales are chained together by following the direction of the tangent to create the maxima curve of the respective edge at different scale levels. Therefore, by comparing the phase of edge points with the phase of an adjacent neighborhood in the next scale, it can give the pixel that corresponds to the maxima curve. It was shown by Mallat that point-wise singularities can be computed by measuring the decay of the slope of the logarithmic function of the coefficient of magnitude along the logarithmic function of scales, and this is termed as the Lipschitz exponent.

Lipschitz Regularity in the Context of Images
The Lipschitz regularity of a two dimensional function (image) f is related to the asymptotic decay of |W 1 f (u, 2 j )| and |W 2 f (u, 2 j )| in the corresponding neighborhood, and this decay is dependent on the magnitude defined by M f (u, 2 j ). It was proven by Mallat that f is uniformly Lipschitz α inside the bounded domain of R 2 , if and only if there exist a constant A > 0 such that for all the u inside this domain and all scales 2 j [27]: However, we are concerned with the pointwise regularity in 2D. Therefore, by using the local maxima lines, we have computed the Lipschitz regularity by using Equation (9). In Figures 6 and 7 from the histogram of the Lipschitz estimation, it is estimated that, in the present case, the Lipschitz value (α) lies in a certain range.

Restoration
In order to obtain a noise-free image from the background, the method utilizes data samples and performs non linear functioning to estimate the best fit. We defined an adaptive restoration method for two classes: one for the edges and one for the background. The concept of the class definition is given in Figure 8a. The edge points belong to class 1 and rest of the image belong to class 2. While denoising, we centered the pixel on a 3 × 3 window.
The weighting factor depends on the adjacent neighborhood of the respective class. The method restores each pixel by taking into account extracted edge and background pixel based on their respective Lipschitz exponents. The method preserve edges and restores the background by utilizing noisy samples and their smoothness to estimate the best fit. Initially, we defined the deterministic function in terms of the smoothness and mean square error estimation as: where i, j = 1, . . . , n corresponds to the size of the image and k defines the iteration step.
is the mean square error estimation of the restored subset with the original signal of respective subset (y oi,j ) such that defines the smoothness of the reconstructed subset signal. In this context, each pixel will generate the equation according to the adjacent neighborhood (as explained in Figure 8) for smoothing: where A g for each pixel depends on if the current pixel corresponds to the background or an edge, and hence, based on that, the assigned value would be "1" or it would correspond to "0" as illustrated in Figure 8b. In order to simplify Equation (10), we consider the Taylor series expansion to find the value of λ k i and γ k i and for the given series, the minimum mean square error correspond to f (x + dx) = 0, therefore, by simplifying the Taylor series expansion: and we know that: As x → y k+1 i,j and f → C MSE,k i,j , the variables in Equation (15) can be replaced such that Thus, and we can simplify the equation to Similarly by using the Taylor series on the smoothing criteria, γ k i,j is as follows: In order to determine the best estimation of the deterministic function and minimum error between the restored pixel and the original pixel of either background or an edge, the algorithm works in the iterative mode. We assume the least smoothness of the respective pixels for initialization. Two optimization factors, w M for noisy samples and w S for smoothness, are introduced (in Equation (10) to investigate the most optimal compromise between the data samples containing noise and their smoothness such that Equation (10) will become: Parametric Estimation In order to estimate w M and w S , we can write Equation (20) in the discrete case by simplifying as, where y adj correspond to the adjacent neighbors, The value of w M or w S is dependent on whether the current pixel belongs to the edge or background but, in any case, should satisfy the stability conditions. The term that corresponds to the mean square error in Equation (20) is stable due to the presence of the original noisy pixel, and hence these factors were computed in terms of the smoothness criteria, and Equation (21) can be rewritten as: The Equation (22), should satisfy the following condition in order to ensure the stability of the signal: The stable range to define the constant w S term in numerical form between 0 to 1 64 is obtained by considering the above conditions. Hence, based on these findings, we define two constant factors in term of an equation: All pixels were restored individually taking into account the neighbor pixels of the same class. The complete algorithm is given in the Appendix A.

Results and Analysis
This section presents the performance evaluation of the proposed algorithm. The set of set of images used for experimentation consisted of standard test images, including House, Lena, Circuit, Pepper, and Hand, and were tested with various noise levels. The images were corrupted with the addition of Gaussian noise of zero mean and standard deviation ranging from 15 to 35.
The peak signal to noise ratio (PSNR) was employed as the measure of quantitative performance as explained in Section 2. We computed Lipschitz exponents to identify the edges or image structure. These exponents were estimated using the modulus maxima lines. We analyzed these estimations with the addition of Gaussian white noise (28 dB). In order to illustrate more detailed analysis, we show in Figure 9, a single row (125) of the image as a curve without and with the addition of noise. It can be seen from the results that sharp transitions preserved their Lipschitz estimation even in the presence of noise.
To illustrate the robustness of the estimation (α) in the case of natural images, we analyzed our test samples in the case of a constant affine transformation, as shown in Figure 10. The original image has been deformed by the 10 • of rotation. Lipschitz exponents were estimated, and, from the histogram, we observed that approximately 92% of the Lipschitz regularity values lay in between the same range. The regularity (α) was preserved even when a constant affine deformation was applied to an image. Therefore, we concluded from our findings that the Lipschitz exponents can be used as a significant measure to study natural images. At the same time, the problem of preserving edges in any linear or nonlinear denoising algorithm can be handled with this estimation. The transitions identified based on their Lipschitz exponents represent the main structure of the image; therefore, we preserved these points and performed smoothing on the rest of the image, as explained in the Section 4.  We defined an adaptive restoration method and performed data samples based nonlinear functioning to estimate the best fit. The method restores therest of the image by utilizing noisy sample and their smoothness. As an example, Figure 11 illustrates the focused edges of the image overlaid with the results of edge extraction. Noise level of 24.59 dB PSNR were added and by performing Lipschitz analysis and restoration process as explained in previous sections, the PSNR is improved to 29.14 dB. The sharpness of the edges were significantly preserved, and the smoothing process on the rest of the image results in increasing PSNR of the image.
The obtained results in Figure 11c shows that the noise level has reduced and preserved the image structures without introducing spurious oscillations. To illustrate the results of denoising more precisely, we show in Figures 12 and 13, a single row (200) of House and Lena image with noise (22.11 dB white Gaussian noise), plotted as a curve. It can be seen from the figure that the method has filtered most of the noise elements while preserving the edge points. The sharpness and magnitude of these edges were not altered during denoising and the noise elements were removed from the homogeneous regions. The results in Figure 13 also illustrate that the sharp transition or edges (near to step function or discontinuity), are equally preserved and hence the main information about the structure of the image remain unchanged. We performed the multi scale analysis based splitting and heuristic approach for the reconstruction of House, Lena, circuit, Pepper and Hand images. We generated noisy data from clean image by adding pseudo random numbers (Gaussian White noise) with different peak signal to noise ration (PSNR), the qualitative analysis of the method on both images is summarized in Table 1. We observed improvement of ≈5 at the lowest sigma and relatively better results at high variance of noise with the ≈10 improvement in PSNR in House and ≈8 in the case of Lena image. Figure 14 shows the corresponding correlation curves with different noise levels (sigma) on House image where cross and circle represent the true values computed by the method.
We observed distant correlation with input at high sigma which gradually reduces with decrease in sigma. We generated noisy data from House image by adding pseudo random numbers (Gaussian White noise) resulting in peak signal to noise ratio (PSNR) of approximately 19 dB. Figure 15b,c presents the image with the addition of Gaussian white noise and the denoising result with the proposed method. The obtained PSNR is 26 dB with Lipschitz based method and also the visual evaluation emphasizes the proposed method has not introduced any spurious edges as highlighted in Figure 1.     Similarly, Figure 16b,c also presents the image with the addition of Gaussian white noise and the denoising result with the proposed method on Lena image. The improved PSNR is 29.82 dB as compared to the input 22 dB. In Figure 17, the method improved the PSNR on the Circuit image to 26.12 dB. The results obtained on the Hand and Pepper images are also shown in Figures 18 and 19. The proposed Lipschitz estimation-based method significantly improved not only the statistical results but also the visual evaluation to emphasize that the proposed method has not introduced any spurious edges.
We compared the results on the Lena test images with other methods, including the SURE LET, Sure Shrink, VisuShrink, and Bivariate shrinkage functions (Table 2), as explained [29]. We noted that the results obtained at higher noise level were comparable and better than those obtained with the other methods from the literature. The obtained results were approximately similar to the SURE LET and Bivariate shrinkage functions and outperformed the SureShrink and VisuShrink methods.     Tables 3 and 4 present the qualitative analysis in terms of the Structural SIMilarity (SSIM) of the proposed method with the other denoised images obtained from the waveletbased denoising methods [22]. The SSIM can estimate the similarity index between two images and, therefore, can be viewed as a quality measure of one of the image degraded or altered with the other which is regarded as of perfect quality. We performed analysis on the 'Lena' and 'Peppers' images. It can be observed that the denoised images obtained better results than the Bishrink and NeighSure method and were competitive to the SURE LET denoising method in terms of SSIM.
We also analyzed the performance of the method using another criterion, which focused only around the edges. At first, we applied a canny edge detector as shown in Figure 20a to highlight strong edges. Then, we applied morphological dilation of 5 × 5 to generate a binary mask as shown in Figure 20b. The Mask represents the dilated binary image of edges, and the estimated error is: We estimated the error corresponding to the area of the edge points with the SURE LET denoising method with the proposed method based on the splitting of edges as shown in Figure 20, and we observed increased improvement in SNR with both methods with the difference of ≈1 dB. Table 3. Comparison of the proposed methods with the other wavelet-based image denoising methods on the 'Lena' image in terms of the structural similarity (SSIM) for a range of input noise levels σ = 10 to σ = 50 [22].

Conclusions
This paper introduced a method to remove the noise from the signals or images during the acquisition phase of theNDT process. The method preserved the edges and transitions that corresponded to the most significant information about the nature of the signals or images while denoising. The algorithm removed the noise in the homogeneous areas but preserved all structures like the edges or corners by computing the Lipschitz exponents estimated from the modulus maxima lines.
It is shown that the Lipschitz estimation of the anatural image lay within a certain range of values even in the presence of noise and also did not change significantly even in the presence of affine deformation. Based on the detection of edges from the modulus maximas, we classified each image into the edge and background pixels-the edges where the intensity changed abruptly preserved their Lipschitz estimation even with the addition of noise.
The statistical results illustrated that the proposed method had improved PSNR and maintained the sharpness of the edges without introducing any artifacts. The Lipschitz estimation can highlight the discontinuities present in the images, and in more practical applications (e.g., NDT, medical applications), it is of extreme importance to localize them. The results are encouraging to further improve this method by working on multi-scale basis and employing practical examples e.g., surface inspection.