Adaptive High-order Variation De-noising Method for Edge Detection with Wavelet Coefficients

Abstract

This study discusses the high-order diffusion method in the wavelet domain. It aims to improve the edge protection capability of the high-order diffusion method using wavelet coefficients that can reflect image information. During the first step of the proposed diffusion method, the wavelet packet decomposition is a more refined decomposition method that can extract the texture and structure information of the image at different resolution levels. The high-frequency wavelet coefficients are then used to construct the edge detection function. Subsequently, because accurate wavelet coefficients can more accurately reflect the edges and details of the image information, by introducing the idea of state weight, a scheme for recovering wavelet coefficients is proposed. Finally, the edge detection function is constructed by the module of the wavelet coefficients to guide high-order diffusion, the denoised image is obtained. The experimental results showed that the method presented in this study improves the denoising ability of the high-order diffusion model, and the edge protection index (SSIM) outperforms the main methods, including the block matching and 3D collaborative filtering (BM3D) and the deep learning-based image processing methods. For images with rich textural details, the present method improves the clarity of the obtained images and the completeness of the edges, demonstrating its advantages in denoising and edge protection.

1. Introduction

Images, a crucial component of information sources, are often disturbed by noise during formation, acquisition, and transmission. The existence of noise not only affects the visual effect of the image but also poses some challenges for image compression, edge detection, and image fusion. Therefore, image denoising is a crucial step in image processing and has significant research implications. However, edge blurring inevitably occurs while removing noise because there is a contradiction between suppressing noise and maintaining image features. Human vision is sensitive to the high-frequency components of an image (details and edges), and the important information of an image is primarily found in the edges and contours. Therefore, several studies have focused on the algorithms that can maintain small-scale and edge features while also denoising.

Recent restoration techniques based on variation and partial differential equation (PDE) methods have become popular for image denoising. For instance, the Perona-Malik model (PM) [1], total variation model [2], adaptive total variation method [3], and newly anisotropic variation models [4-11] are some of the models. These schemes have been demonstrated to be effective at removing noise and maintaining edges; however, they produce so-called blocky effects in the piecewise smooth transition region in denoised images [12]. It makes some detailed areas in the image blurry and difficult to identify. Researchers have studied high-order PDE filter techniques to address these shortcomings. For example, the You and Kaveh (YK) [13], Lysaker, Lundervold, and Tai (LLT) [14], high-order regularization [15,16] models, and fourth-order anisotropic diffusion [17-22] are some of the techniques. These models demonstrated that the fourth-order PDEs are superior to the second-order PDEs in some aspects [15]. Although these fourth-order models can reduce the blocky effects and have the advantage of maintaining smoothness in flat areas, they lose clarity in critical geometric structures such as image edges. To overcome the blocky effects and avoid creating blurred edges, hybrid denoising methods combined with variation have been proposed. For example, machine learning-based PDE models [23, 24] lack interpretability and stability, which also produce artificial effects; combining TV and fourth-order PDE filter schemes [25-28], and a high-order non-convex variation method [29]. The edge protection effect and time complexity of these methods must be improved. In many nonlinear diffusion models, image edges are detected by the diffusion function, which is composed of the module of the image gradient or the module of the gradient of the image convolved with a Gaussian function. When processing images with rich textures, the gradient is easily disturbed by noise, and essential information is easily lost. Wavelet analysis has good time-frequency localization properties, which can present the texture and structure information of the image at different resolution levels and can preserve the details of the image while realizing the separation of signal and noise [30-36]. The module of the wavelet coefficients, which reflects changes in image grayscale, is particularly good at resisting noise [37]. Wavelet coefficients can be used to calculate the diffusivity on a single scale or several scales, which accurately reflect the image information at different scales.

In this study, we construct the edge detection function using the module of the wavelet coefficients rather than the module of the gradient to reduce the impact of noise on diffusion coefficients. Conversely, we use the above edge detection function to guide the high-order diffusion through theoretical analysis, which can better maintain the edge information of the denoised image and improve the diffusion filtering effect. Owing to the influence of noise, some wavelet coefficients obtained by the wavelet transform cannot reflect the information of the image itself. We use the concept of state weight to obtain the denoised wavelet coefficients, which are more accurate wavelet coefficients. We normalize the coefficients to obtain the corresponding continuous state weights. The critical step of the method is to apply an anisotropic filter to the state weights, which can ensure the accuracy of the weights. By applying the accurate state weights to the corresponding wavelet coefficients, we can obtain the wavelet coefficients that accurately reflect the image information. This method provides a strong guarantee for the follow-up work.

The rest of this paper is structured as follows. In Section 2, we briefly review the variational diffusion methods and the wavelet shrinkage. Section 3 describes a fourth-order variation denoising method based on the edge detection function constructed by wavelet coefficients. The experimental results for demonstrating the efficiency of the proposed method are provided in Section 4. Finally, the conclusions are presented in Section 5.

2. Mathematical Framework

2.1 Variation and PDE Diffusion Method

Nonlinear diffusion has a significant influence on image processing. It approaches the actual signal in the image with a piecewise continuous function through the PDE. Unlike general linear diffusion, nonlinear diffusion preserves the edges and details while also smoothing the appearance. Numerous studies have carried out extensive research based on this property. Nonlinear PDE [1] and variation-based nonlinear models [2] have been widely used for image denoising. These studies have paved the way for many new fields in the theory and application of image processing. One of the classical denoising models is the total variation (TV) minimization model proposed by Rudin et al. [2], which can be presented as follows:

\(\begin{aligned}\min \int_{\Omega}|\nabla u|+\frac{\lambda}{2} \int_{\Omega}\left(u-u_{0}\right)^{2} d x\end{aligned}\)       (1)

whereΩ denotes an open subset with Lipschitz boundary, ∇ denotes the gradient operator. u and u₀ represent the original image and the observed image with noise, respectively. The first part in (1) is called the regularisation term, which is the TV norm of u . The regularisation term is crucial in the TV model because it causes the image to have discontinuous parts but no oscillations. Thus, noise can be removed. The second part in (1) is the approximation term, which controls the difference between u and u₀ . λ > 0 is the regularisation parameter, which serves to balance the regularisation and approximation terms. The corresponding Euler-Lagrange equation is expressed as follows:

\(\begin{aligned} \frac{\partial u}{\partial t}= div (\left(\frac{\nabla u}{|\nabla u|}\right)-\lambda\left(u-u_{0}\right)\end{aligned}\)       (2)

To improve the edge detection effect, Strong and Chan [3] introduced an adaptive TV variational model as follows:

\(\begin{aligned}\min \int_{\Omega} g(x)|\nabla u|+\frac{\lambda}{2} \int_{\Omega}\left(u-u_{0}\right)^{2} d x\end{aligned}\)       (3)

where g(x) denotes an edge-stopping function that can adaptively control the degree of diffusion. In [3], g(x) is defined as follows:

\(\begin{aligned}g(x)=\frac{1}{1+K\left|\nabla G_{\sigma}^{*} u_{0}\right|^{2}}\end{aligned}\)       (4)

where G_σ denotes the Gaussian filter, σ denotes the standard deviation of the noise, and K denotes a threshold parameter. |∇G_σ^*u₀| has a larger value near the edges, and the corresponding g(x) receives a smaller value, where the adaptive TV model implements a weaker smoothing to protect the edge information. These schemes have been demonstrated to be effective at removing noise while preserving the edges. Unfortunately, the so-called blocky effects are generated in the smooth region of the recovered image because the image gradient magnitude is used as a regularisation term. To overcome the blocky effects, high-order PDEs (typically, fourth-order PDEs) have been introduced into image restoration [13-22]. You and Kaveh [13] proposed a family of fourth-order PDEs (the Y-K model). Lysaker, Lundervold, and Tai [14] proposed another classical fourth-order model called the LLT model. The basic concept of the LLT model is described as follows:

\(\begin{aligned}\min \int_{\Omega}\left|\nabla^{2} u\right|+\frac{\lambda}{2} \int_{\Omega}\left(u-u_{0}\right)^{2} d x\end{aligned}\)       (5)

where ∇² denotes the Laplacian operator and \(\begin{aligned}\left|\nabla^{2}\right|=\sqrt{u_{x x}^{2}+u_{x y}^{2}+u_{y x}^{2}+u_{y y}^{2}}\end{aligned}\). The LLT model has the advantage of effectively suppressing the generation of blocky effects while protecting the smooth region using the Laplacian operator and effectively protect important texture information. However, the second derivative responses to isolated noise more strongly, which is not conducive to preserving the boundary. The LLT model’s weakness is that it has poor edge protection ability for the restored image.

2.2 Wavelet Transform and Denoising Method

The wavelet shrinkage method is widely used in image restoration technology. The primary purpose of the method is to accurately determine the wavelet coefficients corresponding to the noise. The wavelet shrinkage method works on the principle of performing the wavelet transform on the image to obtain the wavelet coefficients. The wavelet coefficient of the smaller amplitude is provided by the smaller noise energy and the larger amplitude is provided by the larger signal Then, select an appropriate threshold value to remove the noise by specifying the wavelet coefficients lower than the threshold value to achieve the denoising goal [30-32].

Wavelet transform is used to produce the wavelet coefficients of different sub-bands, allowing for the representation of the on different scales. In practical applications, the signals are usually processed by computers, and the obtained signals are all discrete. Therefore, the study of discrete wavelet transform (DWT) is practically significant. f denotes a one-dimensional signal, φ denotes a scaling function, and ψ denotes a wavelet function. Given the orthogonal wavelet, the discrete wavelet transform is expressed as follows:

\(\begin{aligned}f=\sum_{i \in Z}\left\langle f, \varphi_{i}^{n}\right\rangle \varphi_{i}^{n}+\sum_{j=-\infty}^{n} \sum_{i \in z}\left\langle f, \psi_{i}^{j}\right\rangle \psi_{i}^{j}\end{aligned}\)       (6)

where ψ_i ^j(x):=2^−j/2ψ(2^-jx - i) and f ∈ L²(R). If the test signal f contains Gaussian white noise, all wavelet coefficients 〈f, ψ_i^j〉 following the wavelet transform will also contain noise; therefore, we must determine the effective wavelet coefficients that constitute the signal. Therefore, to demonstrate how to restore the original signal using the effective wavelet coefficients, we provide the following three steps in accordance with the denoising concept of the wavelet shrinkage method:

1. Analysis: apply a wavelet transform to the noisy signal f , and based on (1) obtain wavelet coefficients d_i^j = 〈f, ψ_i^j〉 and scaling coefficients c_iⁿ = 〈f, φ_iⁿ〉 ;

2. Shrinkage: perform nonlinear threshold processing on wavelet coefficients using shrinkage function S_θ with threshold parameter θ , that is, \(\begin{aligned}\tilde{d}_{i}^{j}=S_{\theta}\left(d_{i}^{j}\right)=S_{\theta}\left(\left\langle f, \psi_{i}^{j}\right\rangle\right)\end{aligned}\)       (6)

3. Synthesis: reconstruct the scale and wavelet coefficients processed by the threshold to obtain the recovered signal u :

\(\begin{aligned}u:=\sum_{i \in Z}\left\langle f, \varphi_{i}^{n}\right\rangle \varphi_{i}^{n}+\sum_{j=-\infty}^{n} \sum_{i \in z} S_{\theta}\left(\left\langle f, \psi_{i}^{j}\right\rangle\right) \psi_{i}^{j}\end{aligned}\)       (7)

Choosing the appropriate threshold in the wavelet shrinkage method directly affects the results of signal noise removal. If the selected threshold is small, some noise is preserved as image information. The denoising effect is not ideal because the image still contains more noise after the process. Although more noise can be removed with a large threshold, some of the detail features in the image are also filtered out, which blurs the vision. Therefore, the threshold must be determined according to the magnitude of the noise. The hard and soft threshold functions are the most important shrinkage functions.

A. Hard shrinkage

\(\begin{aligned} \tilde{d}_{i}^{j}=S_{\theta}(d_{i}^{j}= \begin{cases} d_{i}^{j} & |d_{i}^{j}| \geq \theta \\ 0 & |d_{i}^{j}|<\theta \end{cases} \end{aligned}\)       (8)

B. Soft shrinkage

\(\begin{aligned} \tilde{d}_{i}^{j}=S_{\theta}(d_{i}^{j})= \begin{cases} sign(d_{i}^{j}(|d_{i}^{j}|-\theta) & |d_{i}^{j}| \geq \theta \\ 0 & |d_{i}^{j}|<\theta \end{cases} \end{aligned}\)       (9)

where \(\begin{aligned}\theta=\sigma \sqrt{2 \ln \mathrm{N}}\end{aligned}\), σ is the standard deviation of the noise, N is the scale or length of the signal.

3. High-order Variation Denoising Method for Edge Detection

3.1 High-order Variation Model for Edge Detection with Wavelet Coefficients

In many higher-order diffusion models, the diffusion function generally detects edges using the module of the image gradient or the module of the gradient of the image convolved with a Gaussian function. The diffusion function is unable to fully extract the image edge information corresponding to different scales because the gradient is easily disturbed by noise. The difficulty of constructing a single high-order model is a reasonable trade-off between eliminating the blocky effects and protecting the edges. The wavelet transform can present the texture and structure information of the image at different resolution levels. It can preserve the details of the image while achieving the separation of the signal and noise [32-36]. Wavelet coefficients on a single scale or several scales can be used to calculate the diffusivity, which accurately reflects the image information at different scales because Gaussian noise has less interference with wavelet coefficients. To reduce the impact of noise on diffusion coefficients, we construct the edge detection function using the module of the wavelet coefficients rather than the module of the image gradient. Using the above edge detection function to guide the high-order diffusion can better maintain the edge information of the denoised image. Thus, we can improve the effect of diffusion filtering.

According to Section 2.2, we can obtain the wavelet coefficients d_j,k by discrete wavelet transform of the noise image u(x,y) . j represents scale. Because the edge information and noise are mainly concentrated in the high-frequency sub-bands, we set the horizontal, vertical, and diagonal high- frequency wavelet coefficients to be d_j,H , d_j,V and d_j,D , respectively. Thus, the module of the wavelet coefficients is defined as follows:

Mu(x, y) = (d_j,H² + d_j,V² + d_j,D²)^1/2       (10)

Therefore, we represent the high-order variation regularization scheme for edge detection with the wavelet coefficients:

\(\begin{aligned}\min _{u \in W^{2,1}(\Omega) \cap L^{2}(\Omega)} E(u)=\min _{u \in W^{2,1}(\Omega) \cap L^{2}(\Omega)} \int_{\Omega} g(M u(x, y))\left|\nabla^{2} u\right|+\frac{\lambda}{2} \int_{\Omega}\left(u-u_{0}\right)^{2} d x\end{aligned}\)       (11)

where g(Mu) denotes the adaptive diffusion function, which is determined using in the following equation:

\(\begin{aligned}g(M u)=\frac{1}{1+(|M u| / k)^{2}}\end{aligned}\)       (12)

where k denotes a parameter to control the diffusion extension, which can be set in advance or change as a result of each iteration. Mu(x, y) is small in the flat region of the image, the diffusion function g (Mu) is frequently chosen such that g(Mu) → 1 as Mu → 0 , which results in fast diffusion and more smoothing. This is particularly effective for noise removal. Contrary, Mu(x, y) obtains the local maximum at the edge of the image, the diffusion function g(Mu) is chosen such that g(Mu) → 0 as Mu → ∞ , which results in slow or stopped diffusion and less smoothing. This allows details such as edges to be protected.

Because the image is noisy, which affects the feedback of the wavelet coefficients to the image edge information, we can preprocess the wavelet coefficients to accurately reflect the image edge information. In this study, we propose the wavelet coefficients diffusion algorithm based on the continuous state wavelet threshold, which has a better theoretical explanation and experimental results.

3.2 Wavelet Coefficient Pre-processing

We perform DWT on noisy images to obtain the wavelet coefficients d_j,k . We can determine the state of each wavelet coefficient based on the size of the wavelet coefficients d_j,k . When |d_j,k| ≥ λ is satisfied, the current state of the wavelet coefficient is one, and when |d_j,k| < λ is satisfied, the current state of the wavelet coefficient is zero. Thus, we obtain the state quantity v_j,k^λ corresponding to all the wavelet coefficients.

\(\begin{aligned} v_{j, k}^{\lambda}= \begin{cases} 1 & |d_{j, k}| \geq \lambda \\ 0 & |d_{j, k}|<\lambda \end{cases} \end{aligned}\)       (13)

where λ denotes the threshold. According to the state of the wavelet coefficients, the weights w_j,k^λ can be determined as follow:

\(\begin{aligned} w_{j, k}^{\lambda}= \begin{cases} 1 & v_{j, k}^{\lambda}=1 \\ 0 & v_{j, k}^{\lambda}=0 \end{cases} \end{aligned}\)       (14)

By applying weights to the wavelet coefficients, noisy images can be effectively denoised.

\(\begin{aligned} w_{j, k}^{\lambda}(d_{j, k}= \begin{cases} d_{j, k} & |d_{j, k}| \geq \lambda \\ 0 & |d_{j, k}|<\lambda \end{cases} \end{aligned}\)       (15)

We refer to the wavelet hard threshold method described above as the two-state weight denoising method. The two-state weighted denoising method can remove noise, but there are only two states for the wavelet coefficients, which belong to the discrete state. The discrete state acting on the wavelet coefficients produces a discontinuous function. Thus, it is likely to cause ringing and pseudo-Gibbs phenomena, which distort the visuals. Therefore, the continuous state weights wavelet denoising method is proposed. To obtain the continuous state v_j,k , we must normalise the high-frequency wavelet coefficients.

\(\begin{aligned}v_{j, k}=\frac{d_{j, k}-\min d_{j, k}}{\max d_{j, k}-\min d_{j, k}}\end{aligned}\)       (16)

Conversely, ν_j,k is obtained owing to normalising wavelet coefficients and has some structural properties. However, noise causes the state variables to be estimated inaccurately. The nonlinear diffusion equation can effectively eliminate noise if the structure of the state variables is preserved. Then, we obtain more accurate state variables. We use the TV diffusion method to diffuse the state quantity ν_j,k . According to (2), the x-direction and y-direction are equal steps and ∆x =∆y = h , ∆t = τ , ν_i,jⁿ is an approximation of ν(ih, jh, nτ). \(\begin{aligned}g_{i, j}=g\left(\frac{1}{|\nabla v|}\right)_{i, j}\end{aligned}\).

The difference equation is obtained as follows:

\(\begin{aligned}v_{i, j}^{k+1}=\frac{\tau}{h^{2}}\left[g_{i+1, j} v_{i+2, j}^{k}-\left(g_{i+1, j}+g_{i, j}\right) v_{i+1, j}^{k}+2 g_{i, j}+g_{i, j+1} v_{i, j+2}^{k}-\left(g_{i, j+1}+g_{i, j}\right) v_{i, j+1}^{k}\right]+v_{i, j}^{k}\end{aligned}\)       (17)

This effectively removes noise while producing a more accurate state variable \(\begin{aligned}\tilde{v}_{j, k}\end{aligned}\). The weight can be calculated according to the value of \(\begin{aligned}\tilde{v}_{j, k}\end{aligned}\) and the weight coefficient \(\begin{aligned}\tilde{w}_{j, k}\end{aligned}\) is determined by \(\begin{aligned}\tilde{v}_{j, k}\end{aligned}\). \(\begin{aligned}\tilde{w}_{j, k}\end{aligned}\) is the continuous state weight coefficient’s value.

\(\begin{aligned}\tilde{w}_{j, k}=\tilde{v}_{j, k}\end{aligned}\)       (18)

From the above analysis, the corresponding weights \(\begin{aligned}\tilde{w}_{j, k}\end{aligned}\) of each high-frequency sub-band for the wavelet decomposition layer are obtained. By applying the \(\begin{aligned}\tilde{w}_{j, k}\end{aligned}\) to the wavelet coefficients of corresponding high-frequency sub-bands, the wavelet coefficients of denoising in each layer is obtained and the structure of these wavelet coefficients is preserved. The above methods can be expressed as follows

\(\begin{aligned}\tilde{w}_{j, k}\left(d_{j, k}\right)=\tilde{v}_{j, k}\left(d_{j, k}\right)\end{aligned}\)       (19)

The processed wavelet and scale coefficients are reconstructed to obtain the denoising image.

The wavelet transform can only perform a time-frequency decomposition on the lowfrequency wavelet coefficients decomposed at the upper level, and can no longer decompose the high- frequency wavelet coefficients. In this manner, some of the image’s high-frequency information is lost while noise is filtered. Wavelet packet decomposition-based image processing techniques have increased recently [38, 39]. The wavelet packet transform can decompose both the low and high frequency parts, which has higher time-frequency resolution characteristics. Compared with the wavelet decomposition, the wavelet packet decomposition is a more refined decomposition method. The wavelet packet transform has a better effect on noise suppression when used with various resolutions if one intends to obtain the high-frequency wavelet coefficients of each frequency sub-band, particularly for noisy images. The two-dimensional wavelet packet decomposition is shown in Fig. 1.

Fig. 1. Schematic diagram showing a two-layer wavelet packet decomposition

We can process the image’s rows and columns with the one-dimensional wavelet transform to achieve the wavelet packet transform, in accordance with the separability of the scale function of the wavelet transform. According to the study of wavelet basis functions, Daubechies (dbN) wavelet has advantages in symmetry, vanishing moment, and orthogonality. The support length of the dbN wavelet is 2N-1, while that of a wavelet in most applications is generally between five and nine. Therefore, we choose db3 with N=3 as the wavelet base for the experiments.

We conduct experiments to demonstrate the superiority of the module of the wavelet coefficients in edge detection. Simultaneously, the experiment’s wavelet coefficients are pre-processed using the wavelet coefficient diffusion method proposed in this study. By comparing with other edge detection methods, we fully demonstrate the effectiveness of the proposed wavelet coefficient pre-processing edge detection method. We compare our method with the following methods, including the gradient detection method, Gaussian smoothed gradient detection method, and the module of the wavelet coefficients detection method. We use standard images ‘Toys’ and ‘House’ for experiments and provide detailed experimental results, as shown in Figs. 2 and 3. In the experiment, for the module of the wavelet coefficients detection method, we use db3 for three-layer wavelet packet decomposition and apply our proposed method to noisy images. In (17), we set the experiment’s iteration stopping parameter ε=10^-3, τ=0.5 , and h = 2 .

Fig. 2. Results of edge detection obtained with ‘Toys’ image using four methods: (a) original image, (b) noisy image, (c) gradient detection method, (d) Gaussian smoothed gradient detection method, (e) the module of the wavelet coefficient detection method, (f) our proposed method.

Fig. 3. Results of edge detection obtained with ‘House’ image using four methods: (a) original image, (b) noisy image, (c) gradient detection method, (d) Gaussian smoothed gradient detection method (e) the module of the wavelet coefficient detection method, (f) our proposed method.

We conduct the experiment using the original Toys and House images, as shown in Figs. 2(a) and 3(a). The Gaussian white noise with standard variances of 10 and 15 is added to the Toys and House images, respectively, as shown in Fig. 2(b) and Fig. 3(b). Figs. 2(c) and 3(c) show the edge detection results obtained using the gradient. The presence of considerable noise affects edge detection, and the edges of the Toys and House images are blurred in many places. Figs. 2(d) and 3(d) show gradient detection images smoothed by the Gaussian function with a variance of two. We observe that noise is suppressed, and image edge detection performance is improved (for example the body of the toy puppy). According to Figs. 2(e) and 3(e), the influence of noise on wavelet coefficients is less than that on the gradient. The module of the wavelet coefficients detection method is superior to the Gaussian smoothed gradient detection method in terms of suppressing noise. The effect can be seen from the surface of the table and the house walls. Our proposed wavelet coefficient diffusion algorithm is used to pre-process the wavelet coefficients, as shown in Figs. 2(f) and 3(f). Therefore, the diffused wavelet coefficients significantly reduce the noise, which can be seen from the toys, tabletop, and house walls. Additionally, the edges of the toy bear's clothes, the puppy's body, and the roof are clearer than the other methods in Figs. 2(f) and 3(f).

In summary, owing to the use of the wavelet packet transform for image decomposition, the image features can be extracted more precisely and the wavelet coefficients reflecting the image information can be obtained. The pre-treatment method proposed in this study can reduce the impact of noise on wavelet coefficients, and experimental results show improved noise suppression and edge protection, which lays a foundation for the subsequent steps.

3.3 Preliminaries

In this section, we provide the notions and definitions related to the proposed model and theoretically demonstrate that the proposed model has a unique minimal solution. According to the references [21, 25], we define the BV² space as follows:

Definition 3.3.1 Let Ω ⊂ Rⁿ be an open bounded subset with a Lipschitz boundary. Given a locally integrable function u ∈ L¹(Ω) , the BV² semi-norm of u is determined as follows:

\(\begin{aligned}\int_{\Omega}\left|D^{2} u\right|=\sup _{\varphi \in C_{c}^{2}\left(\Omega, R^{n \times n}\right)}\left\{\int_{\Omega} \sum_{i, j=1}^{n} u \partial_{j} \partial_{i} \varphi^{i j} d x \| \varphi(x) \mid \leq 1\right\}\end{aligned}\)       (20)

where φ(x) denotes the vector function and \(\begin{aligned}|\varphi(x)|=\sqrt{\sum_{i, j=1}^{n}\left(\varphi^{i j}\right)^{2}}\end{aligned}\).

remark with the norm ||u||_BV²(Ω) = ∫_Ω|D²u| + ||u||_L¹(Ω) , BV² is a Banach space.

Definition 3.3.2 Let Ω ⊂ Rⁿ be an open bounded subset with a Lipschitz boundary and u ∈ L¹(Ω) .

Let g(x) denote a nonnegative continuous function. The BV² semi-norm of u is determined as follows:

\(\begin{aligned}\int_{\Omega} g\left|D^{2} u\right|=\sup _{\varphi \in C_{c}^{2}\left(\Omega, R^{n \times n}\right)}\left\{\int \sum_{\Omega i, j=1}^{n} u \partial_{j} \partial_{i} \varphi^{i j} d x \| \varphi(x) \mid \leq g\right\}\end{aligned}\)       (21)

and the g − BV² norm is ||u||_g-BV²(Ω) = ∫_Ωg|D²u| + ||u||_L¹(Ω) .

Proposition 3.3.1 (lower semi-continuity) If {u_i}_i=1^∞ belongs to BV²(Ω) and \(\begin{aligned}\lim _{i \rightarrow \infty} u_{i}=u^{*}\end{aligned}\), u_i, u* ∈ L¹(Ω) then

∫_Ωg|D²u*| ≤ lim inf∫_Ωg|D²u_i| (22)

Proposition 3.3.2 (compactness) If u* belongs to BV² (Ω) and g(x) is a nonnegative continuous function, then there is a minimal sequence \(\begin{aligned}\left\{\left\{u_{i}\right\}_{i=1}^{\infty} \mid \lim _{i \rightarrow \infty}\left\|u_{i}-u^{*}\right\|_{L_{(\Omega)}^{1}}=0, u_{i}, u_{i} \in B V^{2}(\Omega)\right\}\end{aligned}\), which satisfies the following

\(\begin{aligned}\lim _{i \rightarrow \infty} \int_{\Omega} g\left|D^{2} u_{i}\right|=\int_{\Omega} g\left|D^{2} u^{*}\right|\end{aligned}\)       (23)

Proposition 3.3.3 If u ∈ W^2,1(Ω) ⊂ BV²(Ω) , then ∫_Ω|D²u| = ∫_Ω|∇²u| .

According to the above-mentioned theoretical knowledge, we can demonstrate the existence and uniqueness of the solution of optimisation energy (11).

Theorem 3.1 The optimisation energy (11) has a unique minimiser in W^2,1(Ω) ⋂ L²(Ω).

Proof The weighted function g(Mu) is a nonnegative continuous function. Obviously, the energy function E(u) in (11) is strictly convex and coercive, while satisfying the low semi-continuity in W^2,1(Ω) ⋂ L²(Ω) . Therefore, from the convex analysis, the minimization problem (11) has a unique solution to the standard arguments.

3.4 Computational Methods

At present, the Euler-Lagrange equation serves as the primary solution to the minimisation problem (11). The Euler equation corresponding to (11) is expressed as follows:

\(\begin{aligned}\left(\nabla\left(\frac{g \nabla u_{x}}{\left|\nabla^{2} u\right|}\right)\right)_{x}+\left(\nabla\left(\frac{g \nabla u_{y}}{\left|\nabla^{2} u\right|}\right)\right)_{y}+\lambda\left(u-u_{0}\right)=0\end{aligned}\)       (24)

under the following boundary conditions:

\(\begin{aligned}\nabla u_{x} \cdot N=0 \quad \nabla u_{y} \cdot N=0 \quad \nabla\left(\frac{g \nabla u_{x}}{\left|\nabla^{2} u\right|}\right) \cdot n_{1}=0 \quad \nabla\left(\frac{g \nabla u_{y}}{\left|\nabla^{2} u\right|}\right) \cdot n_{2}=0\end{aligned}\)       (25)

where N = (n₁, n₂) denotes a unit outward normal vector of ∂Ω . Equation (24) can be solved by applying the steepest descent method to (26)

\(\begin{aligned}\frac{\partial u}{\partial t}=-\left(g \frac{u_{x x}}{\left|\nabla^{2} u\right|}\right)_{x x}-\left(g \frac{u_{x y}}{\left|\nabla^{2} u\right|}\right)_{y x}-\left(g \frac{u_{y x}}{\left|\nabla^{2} u\right|}\right)_{x y}-\left(g \frac{u_{y y}}{\left|\nabla^{2} u\right|}\right)_{y y}-\lambda\left(u-u_{0}\right)\end{aligned}\)       (26)

with the boundary conditions (25) and the initial condition u|_t=0 = u₀ on Ω . To obtain the numerical solution of (26), we adopt the following finite difference scheme [2, 14]. The time step and the space step are represented by τ and h , respectively. Suppose h = 1 and u_i,jⁿ is an approximation of u(i, j, nτ), then

D_x^±(u_i,jⁿ) = ±(u_i±1,jⁿ - u_i,jⁿ) D_y^±(u_i,jⁿ) = ±(u_i,j±1ⁿ - u_i,jⁿ)       (27)

D_xx(u_i,jⁿ) = D_x⁺(u_i,jⁿ) - D_x⁺(u_i-1,jⁿ) D_xy^±(u_i,jⁿ) = ±(D_x^±(u_i,j±1ⁿ) - D_x^±(u_i,jⁿ))      (28)

D_yx^±(u_i,jⁿ) = ±(D_x^±(u_i±1,jⁿ) - D_x^±(u_i,jⁿ)) D_yy(u_i,jⁿ) = D_x⁺(u_i,jⁿ) - D_x⁺(u_i,j-1ⁿ)       (29)

\(\begin{aligned}\left|D^{2}\left(u_{i, j}^{n}\right)\right|=\sqrt{\left(D_{x x}\left(u_{i, j}^{n}\right)\right)^{2}+\left(D_{x y}^{+}\left(u_{i, j}^{n}\right)\right)^{2}+\left(D_{y x}^{+}\left(u_{i, j}^{n}\right)\right)^{2}+\left(D_{y y}\left(u_{i, j}^{n}\right)\right)^{2}+\delta}\end{aligned}\)       (30)

where δ > 0 is a small parameter. Thus, by performing (n + 1) iterations on (26), the corresponding difference equation is characterized by the following:

\(\begin{aligned} u_{i, j}^{n+1}= & u_{i, j}^{n}-\tau\left[D_{x x}\left(g \frac{D_{x x} u_{i, j}^{n}}{\left|D^{2} u_{i, j}^{n}\right|}\right)+D_{y x}^{-}\left(g \frac{D_{x y}^{+} u_{i, j}^{n}}{\left|D^{2} u_{i, j}^{n}\right|}\right)+D_{x y}^{+}\left(g \frac{D_{y x}^{-} u_{i, j}^{n}}{\left|D^{2} u_{i, j}^{n}\right|}\right)\right. \\ & \left.+D_{y y}\left(g \frac{D_{y y} u_{i, j}^{n}}{\left|D^{2} u_{i, j}^{n}\right|}\right)\right]-\tau \lambda\left(u_{i, j}^{n}-u_{0}^{n}\right)\end{aligned}\)       (31)

4. Numerical Experiments

In this section, we consider subjective and objective image quality assessments to evaluate the performance of the proposed algorithm. We further demonstrate the effectiveness of our method by comparing it with the experimental results of four well-established denoising methods, including Hajiaboli proposed anisotropic fourth-order filter method [17], machine learning-based PDE models (LPDE) [24], hybrid high-order anisotropic diffusion model (SFAD) [28], and Block Matching and 3D collaborative filtering (BM3D) [40]. We select the images shown in Fig. 4 as the test images, which have a rich texture and edge details. The noise level affects the degree of smoothing, in which the stronger the noise, the greater the degree of smoothing. All test images are added with Gaussian white noise with standard deviations of 10, 20, 30, and 40 for denoising experiments. In our proposed method, we use db3 to perform a three-layer wavelet packet decomposition for the module of the wavelet coefficient detection method. To measure the quality of the numerical results of the five methods, we introduce the three concepts: the mean squared error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity index measure (SSIM). The MSE and PSNR are defined as follows:

Fig. 4. Testing images. The first row: Clock, Peppers, Lena, Toys; the second row: Cameraman, CoCo, Mandrill, House

\(\begin{aligned}\operatorname{MSE}(u, \bar{u})=\frac{1}{M \times N} \sum_{i}^{M} \sum_{j}^{N}\left(u_{i, j}-\bar{u}_{i, j}\right)^{2}\end{aligned}\)       (32)

\(\begin{aligned}P S N R=10 \log _{10}\left(\frac{255^{2}}{M S E}\right)\end{aligned}\)       (33)

where u and \(\begin{aligned}\bar{u}\end{aligned}\) denote the restored image and the observed image, respectively. We use the SSIM as an evaluation index defined as follows:

\(\begin{aligned}\operatorname{SSIM}(u, \bar{u})=\frac{\left(2 \mu_{u} \mu_{\bar{u}}+c_{1}\right)\left(2 \sigma_{u \bar{u}}+c_{2}\right)}{\left(\mu_{u}^{2}+\mu_{\bar{u}}^{2}+c_{1}\right)\left(\sigma_{u}^{2}+\sigma_{\bar{u}}^{2}+c_{2}\right)}\end{aligned}\)       (34)

Here, µ_u denotes the mean, σ_u² denotes the variance for the image u , \(\begin{aligned}\sigma_{u \bar{u}}\end{aligned}\) denotes the covariance of u and \(\begin{aligned}\bar{u}\end{aligned}\), and c₁ and c₂ are two constants. In this study, the stopping condition required for the experiments is expressed as follows:

\(\begin{aligned}\frac{\left\|u^{n+1}-u^{n}\right\|_{2}^{2}}{\left\|u^{n+1}\right\|_{2}^{2}} \leq \varepsilon\end{aligned}\)       (35)

where uⁿ and uⁿ⁺¹ denote denoising results at nth and (n + 1)th iteration, respectively, and ε denotes a given positive number. We set ε=10^-3 in the experiments.

We select the images from Fig. 4 to demonstrate our experiment and provide the detailed experimental results. The time step is 0.05 seconds in the experiment. The parameter values in the experiment are selected according to the related literature content. In our experiment, λ = 0.1 . Then, we add various levels of additive Gaussian noise to the original images. According to the following figures and tables, we compare and analyse several image denoising effects in terms of visual aspects, evaluation indicators.

4.1 Visual Comparison

We select the images ‘Clock’, ‘Lena’, ‘House’, ‘Toys’, and ‘finger’ for experimental display. The experimental results are shown in Figs. 5-9 and Table 1. To see the details of the experimental results more clearly, we select a portion of each practical picture to enlarge to better reflect the experimental effect. Figs. 5(a), 6(a), and 7(a) show the addition of Gaussian noise with standard deviations of 10, 20, and 30, respectively.

Fig. 5. Results of denoising obtained with ‘Clock’ image using five methods: (a) noisy image ( σ = 10 ), (b) Hajiaboli method, (c) LPDE model, (d) SFAD model, (e) BM3D model, (f) our model.

Table 1. Comparison of PSNR (dB) and SSIM for different models on four levels of Gaussian noise

The method proposed by Hajiaboli overcomes the blocky effects while reducing edge blurring. Nevertheless, there is an uneven area of the recovered image that is unevenly smooth, as shown in Fig. 5(b) (the edge of the number 12) and Fig. 7(b) (the edge of the eaves).

The LPDE method can preserve some details during the denoising process, but the results produced by this method still have considerable noise and blurry image edges. The decorations on the hat and lines of the eaves in Figs. 6(c) and7(c) reveal the image’s blurred edges.

Fig. 6. Results of denoising obtained with ‘Lena’ image using five methods: (a) noisy image ( σ = 20 ), (b) Hajiaboli method, (c) LPDE model, (d) SFAD model, (e) BM3D model, (f) our model.

The restored images obtained by the SFAD method are less noisy compared with those obtained by Hajiaboli. However, the edge-blurring phenomenon is still present in the recovery image (Figs. 5(d), 6(d), and 7(d)).

The BM3D method can effectively remove noise, but it is prone to produce an excessive smoothing phenomenon, Figs. 5(e) and 6(e) show that the enlarged numbers on the clock and feathers on Lena’s hat have smoothed, and Fig. 7(e) shows blurred eave lines. Thus, the BM3D method is not satisfactory in preserving edges.

Fig. 7. Results of denoising obtained with ‘House’ image using five methods: (a) noisy image ( σ = 30 ), Hajiaboli method, (c) LPDE model, (d) SFAD model, (e) BM3D model, (f) our model.

Our proposed method reduces the uneven effects in the smooth area of the image and prevents edge blurring (Figs. 5(f), 6(f), and 7(f)). Meanwhile, from the visual effect, the images recovered by our proposed method are clear, and the noise-suppression effect is obvious. Compared with other methods, the proposed method yields results that are smoother in the uniform regions and more definite in the edge regions.

This study also conducts experiments on images with rich texture details. We use the standard images ‘Toys’ and ‘fingers’ for our investigations and provide detailed experimental results. We add Gaussian noise with a standard deviation of 30 to the Toys image, as shown in Fig. 8(a). However, because the finger image already contains noise at the time of acquisition, we do not add noise to the finger image, as shown in Fig. 9(a).

Fig. 8. Results of denoising obtained with ‘Toys’ image using five methods: (a) noisy image ( σ = 30 ), (b) Hajiaboli method, (c) LPDE model, (d) SFAD model, (e) BM3D model, (f) our model.

Fig. 9. Results of denoising obtained with ‘fingers’ image using five methods: (a) original image, (b) Hajiaboli method, (c) LPDE model, (d) SFAD model, (e) BM3D model, and (f) our model.

Note that the bear clothes have rich texture information in the Toys image. The texture information recovered by the LPDE and SFAD methods contain more noise, and the texture clarity is not high, according to the comparison of the restored enlarged images. The BM3D method has a good denoising effect, but the texture of the image is not effectively recovered. Both the Hajiaboli method and our proposed method significantly reduce noise in the restored texture images, but the Hajiaboli method results in blurred texture details, as shown in Fig. 8(b). However, the image recovered using our proposed method has clear texture information and effectively protects image edges, such as the edge of the toy dog's body, as shown in Fig. 9(f).

The texture information recovered by the LPDE method still contains more noise in finger images. The SFAD and the BM3D methods retrieve texture information with blurred edges, which affects the visual effect. Both the Hajiaboli and our proposed methods retrieve clear texture images. However, the Hajiaboli method provides poor protection of the texture edges, with multiple discontinuities appearing, as shown in Fig. 9(b). Compared with other methods, our proposed method significantly improves the clarity of texture information and the integrity of texture image edges.

4.2 Evaluation indicator comparison

Although the Hajiaboli method, LPDE model, SFAD diffusion, BM3D method, and our proposed method all produce results that are comparable in terms of their visual effects, the given data can also objectively demonstrate the superiority of our proposed method in terms of recovery effects and edge protection. Table 1 compares the performances of the five denoising algorithms based on the quality of the denoised images. Note that the BM3D method and our proposed method have higher PSNRs than the other three denoising methods at the same noise level. When the noise level is low, the BM3D method can produce some images with the highest PSNR values, but as the noise level increases, the PSNR value produced by our proposed method is better than that of the BM3D method. However, compared with other diffusion methods, the LPDE method only obtains a lower SSIM and its edge protection capability is poor. Owing to the overly smooth image processing effect, the edge protection ability of the BM3D method is inferior to that of the SFAD method and the proposed method.

Objectively, our proposed method balances noise removal and edge protection to achieve higher PSNRs and SSIMs for different input noise variances.

Deep learning has recently performed well in the field of image denoising. This method has limitations even though it can produce a better image denoising effect. First, the method is only applicable to the learned images or image features included in the learned images. Second, the method requires many training image samples and a lengthy training period. This study investigates a non-machine learning image denoising method. Because it has a short running time while ensuring the denoising effect, it is more suitable image restoration method for real-time requirements. Recently, popular deep learning-based denoising algorithms include CDNet [41], DRCNN [42] and DCNN [43]. We conduct experiments on these algorithms, and compare them with our algorithm using the PSNR and SSIM indices. The image sets used for the experiments are set14 and Aliasing15 [44]. The average PSNR and SSIM values for the four methods under different image sets are shown in Tables 2 and 3.

Table 2. PSNR (dB) and SSIM comparisons for different levels of Gaussian noise

Table 3. PSNR (dB) and SSIM comparisons for different levels of Gaussian noise

According to Table 2, for the data set14 with image features learned, the average PSNR values obtained by the DRCNN method are higher than other denoising methods. By contrast, the average PSNR values obtained by the deep learning methods are all moderately higher than the average PSNR values obtained using our method, however, none is a significant advance. For SSIM values, the proposed method outperforms the DCNN, CDNet, and DRCNN, which also demonstrates the effectiveness of our method in terms of texture and edge protection.

protection. In Table 3, for the pseudo-image set Aliasing15, the average PSNR and SSIM values obtained using our proposed method are better than those obtained by the deep learning methods because the DCNN, CDNet, and DRCNN have not learned the features of the image set. Overall, our proposed method outperforms deep learning methods in terms of edge protection.

5. Conclusions

In this study, an image denoising method that combines the advantages of wavelet shrinkage and fourth-order anisotropic diffusion is presented from the perspective of edge protection. We obtain wavelet coefficients that can accurately reflect image information by introducing the concept of the state weight of wavelet coefficients. Then, rather than using the module of the gradient, we construct the edge detection function using the module of the wavelet coefficients. On this basis, we present a high-order diffusion model that can consider both denoising and edge protection.

The proposed hybrid framework with an edge-preserving strategy uses both wavelet shrinkage and fourth-order anisotropic diffusions. Our proposed algorithm can effectively remove noise and prevent the unevenness in smooth areas, while satisfactorily maintaining image edges compared with the other reference methods. Moreover, for untrained datasets, our proposed algorithm outperforms the deep learning-based image processing method in terms of PSNR and SSIM values. This study does not analyse and discuss more random noise. The next challenge is to determine how to improve the model so that it produces a good experimental effect under different noises.