An improved coupled PDE system applied to the inverse image denoising problem

1.   Introduction

Currently, Image restoration is one of the outstanding challenging problems in both image processing and computer vision with numberless applications. Its goal is to remove noise from a degraded image to restore the original one. Generally, given a noisy image function defined on $\Omega$, with $\Omega \subset \mathbb{R}^2$ an open and bounded domain, the image denoising process can be modelled as $f = u+n$ where $f$ represents the observed noisy image, $u$ is the true image, and $n$ is the noise component. To solve this inverse problem, there are several approaches such as PDE based technique ^[3,4], wavelet based technique ^[5], stochastic approach ^[6], patch based technique ^[7].

PDE based denoising models can be classified in two categories, namely variational and non variational PDE models. In variational PDE models, the solution is obtained as the steady-state solution of an evolution equation corresponding to the Euler–Lagrange equation of the energy functional. In non-variational PDE models, the PDE is proposed directly on the diffusion equations (systems), without thinking of any energy. Whatever the used approach, one should consider a certain diffusion equation (system).

In recent years, partial differential equations (PDEs) have become one of the most effective tools for image processing. Among them, the ROF model, which was first introduced by Rudin et al. ^[8], is one of the most popular and given by

where $\lambda > 0$ is parameter regularization, $\Vert u-f\Vert_{L^2(\Omega)}^2$ is a fidelity term, and $\int_{\Omega} \vert \nabla u \vert$ is a regularization term. The corresponding PDE which the evolution of the Euler-Lagrange equation for $E(u)$ is given by following equation

This last model is named TV model, it has been the subject of several studies (see ^[9,10]) and has successfully contributed to preserve edges during the restoration process. Given the success of TV-based diffusion, various modifications have been investigated (see ^[11,12]). Though the TV model performs very well for image denoising and edge protecting, it may also destroy small details, such as textures (see ^[13]). To overcome this drawback, Meyer ^[14] proposed a new minimization method by introducing a weaker norm which is more appropriate to represent the oscillatory patterns and small details in the textured image. Subsequently, Osher et al. ^[13] proposed the following $4$th order denoising model which used the $H^{-1}$ norm with the TV minimization. The minimization energy has the form

The corresponding PDE associated to Euler-Lagrange is given by the following equation :

This model is called the $TV-H^{-1}$ denoising model. This model is a highly nonlinear PDE and its numerical solution is not obvious. To solve this problem, Guo et al. ^[1] suggested a reaction-diffusion system applied to image restoration and image decomposition into cartoon and texture. This idea is introduced by ^[12] where the initial image $f$ is decomposed into carton part $u$ and texture or noise $(f-u)$. The Euler-Lagrange of the functional energy (1.3) becomes

To overcome the difficulty of 4th order degrees in (1.5), the authors transform it into the following two coupled second order equations:

and the evolutionary form of the steady system (1.6) is given by the following PDE system:

Recently, Atlas et al. ^[15] proposed a new model based on a nonlinear fractional reaction-diffusion system. This model is based on the following minimization procedure:

where the norm $H^{-s}, \; s > 0$ is introduced by Giga ^[16]. The corresponding PDE associated of its Euler-Lagrange will be

where $H^{-s}$ is negative Hilbert-Sobolev space with $s \in (0, 1]$ For the nonlinear regularization term, the authors propose an Orlicz-operator and generalize Eq (1.9), which becomes:

where the function $A_{\mu(x)} : \Omega \times \mathbb{R}\rightarrow \mathbb{R}$ is chosen to be monotonically close to $1$ when $t$ tends to infinity (for more details see ^[15]). To solve the problem (1.10), the authors follow the decomposition suggested in ^[1,2] introducing the following coupled system

furthermore, the authors consider the solution $(u, v)$ of system (1.11) as a steady state of an evolutionary fractional reaction-diffusion system investigated in ^[15].

More recently, Halim et al. ^[17] proposed a non-variational denoising model by modifying the $TV-TV2$ model introduced in ^[18]. They proposed the following model:

This model called $TV-L^2-H^{-1}$, it contains the diffusion terms of $TV-L^2$ and $TV-H^{-1}$ model. However, it still suffers from the staircasing effect and the combination between the second and fourth-order diffusion terms leads to over-smooth homogeneous areas.

To overcome these limitations, we propose a nonlinear high order system, using the Weickert model ^[19], which is more efficient in reducing high noise intensity, combined with a decomposition suggested by Guo et al. ^[1] to control the diffusivity in smooth areas. This model corrects the disadvantages (staircasing phenomenon and over smoothing) of the diffusion term proposed in Eq (1.12). We can see that in the Figure 1, where a denoising example is carried-out using the two diffusion operators. However, since the Weickert operator is strongly nonlinear, the computation of the solution using classical numerical approximation generates some blur effect. Hence, the introduction of the decomposition procedure, which will facilitate the treatment of the nonlinear part of the Weickert operator. However this decomposition is in general ill-posed and existence of the looked for solution is not ensured. In this work, we will check the existence and uniqueness of the solution in a suitable functional space inspired by the well-known Schauder fixed point theorem.

The organization of the main paper is given as follows. In section 2, we describe the new proposed model. Section 3 is devoted to the existence and uniqueness results of the solution to the proposed equation using Schauder fixed point theorem. After, we give a brief discretization part of the proposed coupled model PDE. At last, section 4 is devoted to numerical results and comparative experiments to improve our model.

2.   Description of the proposed model

In this section, we have proposed a denoising model of non-variational type by modifying the diffusion term in Eq (1.12) to overcome its drawbacks. To take advantage of both the Perona-Malik model ^[20] and the Weickert model ^[19], we propose the new model given by the following equation:

where $D$ is an anisotropic diffusion tensor and $J_{\rho}$ is the structure tensor defined by

where $*$ is the convolution product and $\otimes$ is the outer product defined such as $\nabla X_{\sigma}\otimes \nabla X_{\sigma} = \nabla X_{\sigma} (\nabla X_{\sigma})^t$. Here $X_{\sigma}$ is constructed using a convolution of $X$ with a Gaussian kernel, while $K_{\rho}$ and $K_{\sigma}$ represent two Gaussian convolution kernels such as $K_{\tau}(x) = \frac{1}{2\pi \tau^{2}}\exp(-\frac{\vert x\vert^{2}}{2 \tau^{2}})$. The function $D$ is calculated using the tensor $J_{\rho}$ eigenvalues and the eigenvectors as follows

where $\upsilon_{+/-}$ and $\theta_{+/-}$ are respectively the eigenvalues and the eigenvectors of the tensor structure $J_{\rho}$, the eigenvalues $\upsilon_{+/-}$ are calculated as

While the functions $\varphi_{+}$ and $\varphi_{-}$ represent the isotropic or anisotropic behavior of the smoothing on the image regions. Recently, an efficient choice of these functions is introduced in ^[21], where the authors consider the behavior of the Weickert model according to the two directions $\theta_{+}$ and $\theta_{-}$. The proposed coefficients take into account the diffusion in the vicinity of the contours and corners where the eigenvalues $\upsilon_{+}$ and $\upsilon_{-}$ are very high. This choice is proposed as follows:

where $k_{1}$ and $k_{2}$ are two thresholds defining the diffusion with respect to the directions $\theta_{+}$ and $\theta_{-}$, respectively.

Equation (2.1) is a highly nonlinear PDE and therefore its numerical solution is a non-trivial task. To overcome this difficulty, we follow the decomposition suggested in ^[1], introducing a splitting into the coupled system:

Moreover, the solution $(u, v)$ of system (2.6) is a stable state of the following evolutionary reaction-diffusion system

In order to establish the existence of the solution $u$ associated to the problem (2.7), we introduce the following set of hypotheses $(H)$ :

$H_1$- The tensors $D$ in $C^\infty(\mathbb{R}^{2\times 2}, \mathbb{R}^{2\times 2})$, positive-definite matrix and coercive with the coercivity constant is $\beta$.

$H_2$- The initial condition $f$ in $L^{2}(\Omega)$.

$H_3$- The constants $\beta$, $\alpha$, $\tau$, $\rho$, $\sigma$, $k_1$ and $k_2$ are positives.

3.   Existence and uniqueness

In this section, we give the variational formulation of the problem (2.7) and a priori estimations of the solution $(u, v)$. Indeed, the variational formulation of the problem (2.7) is stated as follows

Inspired by the techniques employed by Catté et al. ^[22], for a nonlinear equation, we will show the existence result for our coupled problem of two equations. The following lemma gives some a priori estimations of the solution.

Lemma 1. Assume that assumptions $(H)$ are satisfied, then there exists a positive constant $C$, such that the weak solution of the problem (3.1) satisfies the following estimations

Proof. We take $\phi = u$ in the first equation of (3.1) and $\phi = v$ in the second, we get

We integrate for $t\in(0, \tau]$ with $\tau\in (0, T]$, then we have

Multiplying the second equation of (3.3) by $\alpha$ and using the initial data $v(0) = 0$ and $u(0) = f$, then the equation becomes

we add the two equations of (3.4), we get

On the one hand, applying Young's inequality we have

On the other hand, thanks to the coercivity of $D$, we see that

Using (3.7) and (3.6), Eq (3.5) becomes

Setting $K(\tau) = \int_{\Omega}v(\tau)^{2}$, we see that

Using the Gronwall's inequality, we get

which implies that

where the constant depends on $T$, we deduce that $v$ is bounded in $L^{\infty}(0, T;L^{2}(\Omega))$. Hence

By taking Eq (3.8) and by adding the term $\beta \int_{0}^{\tau}\int_{\Omega} u^{2}+\alpha \int_{0}^{\tau}\int_{\Omega} v^{2}$ to both sides we get

which means that:

Let's prove now the estimation of $\dfrac{\partial u}{\partial t}$ and $\dfrac{\partial v}{\partial t}$. For that, let's back to Eq (3.1), we have for all $\phi$

using hypotheses $H_1$ and Hölder inequality we get

or $\mid\mid\nabla\phi\mid\mid_{L^{2}(\Omega)}\leq \mid\mid\phi\mid\mid_{H^{1}(\Omega)}$ and $\mid\mid\phi\mid\mid_{L^{2}(\Omega)}\leq \mid\mid\phi\mid\mid_{H^{1}(\Omega)}$, then the previous equation becomes

this means

integrating over $t\in(0, T]$, we find

and according to estimates in $(11)$, we conclude that

which ends the proof.

The following theorem shows the existence and uniqueness of a solution to the proposed model

Theorem 1. Let $f\in L^{2}(\Omega)$ and $T > 0$, under the assumptions above, the problem (2.7) admits a unique weak solution $(u, v)\in C(0, T;L^{2}(\Omega))\cap L^{2}(0, T;H^{1}(\Omega))$.

3.1.   Existence

For the proof of existence, we use the classical Schauder fixed point theorem [23,theorem 2.A,p. 56]. To define the fixed point operator, we introduce first these spaces

which is a Hilbert space equipped with the norm

and

which is a nonempty, convex and weakly compact subset of $V(0, T)$. The estimations introduced in the functional space $V_0$ are deduced from the precedent lemma 1.

The Schauder's fixed point operator is given as follows

where $u_{w}$ is the solution associated to $(w, v)$ for the following problem

which is now linear in $u$. As a result and using the theoretical existence of parabolic equations results [24,Theorem 3,p. 356], we ensure that the problem (3.17) has a unique solution $(u_{w}, v)$ in $V(0, T)$. The existence of a weak solution for the problem (3.1) is equivalent to the existence of a fixed point for the operator $L$. For this reason, we apply Schauder's fixed point theorem, which requires only to prove that the mapping $L$ is weakly continuous. For that, let $(w_{n})_{n\in\mathbb{N}}$ be a sequence in $V_0$ such that $w_{n}\rightharpoonup w$ in $V_0$, and $u_{n} = L(w_{n})$, we should prove that

Using the compact inclusions of Sobolev spaces ^[25] and by estimations established in lemma 1, there exists a subsequence noted also $(w_{n}, v_{n})_{n\in\mathbb{N}}$ and $(u_{n}, v_{n})$ such that

Using these convergences and by the uniqueness of the solution of (3.17), we have

This proves that $L$ is weakly continuous. Finally from the Schauder fixed point theorem, the operator $L$ admits a fixed point solution to the problem (3.1). Moreover, since $u, v\in L^{2}(0, T;H^{1}(\Omega))$ and $\frac{\partial u}{\partial t}, \frac{\partial v}{\partial t}\in L^{2}(0, T;H^{1}(\Omega)^{'})$, by Aubin's theorem ^[26], we deduce that $u, v\in C(0, T;L^{2}(\Omega))$.

3.2.   Uniqueness

In order to show that the solution of (2.7) is unique, we consider two different solutions $(u_{1}, v_{1})$ and $(u_{2}, v_{2})$ to the problem (2.7). Then we have

and

By subtracting both variational formulation of $(u_{1}, v_{1})$ and $(u_{2}, v_{2})$, we obtain

We take $\phi = u_{1}-u_{2}$ in the first equation of (3.21) and $\phi = v_{1}-v_{2}$ in the second, we get

Now, multiplying the second equation of (3.22) by $\alpha$ and adding the two equations we get

Thanks to the coercivity of $D$, we have

Since the operator $D(J_{\rho})$ is smooth enough, we have

Using (3.23) and hypothesis $(H_{1})$, and Holder inequality, we have

Integrating over $[0, s)$ with $s\in(0, T]$ and applying Young's inequality for enough epsilon, we have

For $\varepsilon < 2\beta$, we have

Using Gronwall's inequality we have

Now back to the second equation of (3.22), we have

which means that $v_{1} = v_{2}$,

Hence

4.   Discrete setting

In this section, we are interested in the numerical approximation of the proposed model. In fact, to compute numerically the problem (2.7), we carry out a fully discretization of the model using finite difference method.

Assume $k$ to be the time step size and $h$ the spatial grid size, we discretize time and space as follows:

Denote $u^{n}_{i, j}$, $v^{n}_{i, j}$ and $f_{i, j}$ the approximations of $u(t_n, x_i, y_j)$, $v(t_n, x_i, y_j)$ and $f(x, y)$ respectively, then we define the discrete approximations :

it remains now to discretize the term of diffusion $div(D(J_{\rho}(\nabla u_{\sigma}))\nabla u)$, to do so, we consider the method represented by Weickert [19,Chapter 3], which consists of acting by convolution on the discrete image in any pixel with kernel, called a non-negative stencil. For that we write the diffusion tensor such as $\begin{pmatrix} a&b\\ b&c \end{pmatrix}$ which is calculated by the structure tensor $J_{\rho}(\nabla u_{\sigma})$, then the term of diffusion can be rewritten as

Hence, the discretization of the diffusion term is given by

then this standard discretization is equivalent to acting by convolution on the image $u$ at any point $(i, j)$ by a matrix of size $3\times3$ called by the Stencil matrix denoted by $S$ [19,Chapter 3,p. 95], if $A(u)$ indicates the matrix containing the elements of term $div(D(J_{\rho}(\nabla u_{\sigma}))\nabla u)$ after the discretization, then an element of $A(u)$ is calculated by

Finally, the discrete explicit scheme of the problem (2.7) could be written as

for all $k\geq 0$, $0\leq i\leq N$ and $0\leq j\leq M$.

5.   Numerical results

In this section, we illustrate the performance of the proposed model for image denoising. For a fair comparison, these methods are implemented using Matlab 2012a on the platform: 3 GHz dual-core CPU and 8 Gbytes RAM. We use the relative error ($E_r$) as a stopping criteria of the following numerical simulations, which is defined by

where $u^{n}$ is the restored image at the iteration $n$.

We have tested the proposed multi-frame SR method on a large database image using different applications, we present in the following experiments only some of them. Two commonly used quality measures are considered to evaluate the denoising performance, i.e., PSNR and the SSIM indexes. For the following tests we fix the parameters of our method such as $k = 0.15$, $\alpha = 0.01$, $\rho = 1.4$ $\sigma = 1.6$ and $maxIter = 500$. While the values of $k_1$ and $k_2$ are chosen with respect to the best PSNR values. For the following tests, we chose $k_1\in [30,150]$ and $k_2\in [10,250]$, respectively.

In order to validate our model in computing the component $u$ and $v$, we begin with validation tests, where we select four noisy images (with different standard deviation $\sigma$). In Figure 2, we have presented the restored image $u$ and the computed component $v$ for each test. From this figure it is clear that our model is able to recover image features even if the noise level is high.

For a convenient comparison, we compare our method with some competitive denoising methods, among them: the TV-H$^{-1}$ ^[13], the telegraph coupled partial differential equation (TCPDE) proposed in ^[27], TV-L$^{2}$-H$^{-1}$ ^[17], the fourth-order PDE (FOPDE) ^[28] and also the nonlinear fractional reaction-diffusion system (NFAD) ^[15]. Note that concerning the chosen parameters for the proposed algorithm, we take $k_1 = 55$, $k_2 = 45$, $k = 0.1$, $\alpha = 0.01$, $\rho = 3.2$, time step size $k = 0.1$, $\sigma = 2.2$ and $maxIter = 1000$. In contrast, the used parameters for the compared methods are chosen inspired from the given parameters in numerical results of their respective papers. In fact, the used parameters for these methods (including the proposed one) corresponding to the following tests are depicted in Table 1.

● Firstly, we consider the " Motif" image and we configure a simulated test. We add a Gaussian noise with parameter $\sigma = 30$. The restored images using different denoising approaches are shown in Figure 3 and the associated 3D surfaces in Figure 4. We can see that the proposed reaction-diffusion equation outperforms the other methods and we can see the difference better in the smooth areas, where the staircasing effect is efficiently avoided.

● The next experiments deals with Drop image, which is contaminated by a Gaussian noise with standard deviation $\sigma = 40$. The aim is to restore the clean image and compare it with the other denoising methods. The obtained results for the three tests are shown in Figure 5. As the previous test, the proposed PDE-constrained equation gives the better result. We can see the difference better in the 3D plots, where we show the image region surfaces of each restored image in Figure 6. As a result, it is quite visible that the proposed method demonstrates clearly it performance as seen in the restored image: edges are efficiently preserved, homogeneous areas are free from noise and low-contrast objects are quite visible. This conforms the fact, that the idea of computing the couple $u$ and $v$ with the tensor diffusion deals with the competitive image denoising techniques.

● For the third test we use Brain image, we added a Gaussian noise with $\sigma = 50$. In Figure 7, we present the restored image using different enhancement methods. Once again, we can see clearly the robustness of the proposed PDE compared with the other ones. This also confirmed in Figure 8, where we present the respective 3D surfaces. A we can see the image regions recovered by our method look very close to the ones of the original image, which is not the case of the compared methods.

To give a fair and quantitative comparison with the other methods, we use the PSNR and SSIM measures. Effectively, we show the evolution of the PSNR and SSIM with respect to the first 1000 iterations for the proposed and compared methods in Figure 9. We can see once again that always our approach reaches the best PSNR and SSIM values for the three tests.

Finally, we propose two real MRI tests downloaded from^*, where the images are corrupted by unknown noise supposed to be of Poisson type. The aim is to recover a clean image with strong edges and considerable contrast. Figure 8 illustrated the recovered MRI Head using different denoising approaches, while Figure 11 presents the obtained clean image $u$ for the real MRI Knee image with comparison to other methods. The final two results obtained using the proposed equation clearly indicates that the significant edges and textures are recovered with less blurring and staircasing. As a conclusion, from these real images, it is easy to observe that the proposed PDE is more efficient to remove almost all noise particles with effective structure preservation without information about the noise type.

^*www.mr-tip.com

6.   Conclusions

In this paper, we have proposed a coupled PDE denoising model. This PDE is based on a diffusive tensor with efficient diffusion properties along contours. The well-posedness study of the proposed PDE is also achieved using the Schauder theorem. Finally, experiment results demonstrate the ability of the proposed model compared with the results of some competitive models. Since the two parameters$k_1$ and $k_2$ affect the smoothing phenomena near strong edges, an optimization procedure will be of interest to compute these parameters, which will help the model to control the diffusion tackled by the tensor term.

Acknowledgement

The authors are grateful to the two reviewers for their valuable comments, which improved the readability of this paper.

Conflict of interest

The authors declare that they have no conflict of interest.