Entropy solutions for an adaptive fourth-order nonlinear degenerate problem for noise removal

1.   Introduction

Noise is the most difficult task in the field of image processing and computer vision. In this work, we focus on removing additive Gaussian noise. The problem is formulated mathematically as: let $u(x, y)$ be a digital image and $u_0(x, y)$ be its observation with random noise $\eta(x, y)$. For $(x, y)\in \Omega$ s.t.

The noise level is approximately known

The goal of denoising is to filter out high frequency signals, while preserving the important features of the image such as edges. Therefore, we search for an image processing model which removes noise and offers better handling of edges.

During the last two decades, the method of partial differential equations (PDEs) for image processing has become a major research topic. A classical PDEs model named ROF (Rudin, Osher, Fatemi) model which was based on the total variation (TV), was first introduced by Rudin et al. ^[18]. The idea in the ROF model is to minimize the the total variation of the image $u$,

In fact, one of the main advantages of using ROF model for image restoration is that the discontinuities are allowed. However, the main drawback of denoising models based on the TV is that they tend to yield piecewise constant images, a phenomenon known as staircase effects. Strong ^[19], introduced an adaptive TV functional ($\text{TV}_{\alpha}$)

for spatially adaptive image restoration. The function $\alpha(x)$ is an edge detector to control the diffusion. The main idea of edge detector is that edges of an image are associated with location of high gradient in a slightly smooth version of the noisy image.

In the recent decades, to overcome the staircase effects that caused by second-order variational model, fourth-order PDEs have been introduced in image restoration, ^{[1,2,3,4,9,10,12,14,16,20,22,24,25,26,27,28]}. You-Kaveh ^[26] proposed the following functional

where $\triangle$ denote the Laplacian operator. Based on the gradient descent method, this second order functional yields a fourth-order PDE

where $g(s) = k^2/(k^2+s^2)$ and $k$ is an image dependent parameter. This equation does avoid blocky piecewise constant solution. However, it produces speckles in the processed image ^[16]. Many other authors have considered image denoising models based on the minimizer of high-order functionals. Laysaker et al. ^[16] proposed the LLT model

and

they try to minimize the TV of $\nabla u$. Minimizing these two functionals is equivalent to solve the following PDEs respectively

and

where $\mid D^2 u\mid = \sqrt{|u_{xx}|^2+|u_{xy}|^2+|u_{yx}|^2+|u_{yy}|^2}$. These equations have proved to be the improved version of (1.6).

The theoretical analysis showed that fourth-order equations have advantages over second-order equations in some aspects. Fourth-order PDEs usually produce the smooth image of the observed image. This is believed to be a better approximation in smooth region. Therefore, the staircase effect is suggested to be reduced and the recovery image will look better. It is reasonable to conclude that fourth-order diffusion performs better than the second-order models in the aspect of the recovery of smooth regions.

In this paper, to address the problem of denoising images contaminated with additive noise, a fourth-order model is suggested. Using the gradient module of the image to design a speed controlling function. This function indicated where is the edge in the image, thus the new model can preserve edge in this region. The motivation for proposing this model is to overcome certain inconsistencies in second order models founded during the process of recovering smooth regions and better preservation of the fine details. The model is based on solving a nonlinear fourth order degenerate equation with the noisy image as its initial data. By use of Roth's method, we proved the existence and uniqueness of the entropy solution. Additionally, the numerical results demonstrate that the proposed model is superior to PM (Perona, Malik) ^[17] and ROF models.

The rest of this article is organized as follows. In Section 2, we give some preliminaries that we will use. Section 3 is devoted to the proposed model and the proofs of existence and uniqueness of its solution. The difference schemes are presented in Section 4. Numerical experiments are presented in Section 5 and the conclusion of this paper is given in Section 6.

2.   Preliminaries

In this section, we recall some necessary definitions and notations, ^[11,13,14]. We begin with some definitions of the space $BV^2$, which consists of functions $u\in W^{1, 1}(\Omega)$ s.t $\nabla u \in BV(\Omega)$, this space is also denoted by $BH(\Omega)$. To know more about space of bounded Hessian, we refer the reader to ^[4,5,8].

Definition 2.1. Let $\Omega\subseteq \mathbb{R}^n$ be a bounded open domain with Lipschitz boundary. Let $u\in L^1(\Omega)$. Then the $BV^2$ semi-norm of $u$ is characterized by

where $C^2_0(\Omega, X)$ is the space of functions from $\Omega$ to X, 2-times continuously differentiable with compact support and $\phi(x)$is a vector valued function, with $|\phi(x)| = \sqrt{\sum_{i, j = 1}^n(\phi^{ij})^2}$. Here we remark that the space $BV^2$ equipped with $\|u\|_{BV^2(\Omega)} = ||D^2 u||+\|u\|_{L^1(\Omega)}$ is a Banach space.

Definition 2.2. Suppose that $\Omega\subseteq \mathbb{R}^n$ be a bounded open domain with Lipschitz boundary, $u\in L^1(\Omega)$, and $\alpha(x)\geq 0$ is continuous and real function. Then we define the weighted $BV^2$ semi-norm of $u$ as

3.   The proposed model

In this section, we propose a fourth-order image denoising model, with some guidance from previous work ^{[6,14,15,16,26,27]}. There are some benefits of fourth-order models. On the one hand, it can remove high frequency oscillation more effectively than second-order models because the evaluation of the second-order becomes weak in the high frequency area. One the other hand, for the fourth-order model, there is flexibility in employing different functional behaviors in the formulation.

Consider the following boundary value problem

where $\alpha(x) = \frac{1}{\sqrt{1+|G_\sigma *\nabla u_0|^2}},$ $G_\sigma(x)$ is the Gaussian filter with parameter $\sigma$, $u_0(x)$ is the original image, $\Omega$ is bounded domain of $\mathbb{R}^2$ with appropriate smooth boundary, $T > 0$ is fixed, $\overrightarrow{n}$ denote the unit outward normal of the boundary $\partial \Omega$.

The term $\alpha(x)$ is used to enhance edges. In fact, it controls the speed of the diffusion: in the smooth region where $\nabla u_0$ is small, the diffusion is strong. Near possible edges, however, where $\nabla u_0$ is large, the diffusion spread is low. The convolution with $G_\sigma$ should smooth away any large oscillations of noise. Therefore, we can get the smooth image and further preserve the edges in a best way.

Definition 3.1. A measurable function $u:\Omega_T\rightarrow \mathbb{R}$ is an entropy solution of (3.1)–(3.4) in $\Omega_T$ if $u\in C([0, T]; L^2(\Omega))\cap L^\infty(0, T; BV^2(\Omega)), \frac{\partial u}{\partial t}\in L^2(\Omega_T)$ and there exist $z$, such that $\alpha z\in L^\infty(\Omega_T)$ with $\parallel \alpha\, z\parallel_{L^\infty(\Omega_T)}\leq 1, u_t+D_{ij}\alpha z_{ij} = 0$ in $D'(\Omega_T)$ such that

for every $v\in L^\infty (0, T; W_0^{2, 1}(\Omega))$.

Before investigating the existence and uniqueness of problem (3.1)–(3.4), let us consider the following approximate evaluation problem: for $1 < p\leq2$ and $u_{0p}\in W^{2, p}(\Omega)$, we construct the following problem

Lemma 3.1. For any fixed $p$, $1 < p\leq2$, the above problem (3.6)–(3.9) admits a weak solution $u_{p} \in L^{\infty}(0, T; W^{2, p}_{0}(\Omega))\cap C([0, T]; L^2(\Omega))$ and $\frac{\partial u_p}{\partial t}\in L^2(\Omega_T)$ such that

and for any $\varphi \in C_0 ^\infty(\Omega_T)$ the following integral equality holds

with

where C is a constant independent of $p$.

Proof. We apply Rothe's method ^[23], to construct an approximate solution sequence. Divide the interval $[0, T]$ into $n$ equal segments and define $h = \frac{T}{n}$. For any $j$: $1\leq j\leq n$, for any positive integer $n$ and a function $u(x, t)$, denote

For fixed $j$, define the following functional on $W^{2, p}_0(\Omega)$

The idea here is to prove that if $u^{n, j-1}_p$ is known and $u^{n, 0}_p = u_{0p}$, then there is a minimizer for (3.13).

Let $u_m\in W^{2, p}_0(\Omega)\cap L^2(\Omega)$ be a minimizing sequence for $E$. Since $\alpha$ is bounded below, then the sequence $u_m$ is bounded in $W_0^{2, p}(\Omega)$ and $L^2(\Omega)$. Therefore, there exists a subsequence $u_{m_i}$of $u_m$ and a function $u_p^{n, j}\in W^{2, p}_0(\Omega)\cap L^2(\Omega)$ such that as $i\rightarrow \infty$,

From this and the weak lower semicontinuity of the norms, we get

Then $u^{n, j}_p$ is the solution of the Euler equation corresponding to $E(w)$

which implies

for any $\eta(x)\in C^{\infty}_{0}(\Omega)$.

Let $\chi^{n, j}(t)$ be the indicator function of $[h(j-1), hj)$ and

We construct an approximation function as

and

By (3.16), we have

for every $\eta(x)\in C^{\infty}_{0}(\Omega) \text{ a.e. } t \in [0, T]$, which implies that

for every $\varphi\in C^{\infty}_{0}(\Omega_{T})$.

Next, we obtain some estimates for $u^{n}_{p}(x, t)$ and $w^{n}_{p}(x, t)$. Notice that, we choose $\eta(x)\in W^{2, p}_0(\Omega)$ as the test function in (3.16). Let $\eta(x) = u^{n, j}_{p}-u^{n, j-1}_{p}$, we have

Using Young's inequality, we have

That is

For any $m$ with $1\leq m\leq n$, summing (3.18) for $j$ from $1$ to $m$

which implies

where $C$ is a constant independent of $p, n$.

Summing (3.18) for $j$ from 1 to $n$ yield

By the definition of $w^{n}_{p}(x, t)$, we have

Thus

By (3.19), we can obtain

Choosing $\eta(x) = u^{n, j}_{p}$ in (3.16), we get

which implies

By Young's inequality, we have

Thus

This implies

Similar to the proof of (3.22), we also see

Choosing $\eta(x) = u^{n, j-1}_{p}$ in (3.16)

Applying Hölder's inequality and the estimate (3.19), we have

By Young's inequality again yields

Thus

Define $B(u^{n}_{p}) = {\alpha(x)}|\ D_{ij}^{2} u^{n}_{p}|^{p-2}\ D_{ij}^{2} u^{n}_{p}$. Combining (3.19), (3.21), (3.22), (3.25) and (3.26), we conclude that there exist subsequences of $u^{n}_{p}, w^{n}_{p}, \frac{\partial w^{n}_{p}}{\partial t}$ and $B(u^{n}_{p})$, denoted by themselves such that, as $n \to \infty$

holds for some $u_{p}, w_{p}, \zeta$. And we also have

Then let $n\to \infty$ in (3.17),

Next, we show that $u_{p} = w_{p}$. By the definition of $u_{p}$ and $w_{p}$, we have

which combined with (3.20), leads to

Now, it remains to show that $\zeta = B(u_p)$. From (3.17) and the convergence sets (3.28), as $n\to \infty$, we can get

For any $g\in L^p(0, T, W^{2, p}(\Omega))$ and for $j$ from $1$ to $n$, we can obtain, by the monotonicity condition, the inequality

Letting $\eta = u^{n.j}_p$ in (3.16), we obtain

Applying Young's inequality on the first term of (3.34) together with the inequality (3.33) and integrating over $((j-1)h, jh)$, we get

Summing up (3.35) for $j$ from $1$ to $n$, we obtain

Recalling the convergence sets (3.28) and letting $n\to \infty$, (3.36) yields

We can rewrite (3.37) in the form

Letting $\varphi = u_p$ in (3.30), we obtain

Then, substituting (3.39) into (3.38) leads to

Choose $g = u_p-ks$ where $k > 0$ and $D_{ij}^{2}s\in L^\infty(0, T; W^{2, p}(\Omega))$. We then have

Sending $k \to 0$, we obtain

Since $s$ is arbitrary, we see that $\zeta = B(u_p)$.

Now, we prove (3.10), we let $\varphi = u_p(x, t)$ and $\varphi = u_p(x, t_1)$ in (3.30), for $0\leq t_1\leq t\leq t_2\leq T$, we obtain

and

Then

From the above equation, we deduce that

and the proof is completed.

Theorem 3.1. If $u_0\in BV^2(\Omega)$ and $u_0 = 0, \frac{\partial u}{\partial n} = 0, x\in \partial \Omega$ in the sense of trace then the problem (3.1)–(3.4) admits one and only one entropy solution.

Proof. By Lemma 3.1, there exists $u_p$, which is a weak solution of the problem (3.6)–(3.9) and a constant $C$ such that

So, from (3.43), there exists a subsequence of $u_p$, denoted by itself and a function $u\in L^\infty(0, T; BV^2(\Omega))\cap C([0, T]; L^2(\Omega))$ with $\frac{\partial u}{\partial t}\in L^2(\Omega_T)$ such that, as $p\to 1^+$,

and

$\frac{\partial u_p}{\partial t}\rightharpoonup \frac{\partial u}{\partial t}$, weakly in $L^2(\Omega_T)$.

We also have $u_p\to u$ strongly in $L^2(\Omega_T)\quad \text{a.e. } t\in (0, T)$ and

Applying the method in ^[7], we next prove that ${\alpha(x)}| D_{ij}^{2} u_{p}|^{p-2}\ D_{ij}^{2} u_{p}$ is weakly relatively compact in $L^1(\Omega_T)$. Employing (3.43) and Hölder's inequality,

where $C$ is independent of $p$. Thus, $\left\{{\alpha(x)}|D_{ij}^{2} u_{p}|^{p-2}D_{ij}^{2} u_{p}\right\}$ is bounded and equi-integrable in $L^1(\Omega_T)$ and is therefore weakly relatively compact in $L^1(\Omega_T)$. Thus we deduce that as $p\to 1^+$,

$\left\{{\alpha(x)}|D_{ij}^{2} u_{p}|^{p-2}\ D_{ij}^{2} u_{p}\right\} \rightharpoonup \alpha z$, weakly in $L^1(\Omega_T)$.

So we get by Lemma 3.1 and the fact that $\frac{\partial u_p}{\partial t}\rightharpoonup \frac{\partial u}{\partial t}$ in $L^2(\Omega_{T})$,

for every $\varphi (x, t)\in C_0^\infty (\Omega_T)$ and $u_t+D_{ij}\alpha z_{ij} = 0$ in $D'(\Omega_T)$.

Now, it remains to prove that $\parallel \alpha z\parallel_{L^\infty(\Omega_T)}\leq 1$.

For any $k > 0$, setting

$A_{p, k} = \{(x, t)\in \Omega_T:\mid D_{ij}^2u_p\mid > k\}$, we have that

As above, there exists a function $g_k\in L^1(\Omega_T)$ such that

$\alpha(x)|D_{ij}^{2} u_{p}|^{p-2}\ D_{ij}^{2} u_{p}\chi_{A_{p, k}}\rightharpoonup g_{k}$, as $p\to 1^+$ weakly in $L^1(\Omega_{T})$,

where $\chi_{A_{p, k}}$ is the indicator function of ${A_{p, k}}$. Now for any $\phi \in L^\infty (\Omega_T)$ with

$\|\phi\|_{L^\infty(\Omega_T)}\leq 1$, by the definition of $\chi_{A_{p, k}}$, we see that

Letting $p\to 1^+$, we have

Since we have that

letting $p\to 1^+$, we obtain that $\alpha(x)\mid D_{ij}^{2} u_{p}\mid^{p-2}\ D_{ij}^{2} u_{p}\chi_{\Omega_T/A_{p, k}}$ weakly converges in $L^1(\Omega_T)$ to some function $f_k\in L^1(\Omega_T)$ with $\|f_k\|_{L^\infty(\Omega_T)}\leq 1$. Since, for any $k > 0$, we may write $\alpha z = f_k+g_k$ with $\|f_k\|_{L^\infty(\Omega_T)}\leq 1$ and $g_k$ satisfies (3.45), it is easily follows that $\|\alpha z\|_{L^\infty(\Omega_T)}\leq 1$.

Next, we verify the solution definition inequality (3.5). For any $v_n\in C_0^\infty(\Omega_T)$ and taking $\varphi = (u_p-v_n)\xi(t)$ in (3.11), we have

Letting $p\to 1^+$,

Then for any $v\in L^\infty (0, T; W_0^{2, 1}(\Omega))$, letting $n\to \infty$,

Since $\xi(t)$ is arbitrary, we have

for every $v\in L^\infty(0, T; W^{2, 1}_0(\Omega))$ and a.e. on $[0, T]$.

Finally, we prove the uniqueness of the entropy solution. Let $u_1, u_2$ both be entropy solution with data $u_{10}, u_{20}$. Then there exists $\alpha z_1, \alpha z_2\in L^\infty(\Omega_T)$ such that

and

for every $v\in L^\infty(0, T; W^{2, 1}_0(\Omega))$ and a.e. on $[0, T]$. Let $u_{1n}, u_{2n}\in L^\infty(0, T; W^{2, p}_0(\Omega))$ be approximates functions, respectively, for $u_1$ and $u_2$, such that

and

a.e. on$[0, T]$. Taking $v = u_{2n}$ in (3.46) and $v = u_{1n}$ in (3.47), adding the two equations and rearranging the result, we obtain

So integrating from $0$ to $t$ and letting $n\to\infty$, we get

The proof is completed.

4.   Difference schemes

In this section, assuming $\tau$ to be the time step size and $h$ the space grid size, we discretize time and space as follows:

where $Ih\times Jh$ is the size of the original image. Let $u_{i, j}^n$ denote approximations of $u(n\tau, ih, jh)$. We define the discrete approximation:

The discrete explicit scheme of the problem can be written as

Here the MATLAB function "conv2" is used to represent the two-dimensional discrete Convolution Transform of the matrix $u_{i, j}$.

5.   Numerical experiments

In this section, we demonstrate the performance of our model in denoising images involving Gaussian white noise. We applied difference equations discussed in section 4 and compared the results with the results of ROF model ^[18] and PM model ^[17]. We used step size $\tau = 0.02$, gride size $h = 1$ and $\lambda = 0$.

At the end of the denoising process, the peak signal to noise ratio (PSNR), mean absolute-deviation error (MAE) and structure similarity index measure (SSIM) values were recorded to measure the denoising performance. The values are given by the following formulas:

and

where $\left| \text{max}\, u_0- \text{min}\, u_0 \right|$ gives the gray-scale range of the original image, $u_0$ and $u$ denote, respectively, the original image and the denoised image, $I\times J$ is the dimension of image.

SSIM, designed by Wang et al.^[21], is a quality used to measure the similarity between any two images. Given any two images $u$ and $u_0$, SSIM is given by the formula

$\text{L}(u, u_0) = \frac{2\mu_u\mu_{u_0}+k_1}{\mu^2_u+\mu^2_{u_0}+k_1},$ compares the two images' mean luminance $\mu_u$ and $\mu_{u_0}.$ The maximal value of $\text{L}(u, u_0) = 1$, if $\mu_u = \mu_{u_0}$, $\text{C}(u, u_0) = \frac{2\sigma_u\sigma_{u_0}+k_2}{\sigma_u^2+\sigma_{u_0}^2+k_2},$ measures the closeness of contrast of the two images $u$ and $u_0.$ Contrast is determined in terms of standard deviation, $\sigma.$ Contrast comparison measure $\text{C}(u, u_0) = 1$ maximally if and only if $\sigma_u = \sigma_{u_0};$ that is, when the images have equal contrast.

$\text{R}(u, u_0) = \frac{\sigma_{uu_0}+k_3}{\sigma_u\sigma_{u_0}+k_3},$ is a structure comparison measure which determines the correlation between the images $u$ and $u_0$, where $\sigma_{uu_0}$ is covariance between $u$ and $u_0$. It attains maximal value of $1$ if, structurally, the two images coincide, but its value is equal to zero when there is absolutely no structural coincidence. The quantities $k_1, k_2$ and $k_3$ are small positive perturbations that avert the possibility of having zero denominators.

Two test images of "Cameraman" and "Peppers" are corrupted by white Gaussian noise with standard deviation (SD) of 30, (Figures 1 and 2). Tables 1 and 2, present the numerical results of restoration of Cameraman image, (Figure 1), and those of the Peppers image, (Figure 2). The comparisons are based on PSNR, MAE and SSIM. The proposed method shows the best performance with respect to PSNR, MAE and SSIM.

Our first example is Cameraman image, which is displayed in Figure 1a and 1b is its degraded version. Furthermore, Figure 1c, 1d and 1e, are portions of the recovered images with the proposed model, ROF model and PM model, respectively. It is clear that our method can overcome the staircase effect that caused by the second order method.

The second example is Peppers image, which is displayed in Figure 2a, its degraded version is showed in Figure 2b. Basically, Figure 2c, 2d and 2e, are portions of the recovered images by the proposed model, ROF model and PM model, respectively. It is evident that, our method yields good results in restoring image since it avoids the staircase effect that caused by the second order method while, at the same time, handle edges in a best way.

Similarly, the two test images are corrupted by white Gaussian noise with SD of 20, (Figures 3 and 4). Tables 3 and 4, present the numerical results of restoration of Cameraman image, (Figure 3), and those of the Peppers image, (Figure 4). The comparisons are based on PSNR, MAE and SSIM. Here again, the proposed method shows the best performance with respect to PSNR, MAE and SSIM. In Figure 3a and 3b we display Cameraman image and the noisy version. Figure 3c, 3d and 3e, are portions of the recovered images with the proposed model, ROF model and PM model, respectively. We display Peppers image and the degraded version in Figure 4a and 4b. Figure 4c, 4d and 4e, are portions of the recovered images with the proposed model, ROF model and PM model, respectively. Here also, the proposed model yields better results in denoising image while handling edges in a best way.

Not surprisingly, although the edges are preserved, the staircase effect is visible for the second order models, and there are some speckles in the processed images, with an example given in Figure 2. Comparing the images processed by our model and the original images, we can observed that, the differences are insignificant. The edges are preserved and no speckles appear in the processed images.

Finally, to illustrate the superiority of the proposed model over other related fourth-order models, we compared our results with YK model ^[26] and LLT model ^[16]. Barbara and Lena images have been corrupted by white Gaussian noise with SD of 30 (Figures 5 and 6) and SD of 20 (Figures 7 and 8). Numerical results for the images are tabulated in Table 4. Besides getting better outcomes, as evident from the results (see Figures 5 and 6), the model tackles the speckles caused by YK model at the same time.

In Figures 6 and 8, the results of Lena Image have been displayed. In Figure 6, the test image Lena and its noisy version degraded by Gaussian noise with SD of 30 are shown in the sections (a) and (b), sections (c) to (e) are the results of the YK model, LLT model and the proposed one. Similarly, in Figure 8, the test image Lena and its noisy version degraded by Gaussian noise with SD of 20 are shown with the same order described above. The last image in section (e) of Figures 6 and 8 are the results of our suggested filter in which the extent of the denoising performance is noticeably better than competitor filter.

6.   Conclusions

In this article, we proposed a fourth-order image denoising model. The model was based on solving a fourth order partial differential equation by defining its corresponding functional. We proved, by use of Rothe's method, the existence and uniqueness of the entropy solution of the equation. Compared with the well known ROF and PM models, numerical results showed that our model perform better image recovery and can overcome staircase effects.

Conflict of interest

The authors declare that they have no conflict of interests whatsoever and do approve the publication of this paper.