Rigorous Mathematical Investigation of a Nonlocal and Nonlinear Second-Order Anisotropic Reaction-Diffusion Model: Applications on Image Segmentation

In this paper we are addressing two main topics, as follows. First, a rigorous qualitative study is elaborated for a second-order parabolic problem, equipped with nonlinear anisotropic diffusion and cubic nonlinear reaction, as well as non-homogeneous Cauchy-Neumann boundary conditions. Under certain assumptions on the input data: f (t, x), w(t, x) and v0(x), we prove the well-posedness (the existence, a priori estimates, regularity, uniqueness) of a solution in the Sobolev space W1,2 p (Q), facilitating for the present model to be a more complete description of certain classes of physical phenomena. The second topic refers to the construction of two numerical schemes in order to approximate the solution of a particular mathematical model (local and nonlocal case). To illustrate the effectiveness of the new mathematical model, we present some numerical experiments by applying the model to image segmentation tasks.


Introduction
For the unknown function v(t, x) (hereafter, v), consider the following nonlinear second-order boundary value problem in Q = (0, T] × Ω, with T > 0 and a bounded domain Ω ⊂ I R 2 of Lebesgue measures |Ω|, whose boundary ∂Ω is sufficiently smooth: where: • t ∈ (0, T], x = (x 1 , x 2 ) varies in Ω, Σ = (0, T] × ∂Ω; Setting • ∆v(t, x) is the Laplace operator-a second-order differential operator, defined as the divergence (∇·) of the gradient of v(t, x) in x; • ∂ ∂t v(t, x) is the partial derivative of v(t, x) with respect to t; • p 1 , p 3 , p 4 , p 5 are positive values. • Φ v x (t, x) is a positive and bounded nonlinear real function of class C 1 (Q) with bounded derivatives (see [1]), having the role of controlling the speed of the diffusion process and enhances the edges (e.g., in the evolving image); • Ψ(v x (t, x)) is the mobility; •q(t, x) is a positive and bounded real function; • f (t, x) ∈ L p (Q) is the distributed control (a given function), where • w(t, x) ∈ W 1− 1 2p ,2− 1 p p (Σ) is the boundary control (a given function); • n = n(x) is the outward unit normal vector to Ω at a point x ∈ ∂Ω. ∂ ∂n denotes differentiation along n; Let us note Then, it is easy to recognize Equation (1) 1 as being quasi-linear with while the boundary conditions (1) 2 are of second type: (see [1] and reference therein).
For the reader's benefit, we write problem (1) in the equivalent form Concerning Equation (5) 1 , we recall that it is of quasi-linear type with principal part in divergence form (see [1]), with a i , i = 1, 2, given by (4) and a(t, x, v(t, x), v x (t, x)) = −p 3 v(t, x) − v 3 (t, x) − p 4 f (t, x).
In our paper, we study the solvability of the problems (1) in the class W 1,2 p (Q), characterized by the presence of some new physical parameters (p 1 , p 3 , p 4 , p 5 , Φ v x (t, x) , Ψ(v x (t, x))), the principal part being in divergence form and by considering the cubic nonlinearity p 3 v(t, x) − v 3 (t, x) , satisfying for n ∈ {1, 2, 3} the assumption H 0 in [21], that is: In Theorem 1, we prove the existence, regularity and uniqueness of solution for (1). (see [15] for a numerical study of Equation (1) corresponding to a linear reaction term v(t, x) − v 0 (x), with homogeneous Neumann boundary condition).
In the following we will denote by C several positive constants.

Well-Posedness of the Solution of (5)
Theorem 1 of this section presents the dependence of the solution v(t, x) of (5) on f (t, x) and w(t, x). In our study, we rely on the following:
In addition, we use the set C 1,2 (Q) (C 1,2 (Q)) of all continuous functions inQ (in Q) having continuous derivatives u t , u x and u xx inQ (in Q), as well as the Sobolev spaces W l p (Ω), W l,l/2 p (Σ) with non-integral l for the initial and boundary conditions, respectively (see [1]).
The main result for the study of the existence, a priori estimates, uniqueness and regularity for the solution of (1) (or (5)) is the next theorem.  |v(t, x)| < M for any (t, x) ∈ Q and for any z(t, x), the map Ψ(z(t, x)) is continuous, differentiable in x, its x-derivatives are measurable bounded, satisfies (6) and is a positive and bounded nonlinear real function of class C 1 (Q) with bounded derivatives and 0 < m 1 ≤ Φ v x (t, x) ≤ M 1 . In addition, for every ε > 0, the functions v(t, x) and Ψ(v x (t, x)) satisfy the relations Then, ∀ f ∈ L p (Q) and ∀v 0 ∈ W 2− 2 p p (Ω), with p = 3 2 , the problem (5) has a solution v ∈ W 1,2 p (Q) and the next estimate holds: where the constant C > 0 does not depend on v, f and w.

The Proof of Theorem 1
To prove this theorem, we use the Leray-Schauder principle. Thus, we consider the Banach space and a nonlinear operator H : where v(u, λ) is the unique solution to the next problem with We shall prove now the following technical lemma Lemma 1. We assume Hypotheses I 1 and I 2 to be valid. Then Proof. Indeed, since u ∈ L 3p (Q), then u L 3p (Q) ≤ Konst and thus i.e., the nonlinear term in (16) [1]). Next, from (10) it is easy to conclude that Thus, to prove that we have to prove that u 2 . Finally, we recall that f (t, x) ∈ L p (Q) and, owing to the above, we easy derive that the statement expressed by (16) is true.

The Proof of Theorem 1 (Continued)
Let us show that the nonlinear operator H(u, λ) defined by (14) satisfies the following Properties A and B.
A. If (15) has a unique solution, then H is well-defined. By the right hand of (15) 1 , using x) ∈ L p (Q) and thus, the same reasoning as in [1] allows us to conclude that , the linear parabolic boundary value problem formulated in (15) has a unique solution, that is (see (14)) v = H(u, λ) ∈ W 1,2 p (Q), ∀u ∈ B and ∀λ ∈ [0, 1]. Next, the embedding W 1,2 (3) and (7)), allows us to conclude that Thus, the operator H is well-defined. B. Let us now show that H is continuous and compact. The sketch of the proof is the same as in [1,15]. However, for reader convenience, we present details in the sequel.
The right-hand side in (17) belongs to L p (Q), since v n,λ n ∈ W 1,2 p (Q). Therefore, the L p -theory of PDE gives the estimate V n,λ n ,λ Owing to Lemma 1 we can derive that (u n ) 3 is bounded in L p (Q), ∀u n ∈ W 0,1 p (Q) ∩ L 3p (Q). In addition, the inequality (10), the working Hypothesis I 2 and the inclusion u n,λ n x i x j ∈ L p (Q), imply the boundedness in L p (Q) of the terms A(t, x, u n , u n x i ) and , it results that the remaining terms on the right-hand side from the above inequality are also bounded in To evaluate the difference H(v n , λ) − H(v, λ), we use again the relations (14), (15), and we obtain where The L p -theory applied to (19), gives us the estimate with a new constant C. From the convergence u n → u in W 0,1 p (Q) ∩ L 3p (Q) and the continuity of the Nemytskij operator (see [19] and references therein), as well as the Making use of the relations (18) and (20), we show the continuity of the nonlinear operator H defined by (14). Moreover, H is compact. Indeed, since µ > 3p, the inclusion is compact (see [12] and reference therein). Furthermore, writing H as the composition the compactness of H immediately follows.

The Proof of the First Part in Theorem 1: The Regularity of v(t, x)
We establish now the existence of a number δ > 0 such that The (see (4), (6) and (15)).
Multiplying the first equation in (22) by |v| 3p−4 v, integrating over Q t := (0, t) × Ω, t ∈ (0, T], we get Owing to Green's first identity, the left inequality in (9) and (12), Assumption I 2 and the boundary conditions (22) 2 , the previous equality leads us to for all t ∈ (0, T]. The Hölder and Cauchy inequalities, applied to the last terms in (23), give us By H 0 , relation (3) and Young's inequality, we obtain Owing to the above inequality as well as (i 1 -i 3 ) and, taking into account the continuous embedding L 3p−2 (Σ t ) ⊂ L p (Σ t ), from (23), we derive the following estimate Taking ε small enough, the previous inequality yields for a positive constant C 1 = C(|Ω|, T, n, p, p 1 , p 3 , p 4 , p 5 , q min , q max , m 1 , M 1 ).
By Lemma 1.1 in [21] and (24), we get and then (25) becomes The continuous embedding W 1,2 which, owing to (26), ensures that a constant δ > 0 can be found such that the property expressed in (21) is true. Denoting relation (21) implies that provided that δ > 0 is sufficiently large. Furthermore, following the same reasoning as in [1,4,11,15,19], we conclude that problem (6) has a solution v ∈ W 1,2 p (Q) (see also [21], p. 195). The estimate (11) results from (26), and the proof of the first part in Theorem 1 is finished.

The Uniqueness of the Solution v(t, x)
Now, we prove (13), which implies the uniqueness of the solution of (1) or (5). By hy-

2, and (following [1]) we write the increments of a ij in the form
where Now, we subtract Equation (1) 1 for v 2 (t, x) from Equation (1) 1 for v 1 (t, x), and making use of (27), (28), we obtain the following linear equation wherê Due to (9) and the working hypotheses on v 1 and v 2 , i.e., v 1 p (Q) ≤ M 4 , the conditions on linear equations are fulfilled and, given this, it follows from (29) that estimate (13) is valid for V, which finishes the proof of Theorem 1.
As a consequence, it results the uniqueness for the solution of (5).

Corollary 1.
For the same initial conditions, the problem (5) possesses a unique solution v(t, x) ∈ W 1,2 p (Q).

A Novel Nonlinear Second-Order Anisotropic Reaction-Diffusion Model in Image Segmentation
The nonlinear parabolic second-order PDE problem (5) can be applied for image denoising, enhancement, restoration and segmentation. Here we consider a particularization of this mathematical model by setting the functions Φ v x (t, x) and Ψ v x (t, x) as follow where ϕ, η, α ∈ (0, 4], while the parameter c is the conductance (see [15], p. 177 and [14], p. 633). Therefore, the following PDE scheme with non-homogeneous Cauchy-Neumann boundary conditions is acquired: The edge-stopping (diffusivity) function in (30) 2 is positive, monotonically decreasing and converging to zero (see [28,30]) thus satisfying the conditions imposed by a proper diffusion. Moreover, it is easy to check that Ψ and Φ in (30) satisfy Assumptions I 1 and I 2 in Theorem 1 and thus the nonlinear anisotropic reaction-diffusion model (31) is wellposed, as proved in the previous section. Consequently, it admits an unique classical solution v(t, x) ∈ W 1,2 p (Q), that represents the evolving image of the observed image v(0, x) = v 0 (x).
The corresponding nonlocal anisotropic reaction-diffusion model to (31) can be written as follows: where • K : I R → I R is a real function, symmetric, continuous, nonnegative and it's compactly supported in the unit sphere, such that Details on certain interpretations of the terms K(x − y),

y)dy and
−v(t, x) Ω K(x − y)dy in the mathematical model (32), can be found in the works of P. W.
Bates, S. Brown and J. Han [3] and J. Rubinstein and P. Sternberg [27] and references therein. The solution behavior for the nonlocal model (32) on rescaling the kernel K considering are studied in [33] and for the numerical solutions we refer to [3,40] and references therein.
In what follows, we will approximate the solution v(t, x) in (31) and (32) using the finite-difference method (of second-order in time, see (36)).
To approximate ∂ ∂t v(t, x), we employ a second-order scheme (see [16,41] and references therein): We write Equation in (32) as: where we denote the nonlocal diffusion term by: and the reaction term by: The left-side term in (37) is approximated by i,j and the right side terms are discretized using central differences (see [16] and references therein). We also denote Φ i,j = Φ( ∇v i,j ) and for all i = 2 . . . , I − 1, j = 2, . . . , J − 1. To complete the discretization schema we need to approximate NlD(t, x, v, v x ) and R(t, x, v, v x ) terms as follows: Continuing the discretization by using the Riemann sums to approximate the integral terms, we have: For the second integral on ∂Ω, we have: For the reaction term discretization, we use the following scalar product approximation: which leads to Further, since the second-order derivatives do not vary too much, we can use where v x1 = ∂v/∂x1, v x2 = ∂v/∂x2 and v x1x2 = ∂ 2 v/∂x1∂x2 are discretized by applying the finite difference method (see [15,28]). To conclude we obtain the following explicit numerical approximation for reaction term: and thus we get the following explicit numerical approximation scheme for (32): In a similar manner one obtains the following explicit numerical approximation scheme for (31):
The explicit numerical approximation scheme developed in (47) is consistent to the nonlinear second-order anisotropic reaction-diffusion model given by (32).
Some image segmentation results provided by our proposed model are displayed in Figures 1-4. All the results presented in this section are compared to standard K-means image segmentation model with two clusters [24] and the Chan-Vese image segmentation model presented in [5].
Our model successfully extracts the objects after up to three iterations. One may see multiple objects as well as objects with boundary concavities and blurry boundaries are accurately extracted from the background.

Algorithm 1: Reaction-diffusion based image segmentation algorithm
1 Set m = 1 2 Initialize the unknown function v 1 with the input image to be segmented 3 while v m did not reach stable state do 4 Compute diffusion and reaction terms according to (41), (42) and respectively (46) 5 Evolve function v m in (47) to obtain v m+1 i,j 6 Increase m by 1 7 Segmented image is given by v m Figure 1 shows the segmentation results of our model for a brain CT scan image. The results are satisfactory even after only one iteration. We also see the model reaching stability after two iterations in this case. Compared to K-means segmentation results, we observe the extracted objects edges (brain tissue and cranium bone) are better delimited from the background. Compared to Chan-Vese segmentation results, our model produces more accurate results too. In this example, Chan-Vese model seems to not follow the real object boundaries, especially at the border between cranium bone and brain tissue. Figure 2 shows the segmentation comparison between three cases: first the input image is segmented 'as is', second the input image is contaminated with noise before segmentation and third we double the noise added to the input image. For all three cases, we can also see the results of applying K-means and Chan-Vese segmentation. We see our model successfully removes most part of the noise in Figure 2h,l while still preserving a good approximation for the edges on the leaf object (better than both K-means and Chan-Vese).
In Figure 3, we see the segmentation results for a blurry boundary object as galaxy boundaries are slowly fading. Even after one iteration, our segmentation is superior to K-means and Chan-Vese as the real galaxy boundaries are correctly identified in Figure 3d. Figure 4 (virus microscopy) brings together noise, blur and irregular boundaries. Again, after two iterations, the model successfully identifies all objects of interest and the results, starting with the first iteration, are better than the compared K-means method. The Chan-Vese segmentation does not separate the virus blobs successfully, although it provides a good outer boundary approximation.
Regarding time complexity, due to the integral formulation of NlD term in (41) and (42), the proposed algorithm is slower than the compared K-means or Chan-Vese counterparts. To obtain better performance results, regarding running time, we had to implement the program on parallel architectures such as CUDA [42]. Table 1 shows the time taken by a CUDA implementation for different input image sizes (total number of pixels being I * J).
Using the local scheme in (48), we obtained promising results for image restoration tasks. Future work will show if we can succeed in mixing the local and nonlocal models for better noise removal before applying segmentation tasks.

Conclusions
The starting point in the elaboration of the present work is the paper by Miranville, A. and Moroşanu, C. [1], which is a major challenge for both theory and applications, focused on finding concrete cases of functions for the general case Φ(t, x, v(t, x), v x (t, x)) and Ψ(t, x, v(t, x), v x (t, x)) introduced in [1]. In this respect, a rigorous mathematical investigation is performed to analyze the well-posedness of the nonlinear anisotropic reaction-diffusion model (1) (in particular, (31)). The Leray-Schauder principle is applied to prove the existence and uniqueness of a unique classical solution v(t, x) ∈ W 1,2 p (Q), while the L p theory is used to derive the regularity properties for the solutions, considering that the initial data and the boundary constraints are compatible with the regularity and compatibility conditions (see (3)). In addition, the a priori estimates are made in L p (Q), which means the approximation for unknown functions v(t, x) are more precise (see [1,[11][12][13]15,[19][20][21]35]).
Using the finite-difference method (of second-order in time), two numerical schemes are constructed see (47) and (48) to approximate the solution v(t, x) of the new mathematical model. Numerical experiments show the model can be successfully applied to image segmentation tasks. We tested on images with multiple objects as well as objects with complex concavities or blurry boundaries and proved our model can accurately extract them, most of the time showing better results than the compared K-means model.
Summarizing, the main contributions in the present work are as follows: • We use novel techniques, such as Leray-Schauder principle, a priori estimates, L ptheory, to elaborate a rigorous qualitative study of the nonlocal and nonlinear secondorder anisotropic reaction-diffusion parabolic problem, endowed with a nonlinearity of cubic type as well as non-homogeneous Cauchy-Neumann boundary conditions, expressed by (1) and (31). We note that, due to the presence of the nonlinear coefficient Φ( v x (t, x)) (see (30)), the proposed second-order nonlinear reaction-diffusion scheme (31) represents a non-variational PDE model. Therefore, it cannot be obtained from a minimization of any energy cost functional, thus this scheme is not a variational PDE model. • Two two numerical schemes (47) and (48) are constructed to approximate the solution of the mathematical models (31) and (32) (local and nonlocal case).
Regarding the second theme, we aim to improve the scheme in (47) and (48), as part of our future research on the topic, by introducing new edge-stopping functions (see [28]) and by taking advantage of non-local image information which will allow us to apply the model to images with inhomogeneity (see [33] and reference therein).
The qualitative results obtained in this current work can be used in quantitative studies of the mathematical models in (1) or (5) as well as in the study of optimal control problems involving such nonlinear problems. We look forward to exploiting all these in our future works.