Some class of parabolic systems applied to image processing

In this paper, we are interested in the mathematical and numerical study of a variational model derived as Reaction-Diffusion System for image denoising. We use a nonlinear regularization of total variation (TV) operator's, combined with a decomposition approach of $H^{-1}$ norm suggested by Guo and al. ([19],[20]). Based on Galerkin's method, we prove the existence and uniqueness of the solution on Orlicz space for the proposed model. At last, compared experimental results distinctly demonstrate the superiority of our model, in term of removing noise while preserving the edges and reducing staircase effect.


1.
Introduction. In recent years, many mathematicians have been attracted by image processing and computer vision. Several field of application of the image processing, we quote that the fundamental problem in image processing is the image restoration. Their methods has grown with the massive production of digital images and movies, often taken in poor conditions, which gives a noisy images. Image denoising refers to the process of the recovering an image contaminated by the noise. The common starting point is an image which is a collection of information about N pixels. Generally, noised images can be modeled as where f (i) represents the i th pixel of the observed noisy image, u(i) is the true image and n(i) the noise which often considered the stationary Gaussian with zero mean and variance σ 2 .
The challenge of image denoising is to remove the noise while maintaining and recovering the features and details of image as much as possible.
Over decades, the nonlinear diffusion and partial differential equation based variational models have become popular and useful tools for image processing (see [3], [5], [9], [29], [30] and references mentioned therein). Generally, the restored (denoised) image u : Ω → R is computed from the following minimization problem: where Ω ⊂ R 2 is an open and bounded domain, X is a Banach space, Ψ : R 2 ×R → R is a given function, λ > 0 is a weight parameter, ||f − u|| 2 X is the fidelity term and Ω Ψ(x, |∇u|)dx is a regularizing term to remove the noise. Notice that a different choice of Ψ and X correspond to different models. In the case where X = L 2 (Ω), the problem (1) is formally equivalent to the following nonlinear equation where for a.e. x ∈ Ω, Ψ(x, s) = s 0 Φ(x, t)dt.
Perona and Malik (case Φ(x, s) = s 1 + s 2 a.e. in Ω ) have presented the first nonlinear partial differential equations model for denoising in a class of models commonly referred to as anisotropic diffusion. By choosing Φ(x, s) = 1 a.e. in Ω, Rudin et al. ( ROF) in [30] have proposed an alternative model called total variation model, that can protect the details of image better in the course of denoising. When the parameter λ is too small, the smaller details of image are destroyed. That is why, Meyer [27] proposed a new minimization problem, changing in (1) the L 2 norm of (f − u) by a weaker norm, more appropriate to present textured or oscillatory patterns, but, the drawback of this model, we cannot express directly the associated Euler-Lagrange equation with respect to u. To overcome this problem, Osher et al. in [28], has proposed an alternative practical approach that changing this weaker norm by H −1 . Since the H −1 norm is defined as ||.|| 2 the minimization problem (1) in this case is formally associated to Euler-Lagrange equation We recall that the Total variation model, is not well defined at locations where |∇u| = 0 due to the presence of the term 1 |∇u| and it is known that the solution of this problem is defined in the non standard variational sense (see [4] and [13]). Although the numerical simulations were done by regularized |∇u| with ε + |∇u|.
In current work, we propose a reaction diffusion system, by using a nonlinear regularization of total variation model combined with a decomposition suggested by Guo et al. ([19], [20]), which is written as : where Φ γ is defined by The parameters ε, λ and γ are nonnegative and A typical example of such B are B(x, t) = t 2 exp(t) or t α(x) log γ (1+t) for a.e. x ∈ Ω, where α is a given continuous function. When γ is close to 0, we have and when [|∇u| ≤ ε], Φ γ (x, |∇u|) ∇u |∇u| behaves like exp(|∇u|)∇u or |∇u| α(x)−2 ∇u .
That is, at edges where |∇u| is large, the proposed operators can be seen as a regularization of TV operator. Our model is well posed in a weak sense and a numerical simulations substantially reduces the staircase effect.
For the mathematical analysis, it is reasonable to study the solutions of our problem in the Orlicz space. For this, we give some basic definitions and preliminaries needed to state the results in the next section. In Section 3, we show the existence and the uniqueness result of solutions for the proposed model (4). Finally, section 4 is devoted to numerical results and comments to improve our model. where a : R + → R + is given by a(t) = sup{s : a(s) ≤ t} (see [2], [23] and [24]).
The N-function M is said to satisfy the ∆ 2 condition if, for some k > 0: when this inequality holds only for t ≥ t 0 > 0, M is said to satisfy the ∆ 2 condition near infinity.
Let Ω be an open subset of R N . The Orlicz class L M (Ω) (resp. the Orlicz space L M (Ω)) is defined as the set of (equivalence classes of) real-valued measurable functions u on Ω such that:  x the distributional derivative on Q of order α with respect to the variable x ∈ R N . The inhomogeneous Orlicz-Sobolev spaces are defined as follows The last space is a subspace of the first one, and both are Banach spaces under the norm We can easily prove that they form a complementary system when Ω satisfies the segment property. These spaces are considered as subspaces of the product space ΠL M (Q) which have as many copies as there is α-order derivatives, |α| ≤ m. We shall also consider the weak topologies σ(ΠL M , ΠE M ) and (Ω)-valued and is strongly measurable. Furthermore the following imbedding holds: . We can easily show as in [17] that when Ω has the segment property then each element u of the closure of D(Q) with respect of the weak * topology , this space will be denoted by W m,x Thus both sides of the last inequality are equivalent norms on W m,x 0 L M (Q). We have then the following complementary system . It is also, except for an isomorphism, the quotient of ΠL M by the polar set W m,x 0 E M (Q) ⊥ , and will be denoted by F = W −m,x L M (Q) and it is shown that This space will be equipped with the usual quotient norm where the infimum is taken on all possible decomposition The space F 0 is then given by and is denoted by To proof our existence Theorem, we need the following corollary introduced in [1].
3. Existence. For simplicity, we use the function a γ : Thanks to (5) and (6), we see that a γ is a Caratheodory function (i.e is measurable in x ∈ Ω for all ξ ∈ R N and continuous in ξ ∈ R N for a.e. x ∈ Ω) with a γ (., 0) = 0. Now, we prove the following main result: Then for every λ > 0, there exists a unique solution for every τ ∈ (0, T ] and every test-functions Where < ., . > resp << ., . >> denotes the duality bracket between L 2 0, T ; (H 1 (Ω)) and L 2 0, T ; (H 1 (Ω)) resp. between W 1,x L A (Q) ∩ L 2 (Q) and W −1,x Proof of Theorem 3.1. To prove the existence of a weak solution, we shall use the nonlinear Galerkin method. The idea is to solve the problem in a finite-dimensional space firstly and we are looking estimates that allows us to pass to the limit. We decompose the proof of Theorem 3.1 into four parts: first, we write the approximate solution, then we give a priori estimates, after we pass to the limit and finally we close the demonstration by proving the uniqueness of the solution.

First
Step. Approximate solution. we choose a sequence with s large enough such that H s (Ω) is continuously embedded in C 1 (Ω). We consider the following sequence for approximating solutions of the problem (4): where c m k , d m k : [0, T ) → R are supposed to be measurable bounded functions. For the initial conditions, we choose the coefficients as where u m (. The coefficients c m k , d m k are obtained from the following system of ordinary differential equations where for i = 1, 2, M i (t) ∈ L 1 (0, τ ) and C i (m, |Ω|) is a constant which depends on m, f and |Ω|. Then, thanks to the existence result of ordinary differential equations (cf. [22]), the system for k ∈ {1, 2, ..., m}, has a continuous solution c m on an interval (0, τ ), τ > 0 and may depend on m. Using a standard arguments, It is not difficult to show, that the local solution constructed above can be extended to the whole interval [0, T ) independent of m. To passing to the limit in (12)- (13) and proving the existence of u(x, t) and w(x, t), we need the following a priori estimates lemma.

Second
Step. A priori estimates Lemma 3.2. Let (u m , w m ) be a solution of the problem (12)- (13). Then, we have where the constant C is independent of m.
Proof. Let τ < T , taking u m χ (0,τ ) (respectively w m χ (0,τ ) ) as a test function in (12) (respectively in (13) ), we obtain Using the initial condition u m (x, 0) = f (x) and w(x, 0) = 0, then the equation becomes On the one hand, applying Young's inequality, we have On the other hand, thanks to the definition of a γ , we see that Using (20) and (21), the equation (19) becomes Now, setting Θ n (τ ) = Ω w m (τ ) 2 dx, we observe that Using Gronwall's inequality, we get where the constant C depends only on T and λ. Consequently, we have Then, we deduce that w m is bounded in L ∞ 0, T, L 2 (Ω) ∩L 2 0, T, H 1 (Ω) . Finaly, using (3) to obtain the assertion of Lemma 3.2.
Now, to estimate the approximate solution, we state this following lemma: Proof. Combining (19) with the estimates (16) of Lemma 3.2, we have where C is a constant not depending on m. Now, taking φ ∈ (E A (Q)) N and ϕ ∈ (L 2 (Q)) N satisfying ||φ|| A,Q = 1 and ||ϕ|| L 2 (Q) = 1, then we have Using (22), we easily see that Therefore, we have a γ (., ∇u m ) is bounded sequence in L A (Q) and so ∂u m ∂t is a bounded in W −1,x L A (Q)+L 2 (Q) and similarly one can prove that ∂w m ∂t is bounded in L 2 0, T ; (H 1 (Ω)) , which complete the proof of lemma.

Third
Step. Passage to the limit. Thanks to the Lemma 3.2, there exist a subsequences of (u m , w m ) will be noted also by (u m , w m ) such that u m u weakly in W 1,x L A (Q) for σ( L A , E A ) and w m w weakly in L 2 0, T ; H 1 (Ω) .
Using lemma 3.3, there exists χ ∈ L A (Ω) N such that for a subsequence ∇w m → ∇w weakly in L 2 (Q) and Consequently, by applying the corollary 1, we get u m → u and w m → w strongly in L 1 (Q).

4.
Numerical aspects and results. In this section, we present the numerical performance of the proposed model applying to image denoising, in addition to comparative study to some denoising technique described in the literature. For simplicity, the problem (4) can be rewritten as follows: where with λ = e log γ (1 + ε) , the function f describes the noisy image and u(., t) is the image with scale parameter t.
For computing numerically problem (36), we attempt to discretize it by finite difference method. Assuming k to be the time step size and h the spatial grid size, we discretize time and space as follows: t n = nk, n = 0, 1, 2, ... where M h × N h is the size of the original image. Denote u n i,j and w n i,j the approximations of u(t n , x i , y j ) and w(t n , x i , y j ) respectively. We define the discrete approximation: We use the following notations for simplicity: We denote by div(p) n i,j the approximations of div(p(t n , x i , y j )) defined by: The discrete explicit scheme of the problem (36) could be written as: Through the above lines, we can obtain u 1 i,j and w 1 i,j by u 0 i,j and w 0 i,j . The program will stop when it achieves our goal. Most algorithm parameters are chosen heuristically for the algorithms to perform their best. We set λ = 0.01, the time step size k = 0.1 and the space step size h = 1. We start by the improvements tests (cf. Figs. 1-4) in the restorations provided by our approach and we choose the parameter γ = 10 −6 . Fig(4) exhibit an example test for the second operator by taking γ = 10 −6 , α(x) = 2.00001 and ε = 1. In the second experiment, we illustrate the influence of γ with respect to the restored image quality(see Fig 5). We remark that when γ is near to 0, the restored image is well improved. Fig (6) illustrate another experiment result of our proposed method compared with Total Variation model. In (Fig 6b), the image restored by the TV-based diffusion reconstructs sharp edges, but the staircasing phenomenon is clearly present ( nose, mouth and cheeks). In our case, we remark that our model reconstructs sharp edges as effectively as TV model and recovers smooth regions and reducing the staircasing effect (Fig 6a).       Figure 6. In this experiment, we present the comparison between restored images with our model and TV model. Prominent staircasing can be observed on Lena's face in the image obtained from the ROF model (b). Staircasing has been successfully reduced in the image obtained from our model.  Figure 7. Here, we compare the texture of restored images between our model and TV model. This results shows that texture and fine details are better preserved using the proposed framework than when using TV model.

5.
Conclusion. This paper describes a model for filtering gray scale images corrupted by independently and identically distributed Gaussian noise. The proposed model, which is based on reaction-diffusion system, reconstructs sharp edges as a TV model, preserves a fine details and reduces the staircasing phenomenon during the image denoising.