Anis Theljani

Zakaria Belhachmi
email: zakaria.belhachmi@uha.fr

A discrete approximation of Blake & Zisserman energy in image denoising and optimal choice of regularization parameters

Keywords:

We consider a multi-scale approach for the discrete approximation of a functional proposed by Bake and Zisserman (BZ) for solving image denoising and segmentation problems. The proposed method is based on simple and effective higher order variational model. It consists of building linear discrete energies family which Γ-converges to the non-linear BZ functional. The key point of the approach is the construction of the diffusion operators in the discrete energies within a finite element adaptive procedure which approximate in the Γ-convergence sense the initial energy including the singular parts. The resulting model preserves the singularities of the image and of its gradient while keeping a simple structure of the underlying PDEs, hence efficient numerical method for solving the problem under consideration. A new point to make this approach work is to deal with constrained optimization problems that we circumvent through a Lagrangian formulation. We present some numerical experiments to show that the proposed approach allows us to detect first and second-order singularities. We also consider and implement to enhance the algorithms and convergence properties, an augmented Lagrangian method using the alternating direction method of Multipliers (ADMM).

INTRODUCTION

Image denoising and segmentation problems have been intensively studied in the last decades and several approaches based on statistical methods, learning, wavelets or PDEs were developed with variable success have been registered, depending on the nature of the images under consideration (geometric features, textures, . . . ) [START_REF] Anis | High-order anisotropic diffusion operators in spaces of variable exponents and application to image inpainting and restoration problems[END_REF][START_REF] Peter | Total variation image restoration: numerical methods and extensions[END_REF][START_REF] Tony | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF][START_REF] Vese Luminita | Image denoising and decomposition with total variation minimization and oscillatory functions[END_REF][START_REF] Peter | Color TV: total variation methods for restoration of vector-valued images[END_REF][START_REF] Wenqi | Implementation of high-order variational models made easy for image processing[END_REF] . Both the difficulties and the advantages of each method are deeply related to a best handling of small singular sets (edges, corners,. . . ). In the PDEs community, a secondorder model, proposed first by Mumford and Shah (M-S) [START_REF] David | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] , turns out to be among the most successful. It consists of minimising the energy:

2 ∫ Ω∖Γ |∇ | 2 + 1 2 ∫ Ω 0 ( -) 2 +  1 ( ), (1)
where is the (reconstructed) image, the initial image and is the set of the contours of .  1 denotes the 1-dimensional Hausdorff measure is the size of , which generalizes the length for regular curves. The variational framework, allows for various approximations of the Mumford & Shah functional in both continuous and discrete settings (see [START_REF] Tortorelli | Approximation of functional depending on jumps by elliptic functional via Γ-convergence[END_REF][START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF][START_REF] Braides | Non-local approximation of the Mumford-Shah functional[END_REF] . In [START_REF] Tortorelli | Approximation of functional depending on jumps by elliptic functional via Γ-convergence[END_REF] , Ambrosio and Tortorelli proposed an elliptic approximation on the Sobolev spaces via the Γ-convergence. In their approach, the discontinuity set is approximated by an auxiliary variable which plays the role of an indicator function intended to take the value 0 at the edges and 1 in the rest. In [START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] , a discrete approximation of the Mumford-Shah energy were proposed and it is based on a finite element discretization and an adaptive mesh strategy.

However, the Mumford-Shah functional fails to capture the "discontinuities of second kinds", i.e., the set of gradient discontinuity. This limitation produces also some staircasing effect. In order to overcome such shortcomings of the Mumford-Shah functional, Blake and Zisserman (B-Z) [START_REF] Blake | Visual reconstruction[END_REF] proposed a new variational model (1987) which exhibits a second order derivatives by considering the following energy:

∫ Ω⧵( ∪ ∇ ) |∇ 2 | 2 +  -1 ( ∩ Ω) +  -1 (( ∇ ⧵ )) ∩ Ω) + 1 2 ∫ Ω 0 ( -) 2 . ( 2
)
The previous energy is defined in the class of functions of

2 (Ω) = { ∈ (Ω), ∇ ∈ (Ω)} .
It depends on the free discontinuities, free gradient discontinuities and second derivatives of (see [START_REF] Michele | A second order model in image segmentation: Blake & Zisserman functional[END_REF][START_REF] Michele | Strong minimizers of Blake & Zisserman functional[END_REF][START_REF] Andrea | A compactness result for a second-order variational discrete model[END_REF][START_REF] Luigi | Variational approximation of a second order free discontinuity problem in computer vision[END_REF] for more details).

The discontinuities of and of ∇ are "a priori "unknown, hence the associated minimization problem turns out to be essentially non trivial. A challenging problem is how to get effective and suitable numerical schemes for the computation of solutions for Blake-Zisserman energy? In the same spirit of the elliptic approximation of the Mumford-Shah energy, Ambrosio et al. (2001) proposed the following functional:

( ) = ∫ Ω 2 |∇ 2 | 2 + ∫ Ω ( 2 + )|∇ | 2 + ( -) ∫ Ω |∇ | 2 + ( -1) 2 4 + ∫ Ω |∇ | 2 + ( -1) 2 4 + 1 2 ∫ Ω 0 ( -) 2 , (3)
defined on Sobolev spaces. The minimization of the previous energy acts not only on the restored image but also on two auxiliary functions: which is a control function for ∇ and which is a control function of the Hessian of , i.e., ∇ 2 (see [START_REF] Blake | Visual reconstruction[END_REF][START_REF] Zanetti Massimo | Numerical minimization of a second-order functional for image segmentation[END_REF] ). Thus, the functional implies a minimization over the variables ( , , ) of a non-linear PDEs system. Besides, one of the limitation of this method is the dependency on threshold parameters, producing thickness around the edges. But a major shortcoming is the knowledge of the exact constants and the difference of which gives exactly the size ⧵ ∇ and which is difficult, if not impossible to have a priori.

CONTRIBUTION AND ORGANIZATION OF THE PAPER

In previous works, we have proposed another approach close to the one followed by [START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF] , based on the construction of a family of discrete energies which converge (in the Γ-convergence sense) to MS -like functional, i.e. with the length term  1 ( ), [START_REF] Anis | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF][START_REF] Zakaria | WEIGHTED HARMONIC AND COMPLEX GINZBURG-LANDAU EQUATIONS FOR GRAY VALUE IMAGE INPAINTING[END_REF] . The method acts essentially with two ingredients: a tight location of the singular set using a posteriori error indicators, and a decreasing of a diffusion function which inhibits the diffusion cross the edges. the method doesn't uses threshold parameters and the length of singular sets is controlled by the Lebesgue measure of the set where the diffusion value is minimal. In this article, we follow this idea and build a Γ-convergent family of discrete energies which aims to obtain, in the limit, the BZ-energy. We consider the following higher order energy

( ) = 1 2 ∫ Ω ( )|∇ 2 | 2 + 1 2 ∫ Ω ( )|∇ | 2 + 1 2 ∫ Ω | -| 2 + 2 2 ∫ Ω |∇ -∇ | 2 , ( 4
)
where 0 < ≤ , 0 < ≤ . ∇ 2 is the Hessian matrix of , 1 > 0, 2 ≥ 0. The first and the second part in the energy (4) are regularization terms. The third part is a data fidelity term which enforce the solution to be close to the initial image (in the 2 or 1 norms). We emphasize that under this expression the relevant information on the singularity set and ∇ is contained in the diffusion functions and , but the set ⧵ ∇ is contained in the set ∈ Ω; = Λ > . We notice also that taking 2 > 0, adds a term which is not in the initial BZ functional; this term is only added to enhance the gradient of .

In order to solve problem (4), we consider a Lagrangian approach , where we transform our problem in a constrained minimization one, by considering as a variables and = ∇ , this has also the advantage of leading to solve only second order problems. Thus, problem (4) takes the form

min ,  1 ( ) +  2 ( ), subject to + = 0,
where  1 (⋅) and  1 (⋅) are (second order) convex functions, and are linear operators (namely ∇ and the identity). Then, we make an alternating process to solve separately  1 (⋅) (with fixed and ) and  2 (⋅) (with fixed , ). In each step, we perform the adaptive strategy with and obtained following the method in [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF][START_REF] Anis | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF][START_REF] Zakaria | WEIGHTED HARMONIC AND COMPLEX GINZBURG-LANDAU EQUATIONS FOR GRAY VALUE IMAGE INPAINTING[END_REF] . The selection of these functions is performed locally and adaptively -using error indicators and an efficient mesh adaptation -, to decrease with explicit formulae their values close to edges and corners. These considerations lead to an overall approach that is numerically easy, low cost (in the sense that homogeneous regions are coarsed in the meshes), and efficient in determining the singular sets without a priori knowledge except that they have some length.

E-I

The paper is organized as follows: In section 2, we prove the existence of solution of the minimization problem (4) under a variational mixed formulation and we introduce the Lagrangian formulation. In section 3, we build a family of discrete approximation energies and perform the analysis in the framework of Γ-convergence approach. We give the main algorithm and the a posteriri error indicators which allows us to perform the diffusion coefficients selection that preserves the singular sets for the function and its gradient. In section 4, in order to improve the convergence rate of the algorithm we present an alternative approach to deal the higher order model ( 4) via an augmented Lagrangian formulation. In particular, we use an ADMM algorithm (alternating direction descent) to enhance the convergence of the Lagrange multiplier variable. In section 5, we present several numerical simulations to show that both singularities of first and second kind are preserved while the diffusion operators act as high order filtering in the homogeneous zones. We give some conclusions and perspectives to the proposed approach.

SPLITTING METHOD

The minimization of the energy (4) leads to a fourth-order partial differential equation. In this work, we consider the following constrained optimization problem:

1 2 ∫ Ω ( )|∇ | 2 + 1 2 ∫ Ω ( )|∇ | 2 + 1 2 ∫ Ω | -| 2 + 2 2 ∫ Ω | -∇ | 2 , ( 5
) . . ∇ = , (6)
Thus, the previous optimization problem splits as follows:

min ( , )∈ 1 (Ω)× 1 (Ω)  1 ( ) +  2 ( ), subject to (6) (7)
where

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩  1 ( ) = 1 2 ∫ Ω ( )|∇ | 2 + 1 2 ∫ Ω | -| 2  2 ( ) = 1 2 ∫ Ω ( )| ( )| 2 + 2 2 ∫ Ω | -∇ | 2 and = 1 2 = 1 2 , ( ) = 1 1 1 2 2 1 2 2 .
Proposition 1. The problem (5) admits a unique solution ( * , * ) ∈ 1 (Ω) × 1 (Ω)

Proof. The functional  (⋅, ⋅), as a sum of norms, is continuous, strictly convex on 1 (Ω)× 1 (Ω). In fact, let

( 1 , 1 ), ( 2 , 2 ) ∈ 1 (Ω) × 1 (Ω) such that ( 1 , 1 ) ≠ ( 2 , 2 )
and ∈]0, 1[, then we have:

 ( 1 , 1 ) + (1 -) ( 2 , 2 ) ≥ min( 2 , (1 -) 2 ) ⎡ ⎢ ⎢ ⎣ ∫ Ω ( )|∇( 1 -2 )| 2 + 1 ∫ Ω | 1 -2 | 2 ⎤ ⎥ ⎥ ⎦ + min( 2 , (1 -) 2 ) ⎡ ⎢ ⎢ ⎣ ∫ Ω ( )| ( 1 -2 )| 2 ⎤ ⎥ ⎥ ⎦ + 2 ∫ Ω | 1 -2 | 2 > 0. (8)
We also have that  (⋅, ⋅) is coercive and weakly lower semi continuous in 1 (Ω) × 1 (Ω). Moreover, the feasible set  = {( , ) ∈ 1 (Ω)× 1 (Ω), ∇ = } of the previous minimization problem is non empty, closed and convex set of 1 (Ω)× 1 (Ω). Thus, the existence and the uniqueness follows from classical results of calculus of variations [START_REF] Belhachmi | Mortar spectral element discretization of nonhomogeneous and anisotropic Laplace and Darcy equations[END_REF] Remark 1. Note that when 2 = 0, we modify in a standard way the space for , such that ∈ 1 ⧵ ℝ.

A Lagrangian algorithm

We introduce the space

0 (div, Ω) = ∈ ( 2 (Ω)) 2 ; div ∈ 2 (Ω); ⋅ = 0, on Ω .
We associate to the primal problem (2.1) -(2.2) the following Lagrangian form:

 1 ( , , ) =  1 ( ) +  2 ( ) -∫ Ω ⋅ + div , ( 9
)
We can show that the there exists a saddle-point ( * , * , * ) ∈ 1 × × 0 (div, Ω) of the Lagrangian  1 (⋅). In fact, let  be the dual functional of  1 (⋅), which is defined by the following concave functional:

( ) = inf ∈ 1 (Ω), ∈ 1 (Ω)  1 ( , , )
Then, the dual problem consists in solving the following maximization problem:

max ∈ 0 (div,Ω) ( ). (10)
Then, since the constraint convex minimization problem (4) has a unique optimal solution ( * ; * ), and the dual concave problem (10) also has an optimal solution * , it follows that ( * ; * ; * ) is a saddle point for  1 (⋅).

DISCRETE APPROXIMATION AND Γ-CONVERGENCE

In the sequel, we follows the adaptive method proposed in [START_REF] Anis | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF][START_REF] Belhachmi | Control of the effects of regularization on variational optic flow computations[END_REF] for the discrete approximation of Mumford-Shah energy (in the scalar case) and we explain how it can be extended to our case for the the variables , and . Notice that the Lagrangian is not convex in and thus this part of the energy should be treated separately within the data term.

Γ-convergence

In [START_REF] Tortorelli | On the approximation of free discontinuity problems[END_REF] , the authors proposed an adaptive approach which is based on two ingredients [START_REF] Belhachmi | Control of the effects of regularization on variational optic flow computations[END_REF] : the usual mesh adaptation and a "functional" adaptation which consists in choosing locally the diffusion coefficients in order to "cut" high gradients of the computed solution. This approach turns out to be a well suited finite element approximation with a family of discrete energies, in the Γ-convergence sense [START_REF] Braides | Gamma-Convergence for Beginners[END_REF] to the Mumford-Shah functional [START_REF] David | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] (see [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] for more details). For a fixed angle 0 < 0 ≤ 2 ∕3, a constant ≥ 6, and for > 0, let  (Ω) =  (Ω; 0 ; ) be the set of all triangulations of Ω whose triangles have the following characteristics:

(i) The length of each of the three edges of is between and .

(ii) The three angles of are greater than or equal to 0 .

Let (Ω) be the set of all continuous functions ∶ Ω ←→ ℝ such that is affine on each triangle of a triangulation ∈  (Ω). For a given , we define  ( ) ⊂  (Ω) as the set of all triangulations adapted to the function , i.e., such that is piecewise affine on . We consider a non-decreasing continuous function ∶ [0, +∞) ←→ [0, +∞), ( = 1, 2) such that:

lim →0 ( ) = 1, lim →+∞ ( ) = ,∞ < +∞.
(one may take for such a concave regularization of ( ) = min( , ∞ )). For any = ( , ) ∈ ( (Ω)) [START_REF] Tony | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] , ( > 1) and ∈  (Ω), following [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF][START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF] we consider the finite dimensional minimization problem:

( ) = min ∈ (Ω) ̃ ( , ), (11)
where

̃ ( , ) = ⎧ ⎪ ⎨ ⎪ ⎩ ∑ ∈ | ∩ Ω| 1 ℎ 1 (ℎ |∇ | 2 ) + | ∩ Ω| 1 ℎ 2 (ℎ |∇ | 2 ), ∈ ( (Ω)) 3 , ∈  (Ω), +∞, otherwise.
When goes to zero and provided 0 is less than some Θ > 0, we have that converges, in the Γ-convergence sense to the functional given by:

( ) = ⎧ ⎪ ⎨ ⎪ ⎩ ∫ Ω |∇ ( )| 2 + |∇ ( )| 2 + 1,∞  1 ( ) + 2,∞  1 ( ), ∈ ( 2 (Ω) ∩ (Ω)) 3 , +∞, ∈ ( 2 (Ω)∖ (Ω)) 3 .
where (Ω) is the space of generalized special functions of bounded variation (see [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF][START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF] for details). This means, in particular, that minimizers of converge to a minimizer of .

Remark 2. Notice the following important property that we will use next. If ∶ ←→ [-∞, +∞] is continuous and ( ) Γ-converges to , then ( + ) Γ-converges to + . Remark 3. Notice also that at this stage, the Γ-limit is different than BZ energy that we want to approximate as the singular sets and have no reason to be disjoints.

In the sequel, we better clarify the iterative algorithm used in this work (which mimic the proof of Γ-convergence), and we will see the role and the way to "compute" the diffusion coefficients and . For this purpose, let be the Legendre-Fenchel transform of ( , ) = 1 ( ) + 2 ( ) (concave for both variables). Then, we have ( , ) = 1 ( ) + 2 ( ) where 1 and 2 are the Legendre-Fenchel transforms of -( 1 ) and -( 2 ), respectively. For a given triangulation , it is readily checked (from convex analysis arguments) that the minimization of is equivalent to the minimization problem

′ ( , 1 , ) = ∑ ∈ | ∩ Ω| |∇ | 2 + 1 ( ) ℎ + ∑ ∈ | ∩ Ω| |∇ | 2 + 2 ( ) ℎ ,
over all = ( , ) and 1 = ( , ) where ∈ (Ω) and = ( ) ∈ (Ω) , ∈ ( (Ω)) 2 and = ( ) ∈ (Ω) , piecewise constant on each ∈  . Now, the minimization over each and which may be explicitly computed by:

= ′ 1 (ℎ |∇ | 2 ) and = ′ 2 (ℎ |∇ | 2 ). (12)
Then, the choice of and depends on the functions 1 and 2 and may be computed analytically without solving any PDE. An example of such choices could be given by :

( ) = 2 ,∞ arctan ( 2 ,∞ ),
which produces smoothly decreasing sequences and (as the inverse of the magnitude of the respective gradients squared). Whereas, the minimizers ( , ) solve, up to the data terms,

( , ) = ∑ ∈ ∫ ∩Ω |∇ | 2 + ∑ ∈ ∫ ∩Ω |∇ | 2 (13)
To complete the Γ-convergence analysis, let us denote the fidelity term of the Lagrangian by

Φ( , ; ) = 1 2 ∫ Ω | -| 2 + ∫ Ω | -∇ | 2 -∫ Ω ⋅ + div .
It is readily checked that Φ is continuous in ( 2 (Ω)) [START_REF] Tony | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] × 0 (div, Ω), hence using remark -2, and Γ-convergence of the sequence to , yields:

Proposition 2. Let be a given function in 2 (Ω) and > 0 be a positive parameter. Then the sequence of functionals ( (⋅) + Φ(⋅, ⋅; )) , Γ-converges for ←→ 0. ( (⋅) + Φ(⋅, ⋅; )).

Remark 4. The Lagrange multiplier in the algorithm makes the data term also dependent on . However, due to the continuity of Φ, for this third variable, in (div, Ω) this dependance does not affect the Γ-convergence of + Φ. More we precisely, we write Φ(., ), as Φ(., ) plus (Φ(., ) -Φ(., )) and use the strong convergence of the last term and the remark 2.

To conclude this section, let us interpret in term of a PDE the Γ-limit of each functional , i.e., corresponding to and . The minimizers ( , ) converge to a minimum ( , ) of while ( ), respectively ( ), converges to 1 Ω⧵ , respectively 1 Ω⧵ . The singular sets are identified as the sets where lim ←→0 = 0 for , respectively lim ←→0 = 0 for . Notice that the two singular sets are not necessarily disjoints. However, the Hausdorff measure of each set is bounded (tightly) by the Lebesgue measure of the corresponding sets = , respectively = . Moreover, the Hausdorff measure of the set ∇ ⧵ is then controlled by | > , and = |. The constants in the BZ functional (or its elliptic regularization) are obtained directly from the data and not chosen a priori.

Alternate descent algorithm and adaptive approach

Following the method in [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF][START_REF] Anis | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF] , we consider an alternate descent minimization scheme over ( , , ), , and which reads:

′ ( , 1 , ) + Φ( , ; ), (14)
The alternate minimization is combined with an adaptive algorithm to approximate B-Z energy, see also [START_REF] Belhachmi | Control of the effects of regularization on variational optic flow computations[END_REF] . Mainly, we use mesh and "functional" adaptation which consists in choosing locally the diffusion coefficients and in order to "cut" high gradients of the computed solution and its gradient . More precisely, we consider the following algorithm Algorithm 1 MINIMIZATION METHOD:

1. Initialization: choose 0 = ( 0 , 0 , 0 ), 0 and 0 .

2. Iterations: for fixed and , find +1 = ( +1 , +1 ) ∈ (Ω) [START_REF] Tony | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] and +1 ∈ (Ω) [START_REF] Peter | Total variation image restoration: numerical methods and extensions[END_REF] , solution of problem (13).

with

+1 and +1 fixed: a) Update with mesh refinement step (in our case using the metric adaptation with FreeFem++).

b) Perform a local choice of ( ) and ( ) on  +1 ℎ to obtain new functions +1 and +1 . (using a posteriori error indicators).

Go to step (1) until convergence.

The step ), plays the role of minimization over a class of triangulations  , and allows us to obtain a tight location of the singular sets as the residual error indicator is high (or equivalently with the metric adaptation). This means exactly that this locations are the candidate to contain the singularities. We emphasize that we obtain also an approximate solution to problem (14).

For a given and and , the step b) will consists in decreasing the diffusion parameters. Contrary to the scalar case, when a part of the singular set is located with the step ) we don't know whether this part is in or or both. Thus, to perform this step, we first consider the following local error indicators

= 1 2 ‖∇ ℎ ‖ 2 ( ) and ′ = 1 2 ‖∇ ℎ ‖ 2 ( ) , ( 15
)
to decide which parameter ( , or both) we decrease. The change in the intensity of the image, respectively the brightness, indicate the set , respectively . Thus we may know which parameter to change. The map of the above errors indicator allows us to update locally and on each triangle using either the formula (12), or the following one (fast decreasing) [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF][START_REF] Anis | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF] and which is given by

+1 = max ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1 + ‖ ‖ ∞ -0.1 + , ℎ ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
, where ℎ is a threshold parameter and is a coefficient chosen to control the rate of decrease in , + = max( , 0). The same formula is used to update the function ( ) and where is replaced by ′ . Notice that equivalently, we can use the residual errors indicators related to the discrete and subproblems.

Let us gives the details for the step 2. of the algorithm

The -subproblem For , and fixed, the -subproblem is quadratic and consists in solving the following minimization problem:

min ∈ 1 (Ω) { ∑ ∈ ∫ ∩Ω |∇ | 2 + Φ( , ; )} (16)
The minimizer +1 satisfies the following first-order optimality condition

-Δ +1 + 1 +1 = 1 -div , in Ω, +1 = 0, in Ω,
The -subproblem For a fixed +1 , and , the -subproblem consists in solving the following minimization problem:

min ∈ 1 (Ω) { ∑ ∈ ∫ ∩Ω |∇ | 2 + Φ( , ; )} (17)
which clearly admits a unique solution +1 which fulfils:

-Δ +1 + 2 +1 = 2 ∇ -, in Ω, +1 = 0, in Ω,

The -subproblem

The multiplier updates linearly as follows: for a given parameter (small)

+1 = + ( +1 -∇ +1 ). (18)
Remark 5.

1. The existence of solutions of -and -problems follows from standard techniques of calculus of variations. In fact, it is readily checked by direct calculation that the functional to be minimized are strictly convex, coercive and weakly lower semi-continuous.

2. Notice that problems ( 16) and ( 17) may be solved in parallel to accelerate the convergence.

3. The role of 2 > 0 appears clearly at this step to deal with a "non degenerate" Neumann problem.

AUGMENTED LAGRANGIAN FORMULATION AND ALTERNATING DIRECTIONS METHOD OF MULTIPLIERS (ADMM)

In order to enhance the quality and the cost of the numerical computations, we use the Augmented Lagrangian formulation, which will be solved using the alternating directions method of multipliers (ADMM). The ADMM is one of the most extensively used algorithms to solve constrained optimization problems, particularly for augmented Lagrangian formulations (see [START_REF] Stephen | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF][START_REF] Yu | Global convergence of ADMM in nonconvex nonsmooth optimization[END_REF][START_REF] Roland | Lectures on numerical methods for non-linear variational problems[END_REF][START_REF] Tom | Fast alternating direction optimization methods[END_REF] and the references therein).

In our case, we associate to the primal problem (2.1) -(2.2) the following augmented Lagrangian form:

 ( , , ) = ∑ ∈ ∫ ∩Ω |∇ | 2 + |∇ | 2 -⋅ + div + 2 | -∇ | 2 (19) = ∑ ∈ ∫ ∩Ω ( + 2 )|∇ | 2 + |∇ | 2 + 2 2 -⋅ ( + ∇ ) + div (20)
for > 0 (usually = 1). The ADMM iterations consists in minimizing  (⋅, ⋅, ) with respect to one variable and keeping the other fixed. The Lagrange multiplier is up dated at each step in the algorithm and takes the same form (18) as for the classical Lagrangian method. The only change comparing to Algorithm is the and subproblems. The overall algorithm consists in the following steps:

Algorithm 2 ADMM ALGORITHM Given , , , .

Starting guess: 0 ∈ 1 (Ω), 0 ∈ 0 (div Ω) and 0 ∈ 1 (Ω). For = 0, 1, 2, ⋯ until convergence:

1. Update +1 by +1 = arg min  ( , , ).

2. Update +1 by solving +1 = arg min  ( , , ).

3. Update +1 by solving +1 = + ( +1 -∇ +1 ).

In the new formulation, problems ( 16) and ( 17) are replaced by

⎧ ⎪ ⎨ ⎪ ⎩ -Δ + 2 +1 + 1 +1 = 1 -div -div , in Ω, +1 = 0, in Ω, and -Δ +1 + 2 +1 + +1 = 2 ∇ + ∇ +1 -, in Ω, +1 = 0, in Ω.

The convergence of the ADMM algorithm

There are many convergence results for ADMM discussed in the literature [START_REF] Eckstein | Augmented Lagrangian and alternating direction methods for convex optimization: A tutorial and some illustrative computational results[END_REF][START_REF] Robert | A general analysis of the convergence of ADMM[END_REF][START_REF] Stanley | Plug-and-play ADMM for image restoration: Fixed-point convergence and applications[END_REF] . The convergence is guaranteed if the following assumptions are satisfied.

• Energies  1 (⋅) and  2 (⋅) are convex, closed and proper.

• The augmented Lagrangian  admits a saddle point ( * , * , * ) ∈ 1 × × 0 (div, Ω).

The first assumption is satisfied as each of the -and -subproblems admits a unique solution. For the second assumptions, we can show that the there exists a unique saddle-point ( * , * , * ) ∈ 1 × × 0 (div, Ω) of the augmented Lagrangian  . In fact, any saddle-point of the augmented Lagrangian  (⋅) is also a saddle-point of the Lagrangian  1 (⋅) in (9), which has a unique saddle point. Remark 6. The use of the ADMM seems to enforce the diffusivity on the function , which may cause loosing its singularities set. Moreover, the analysis of the Γ-convergence is questionable in the new formulation. To keep the adaptive algorithm able to detect the singularities sets of , we set a new diffusivity function

( ) = ( ) +
2 . Then, we adapt the function using the same formula and the we rescale it into the interval [ -2 , -2 ]. Thus, we will be able to detect singularities from the adaptation and the new diffusivity fulfils 0 < ≤ ( ) ≤ , which guarantees the well-posedness of the -problem.

NUMERICAL EXAMPLES

In this section, numerical experiments are conducted on images that are corrupted by Gaussian noise. We test the proposed model for first-and second-order discontinuities detection. For the denoising task, we compare our approach with the total variation model (TV) [START_REF] Vese Luminita | Image denoising and decomposition with total variation minimization and oscillatory functions[END_REF][START_REF] Peter | Total variation image restoration: numerical methods and extensions[END_REF] and a high-order model that we call HOM, which corresponds the energy (4) for (⋅) ≡ ≡ (⋅) = 1. The original image were corrupted with Gaussian noise having zero mean and different standard deviations, i.e., between 10 and 15.

Edge detection

In Fig 1 , we test proposed model for ≠ 0 and = 0 for edge detection. We note that the case = 0 corresponds to the discrete approximation of the Mumford-Shah energy [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF][START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF] . The values of the functions and are 1 (white) on homogeneous areas and they are 0 (black) in correspondence of the edges of and ≈ ∇ . It is clear the map is unable to detect the gradient discontinuity. These are regions where is continuous while ∇ is discontinuous. The output map of our approach is given by meaningful boundaries corresponding to the discontinuity set of its first derivatives.

Lagrangian vs Augmented Lagrangian

In Figure 2 , we show the result of the proposed model, solved using the Lagrangian and ADMM, respectively. for both approaches, We display the denoised images, the Lagrange multiplier = ( 1 , 2 ) and the final mesh. The results of both approaches look visually similar, however, the ADMM is much faster for solving the model ( 5)- (6).

Image denoising and edges detection

In Figure 3 , we test the proposed model, solved by ADMM, in denoising a satellite image which has some fine structures, correspond to the roads. We show the denoised image, the maps of and . The function maps the discontinuity of and the function which maps the discontinuities of , i.e., approximates the gradient of . Form the map of , it clear that latter suitable in detecting fine structures. We also test the model in denoising a medical image, which contains some smooth regions and edges. In these regions, the gradient of the image should discontinuous and they can seen in the map of .

Medical image denoising and edges detection

In the last example, we give in Figure 4 the results obtained with the proposed model solved with the Lagrangian and with ADMM for medical images. We present the plot of and and the final mesh. This example show the efficiency of the method for medical images where in addition to the variety of fine structures to be detected, the high resolution of the images should also be considered.

We compare the proposed model with the total variation (TV) and the high-order (HOM) models. In Figure 5 , we display the restored image we give the PSNR values, which confirm that the proposed model can deliver better denoising quality than the other models. We also display the diffusion functions and which indicates regions that approximate the discontinuity sets of and ∇ , i.e., and ∇ , respectively.

CONCLUSION

In this paper, we proposed a novel variational model with the good performance such as the alleviation of the staircase effect and the preservation of the edges. By using an alternative approach, we obtain two simple subproblems, acting for alternatively for the solution and its first order derivatives ∇ . To solve these subproblems, the ADMM algorithm was used in the numerical computation. The use of the ADMM allowed as to see the initial model, which is of high order, as two first-order variational models with unknowns and ∇ . Then, we performed an adaptive strategy based on the objective selection of the regularization parameters appeared in both energies and mesh adaptation techniques. Analysing the proposed algorithm using Γ-convergence tools draw connections with the Blake and Zisserman energy. The numerical experiments have demonstrated the effectiveness of our proposed method.

FIGURE 1 (FIGURE 2

 12 FIGURE 1Example of first-and second-order discontinuities detection using our model

FIGURE 3

 3 FIGURE 3 Example of denoising and detection of thin structures (roads) using our model. The coefficient (⋅) is clearly more efficient than (⋅) for this task.

  (a) The noisy image (b) The restored image using our model (c) The map of (⋅) (d) The map of (⋅) (e) The adapted mesh

FIGURE 4

 4 FIGURE 4 Our model solved by Augmented Lagrangian and ADMM: Example of medical image denoising and edges detection.

FIGURE 5

 5 FIGURE 5 Medical image denoising and edges detection: Comparison between Total variation model (TV), the high-order model (LLT) and our model.


Anis Theljani
email: thaljania-nis@gmail.com

Zakaria Belhachmi
email: za-karia.belhachmi@uha.fr

Maher Moez Kallel
email: mz.kallel@gmail.com

Maher Moakher
email: maher.moakher@gmail.com

Combined second and fourth-order PDEs model and associated variational problems for geometry images inpainting and denoising

Keywords: Image processing, Inverse problems, Regularization, Γ-convergence, PDEs AMS subject classifications. 35J50, 35Q68, 49J45, 65M10, 65N30

published or not. The documents may come

Combined second and fourth-order PDEs model and associated variational problems for geometry images inpainting and denoising

Introduction.

Image inpainting is a central problem in image processing. It refers to the problem of recovering damaged/missing parts of a digital image by interpolation from the valuable areas [START_REF] Bertalmio | Inpainting, Encyclopedia of Computer Vision[END_REF][START_REF] Cao | Geometrically guided exemplar-based inpainting[END_REF][START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF][START_REF] Bertalmio | Image inpainting[END_REF] For the geometry inpainting, Partial Differential Equations and energy methods are extensively used and turn out to be very efficient (see [START_REF] Ambrosio | A direct variational approach to a problem arising in image reconstruction[END_REF][START_REF] Bertalmio | Navier-Stokes, fluid dynamics, and image and video inpainting, reconstruction[END_REF]7,8,[START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF][START_REF] Esedoglu | Digital image inpainting by Mumford-Shah-Euler model[END_REF][START_REF] Grossauer | Using the complex Ginzburg-Landau equation for digital inpainting in 2D and 3D, Scale Space Methods in Computer Vision[END_REF][START_REF] Kallel | The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting[END_REF][START_REF] Masnou | Level-lines based disocculusion[END_REF][START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF] and the references therein). The main difficulty in geometry inpainting is to accurately recover the singularity set components (edges, corners, . . . ) and to preserve the curvy features of the level lines of the image. When there is no missing area, the problem reduces to a segmentation/filtering problem and the same difficulties may persist. Cconceptually, the main difference between this two situations is that available informations come only from one side in the inpainting case.

Second-order PDEs are used in this field since its beginning leading to well established models, such as TV inpainting [START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF], anisotropic diffusion [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF] and (weighted) harmonic method [START_REF] Belhachmi | Weighted harmonic and complex Ginzburg-Landau equations for gray value image inpainting[END_REF], that provide good interpolation properties in geometry inpainting. However, as they are low-order methods, they usually disconnect edges over large distances (violating the connectivity principle) and they fail to reproduce some geometric features of higher-order (curvature, some corners, see e.g., [START_REF] Allard | Total variation regularization for image denoising, III. examples[END_REF][START_REF] Shah | A common framework for curve evolution, segmentation and anisotropic diffusion[END_REF]). These shortcomings gave rise to a new class of higher-order diffusion models which perform generally better both in image restoration and in geometry inpainting problems (see [START_REF] Bertalmio | Image inpainting[END_REF][START_REF] Ambrosio | A direct variational approach to a problem arising in image reconstruction[END_REF]8,[START_REF] Esedoglu | Digital image inpainting by Mumford-Shah-Euler model[END_REF][START_REF] Masnou | Level-lines based disocculusion[END_REF][START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF][START_REF] Shah | A common framework for curve evolution, segmentation and anisotropic diffusion[END_REF][START_REF] Bertozzi | Inpainting of binary images using the Cahn-Hilliard equation[END_REF]9,[START_REF] Cherfils | Finite-dimensional attractors for the Bertozzi-Esedoglu-Gillette-Cahn-Hilliard[END_REF][START_REF] Cherfils | A Cahn-Hilliard system with a fidelity term for color image inpainting[END_REF][START_REF] Chung | A fast algorithm for Euler's elastica model using augmented Lagrangian method[END_REF][START_REF] Scherzer | Denoising with higher order derivatives of bounded variation and an application to parameter estimation[END_REF]). In fact, they damp the oscillations and high frequencies (noise) in the homogeneous areas faster than any secondorder model [START_REF] Lysaker | Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[END_REF], and they preserve more efficiently the curvy features of the edges as well as the corners [START_REF] Shah | A common framework for curve evolution, segmentation and anisotropic diffusion[END_REF][START_REF] Blake | Visual Reconstruction[END_REF]. Several models are proposed in the literature and most of them resort to a strong prior on how to interpolate the solutions with high order operators. One of the very advanced approaches is due to Masnou and Morel who adapted G. Kaniza principle to the interpolation of the level lines of BV-images [START_REF] Masnou | Level-lines based disocculusion[END_REF]. It consists of minimizing an energy which combines a length term and a given power of the curvature (e.g Willmore energy, . . . , [START_REF] Masnou | Level-lines based disocculusion[END_REF][START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF][START_REF] Brito-Loeza | Fast numerical algorithms for Euler's elastica inpainting model[END_REF][START_REF] Mumford | Elastica and computer vision[END_REF][START_REF] Shah | Elastica with hinges[END_REF]).

Our contribution: In a previous work [START_REF] Theljani | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF], we considered a multiscale fourth-order model based on an adaptive method which uses a weighted bilaplacian operator; the weight being a spatially varying diffusion function. The model is updated dynamically, by choosing the weight to fit the geometry of the computed solution and by allowing a tight location of the singularity set. We have shown that this approach is based on the construction of a family of discrete energies which Γ-converges to a Mumford Shah-H -1 limit. The size of the singularity set is bounded by the Lebesgue measure of the set where the diffusion is minimal. The curvature feature is preserved thanks to an H -1 -filtering. In this article, we improve this approach in two directions: we add a length term to the functional (mimicking the Euler elastica functional [START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF][START_REF] Brito-Loeza | Fast numerical algorithms for Euler's elastica inpainting model[END_REF][START_REF] Mumford | Elastica and computer vision[END_REF][START_REF] Shah | Elastica with hinges[END_REF]) and we extend the adaptive strategy to two diffusion functions, which allows us to balance more efficiently both shortness and the curvature of the level lines in the inpainted region. Moreover, the family of the discrete energies associated to this approach converges in the Γ-convergence sense to a generalized Mumford-Shah-H -s , (0 < s < 1) functional. We emphasize that the proposed model allows us to capture the singular set of the solution and, as byproduct, the enhancement of the filtering of the data term with respect to a well-suited Sobolev norm. Besides, we obtain two gradient flows associated to the problem leading to PDEs systems of second-order that we solve with mixed formulations finite element methods.

The article is organized as follows: In Section 2, we prove the existence of H 1 -solution of the evolution equation (2.1) and its stationary counterpart. In Section 3, we give a simple and reliable discretization based on mixed finite elements. Then we present the adaptive strategy for the choice of α and β and the limiting behavior of the proposed iterative algorithm. Finally, in Section 4, we present several numerical examples which illustrate the efficiency of our model in both geometry inpainting and restoration problems.

Variational formulation in H

1 (Ω) -H -1 (Ω).
2.1. Model problem. We define the image as a mapping from a bounded open set Ω of R 2 , with piecewise smooth boundary ∂Ω, into R. We consider the following equation:

(2.1)      ∂ t u + a∆ β (∆ α u) -b∆ β u + λ D (u -f ) = 0 in R + × Ω, ∂u ∂n = ∂∆αu ∂n = 0, on R + ×∂Ω, u(0, x) = f, in Ω,
where a, b > 0 are two constants, ∆ β u = div(β(x)∇u). We assume that we are given two constants α 0 , α 1 such that:

0 < α 0 ≤ α(x), β(x) ≤ α 1 , ∀x ∈ Ω.
We denote by D ⊂⊂ Ω the missing area and we set λ D = λ 0 1 in the image restoration problem (the damaged region is then empty), and λ D = λ 0 χ Ω\D where χ Ω\D is the characteristic function of the domain Ω \ D.

We emphasize that the choice of the differential operator in (2.1) may appear somewhat surprising as the operator ∆ β and ∆ α do not commute and the operator, without further assumptions on α and β, is not self-adjoint which might results in "unwished" difficulties for solving the problem. However, the main goal of the proposed approach is to capture the singular set of the solution, with preserving long range curvatures, which means that α and β should coincide (say equal 1) except on a set of zero Lebesgue measure (characterized by either α or β, or both, is equal to 0). We do not make this simplifying assumption because the singular set is unknown a priori and our numerical method gives approximated singular sets which have some non nul Lebesgue measure (even small). In addition, as we use an H 1 -H -1 formulation of the problem we circumvent this apparent difficulty.

Note that the constants a, b could be incorporated in the definition of α and β, however we choose to express them as in (2.1) in order to single out the influence of the second-and the fourth-order terms in the model. Following [START_REF] Theljani | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF], the selection of the functions of α and β is performed locally and adaptively -using residual error indicators and an efficient mesh adaptation -, to decrease, respectively increase, their values close to, respectively far from, the edges and corners. The numerical algorithm which realizes this program is very simple: it consists in solving linear problems (with fixed α and β) and to update the diffusion functions with explicit formulae. This gives a reliable numerical method, ready for the use of highperformance computing techniques of PDEs [START_REF] Belhachmi | A massively parallel multi-level approach to a domain decomposition method for the optical flow estimation with varying illumination[END_REF], low cost (the final meshes are coarsened in the homogeneous areas -usually the largest part of low textured images -) and converging to minimizers of nonlinear functionals (including the Hausdorff measure of the singular set).

The formulation

H 1 -H -1 . Let ∆ -1
β be the inverse of the negative -β-weighted Neumann-Laplacian,∆ β , with homogeneous boundary conditions, i.e., for g ∈ L 2 0 (Ω) = h ∈ L 2 (Ω); Ω h dx = 0 we set u = ∆ -1 β g to be the unique solution in V = L 2 0 (Ω) ∩ H 1 (Ω) of the following Neumann problem:

(2.2) -∆ β u = g, in Ω, ∂u ∂n = 0, on ∂Ω.
In the space L 2 0 (Ω), we define the inner product:

u, v -1,β = β 1 2 ∇∆ -1 β u, β 1 2 ∇∆ -1 β v 2 ,
and associated norm u -1,β = β

1 2 ∇∆ -1 β u 2 .
Throughout the paper, we assume that the domain Ω is partitioned into I disjoint sub-domains (Ω ) such that α and β are given by the piecewise-constant scalar functions:

α = α and β = β , in Ω , = 1, . . . , I.
We denote α m = min

1≤ ≤I α > 0, α M = max 1≤ ≤I α , β m = min 1≤ ≤I β > 0 and β M = max 1≤ ≤I β . Under
the above assumptions, we can prove the existence of the solution for (2.1) which also involves the study of the stationary boundary value problem associated to (2.1). More precisely, we prove the existence of the weak solution for the following stationary equation:

(2.3) a∆ β (∆ α u) -b∆ β u + λ D (u -f ) = 0, in Ω, ∂u ∂n = ∂∆αu ∂n = 0, on ∂Ω.
Up to an additive constant and without loss of generality, we will assume that χ D f ∈ L 2 0 (Ω). Notice that in the inpainting problem the image intensity f is given only in Ω \ D. Let v be a function in L 2 (D) extended by 0 in Ω and such (2.4)

F (v) ∈ L 2 0 (Ω) where F (v) = f in Ω \ D v in D,
and consider the alternative problem

(2.5) -a∆ α u + bu + λ 0 ∆ -1 β (u -F (v)) = 0, in Ω, ∂u ∂n = ∂(∆ -1 β (u-F (v))) ∂n = 0, on ∂Ω.
that we call H 1 -H -1 formulation. Now, consider T : L 2 (Ω) -→ L 2 (Ω) such that, for a given v like in (2.4), T (v) = u the unique solution of (2.5) (we prove below that T is well defined). Thus, if a fixed point u = v of T exists, it will be a solution of the equation (2.3). In fact, such a fixed point verifies:

(2.6)

-∆ α u + bu + ∆ -1 β (λ D (u -f )) = 0, in Ω, ∂u ∂n = ∂(∆ -1 β (λ D (u-f )) ∂n = 0, on ∂Ω.
By applying the operator -∆ β to problem (2.6), we get that the fixed point, equivalently, is a weak solution to the system (2.3), see Remark 2.3. We first prove that for a fixed v, the system (2.5) admits a unique solution u. To do so, we introduce the energy functional:

(2.7) J (u, v) = a Ω α(x) 2 |∇u| 2 dx + b 2 u 2 2 + λ 0 2 u -F (v) 2 -1,β for (u, v) ∈ V × L 2 (Ω)
Notice that a minimizer of u -→ J (u, •) satisfies system (2.5).

Proposition 2.1. There exists R > 0, such that for v ∈ L 2 0 (Ω) ∩ B -1 (0, R), the functional J (•, v) admits a unique minimizer u ∈ V . Moreover, u ∈ B -1 (0, R) Proof. We start by proving the existence of R > 0. Multiplying the first equation of (2.5) by u and integrating over Ω, it follows

a Ω α|∇u| 2 + b Ω u 2 + λ 0 Ω u ∆ -1 β u dx = λ 0 Ω u ∆ -1 β F (v) dx, which means Ω u ∆ -1 β u dx ≤ Ω u ∆ -1 β F (v) dx.
Thus, using the definition of ∆ -1 β • and the boundary conditions on u, we have:

Ω β∇∆ -1 β u • ∇∆ -1 β u dx ≤ Ω β∇∆ -1 β F (v) • ∇∆ -1 β u dx.
Using the Cauchy-Schwarz inequality for the scalar product < ., . > -1,β , we obtain:

(2.8) u 2 -1,β ≤ F (v) -1,β u -1,β .
Thanks to the continuity of ∆ -1 β , we have:

(2.9) < F (v), F (v) > 1 2 -1,β ≤ C F (v) L 2 (Ω) ≤ 2C max( f L 2 (Ω) , v L 2 (Ω) ).
Let R 1 > 0 be a given constant and choose v satisfying (2.4) and such that

v 2 ≤ R 1 . We have by direct computation v -1,β ≤ α 1 R 1 Choosing R ≥ 2C max( f L 2 (Ω) , α 1 R 1 ), we deduce that for v ∈ L 2 0 (Ω) ∩ B(0, R 1 ) (the closed ball of radius R 1 with respect to the L 2 -norm), v belongs to B -1 (0, R) (the closed ball with radius R with respect to the norm • -1,β ). Besides, the inequality (2.8) implies that u ∈ B -1 (0, R), i.e, (2.10) u -1,β ≤ R.
For the existence and uniqueness of the solution, the proof can be handled by standard techniques of calculus of variations. In fact, it is readily checked by direct calculation that the functional J (•, v) is strictly convex, weakly lower semi-continuous (as a sum of norms and thanks to the continuity of ∆ -1 β ), and V -coercive (uniformly with respect to v). Thus, J (•, v) admits a unique minimizer u in the closed convex V ∩ B -1 (0, R). Now, we prove the existence of a fixed-point of T using Schauder fixed-point theorem.

Proposition 2.2. The operator T admits a fixed point u ∈ V ∩ B -1 (0, R). Moreover, u is a weak solution of the formulation H 1 -H -1 (2.6).

Proof. We prove here the existence of a fixed point of

T . Let v ∈ L 2 0 (Ω) ∩ B -1 (0, R). From Proposition 2.1, the energy (2.7) admits a unique minimizer u = T (v) in V ∩ B -1 (0, R). Since the embedding V → L 2 (Ω) is compact for the norm • -1,β (
as composition of a bounded and a compact operators), the operator T then maps L 2 (Ω) → K, where K is a compact subset of L 2 (Ω) for the norm • -1,β . Thus we have:

T : B -1 (0, R) -→ B -1 (0, R) ∩ K = K,
where K is a compact and convex subset of L 2 (Ω) for the norm • -1,β . To apply Schauder's fixed-point theorem, it remains to prove that T is continuous for the norm • -1,β . Let (v k ) k≥0 be a sequence which converges to v ∈ B -1 (0, R) and T (v k ) = u k . The function u k is then the unique minimizer of (2.7) associated with v k , and we have:

J (u k , v k ) ≤ J (0, v k ), i.e., J (u k , v k ) ≤ λ 0 2 F (v k ) 2 -1,β ,
from the boundedness of β and since v k -1,β ≤ R and f -1,β ≤ C for some given constant C > 0, we get the following estimate:

J (u k , v k ) ≤ C, C > 0,
and then (u k ) k≥0 is uniformly bounded in V . Thus, we can consider a convergent subsequence

u k j j→∞ u * ∈ V and u k j -→ j→∞ u * in L 2 (Ω). Hence, the unique (weak) solution T (v k ) = u k of: (2.11) -a∆ α u k + bu k + λ 0 ∆ -1 β (u k -F (v k )) = 0, in Ω, ∂u k ∂n = ∂(∆ -1 β λ 0 ∆ -1 β (u k -F (v k )) ∂n = 0, on ∂Ω,
weakly converges to the unique weak solution of (2.5) by passing to the limit in (2.11)

J(u * , v) ≤ lim inf k J(u k , v k ) ≤ lim inf k J(w, v k ) ≤ J(w, v) ∀w ∈ V.
that is u * = u. From the uniqueness of the minimizer of (2.7), we obtain u = T (v). We then deduce that T is continuous in L 2 (Ω) and the existence of a stationary solution u follows from Schauder's fixed-point theorem.

Remark 2.3. i) The formulation of problem (2.6) and the regularity theory for elliptic equations show that since α and β are discontinuous, the couple (u, ∆ -1

β (λ D (f -u))) belongs to the space [H s (Ω)] 2 , 1 ≤ s < 3
2 (e.g. see [START_REF] Belhachmi | Weighted harmonic and complex Ginzburg-Landau equations for gray value image inpainting[END_REF]) so that even if the operator ∆ β (∆ α •) is of fourth-order, no global regularizing effect can be expected for u. However, its hypo-ellipticity gives that u is smooth in the homogeneous areas where α and β are constants (u is at least

H 3 (Ω) if f ∈ V ).
Thus oscillations are damped strongly in these areas while the singular set is preserved.

ii) It may be checked that problem (2.6) and (2.1) are equivalent in the following sense: the unique solution u of (2.6) is the weak solution of (2.1) written as a uniformly elliptic system on the variables (u,

w = ∆ -1 β (λ D (f -u))).
iii) It is interesting to notice that there are two ways to transform (2.1) into a PDEs system of second order of two variables (u, w), precisely by choosing

-∆ α u = w or -∆ β w = λ(u -f ).
The first choice is more usual in the computer sciences community, as the operator is the usual bilaplacian and thus is obtained as the gradient flow of the H 2 -energy. Beyond some differences in the analysis of this two choices, we show that in the discrete setting as well as with regard to the numerical implementation the two systems are close and yields essentianlly the same results.

2.3. The evolution equation as a gradient flow. The previous existence of a stationary solution applies to show that the following unbounded operator (in L 2 ):

A(u) = a∆ β (∆ α u) -b∆ β u + λ D u,
is maximal. Moreover, we have:

A(u), u 2 ≥ a I =1 β m α m ∆u, ∆u L 2 (Ω ) + b I =1 β m ∇u, ∇u L 2 (Ω ) + λ D u, u L 2 (Ω) ≥ 0,
which means that it is monotone. Thus it follows from the theory of maximal monotone operators [START_REF] Brézis | Opérateurs Maximaux Monotones et Semi-Groupes de Contractions dans les Espaces de Hilbert[END_REF], that the evolution problem (2.1) admits a unique solution u ∈ L 2 ((0, +∞);

V )∩ H 1 ((0, +∞); V ) ∩ C([0, +∞[ ; L 2 (Ω)).
2.3.1. First gradient flow. Since the product (∆ β u, ∆ α v) is not symmetric, the equation (2.1) is not a gradient flow of an energy functional. However, based on the properties of the operator A, we may write the problem (2.1) as "a formal gradient flow" of an "H 2energy", that we can split into two second-order equations by introducing an auxiliary function w ∈ L 2 (R + ; V ) such that:

(2.12)

           ∂ t u -a∆ β w -b∆ β u + λ D (u -f ) = 0, in R + × Ω, -∆ α u = w, in R + × Ω, ∂u ∂n = ∂w ∂n = 0, on R + × ∂Ω, u(0, x) = f (x), in Ω.
Therefore, a weak solution of the previous system is defined as a pair (u, w)

∈ L 2 (R + ; V ) × L 2 (R + ; V ), u(0, x) = f (x) such that: for a.e t > 0 (2.13) ∂ t u, φ + b β∇u, ∇φ 2 + a β∇w, ∇φ 2 + λ D u, φ 2 = λ D f, v 2 , ∀φ ∈ H 1 (Ω), α∇u, ∇ψ 2 -w, ψ 2 = 0, ∀ψ ∈ H 1 (Ω).
Without loss of generality we have replaced the inner product •, • V ,V with the L 2 -inner product. It is easy to check that the pair (u, w = -∆ α u) is a weak solution of (2.13) where u is the solution for the problem (2.1).

To prove uniqueness, we consider another solution (u 1 , w 1 ) ∈ L 2 (R + ; V ) × L 2 (R + ; V ) to the system (2.13), we then have: for a.e. t > 0

     ∂ t (u -u 1 ), φ 2 + b β∇(u -u 1 ), ∇φ 2 +a β∇(w -w 1 ), ∇φ 2 + λ D (u -u 1 ), φ 2 = 0, ∀φ ∈ H 1 (Ω), α∇(u -u 1 ), ∇ψ 2 -(w -w 1 ), ψ 2 = 0, ∀ψ ∈ H 1 (Ω).
Let (ζ ) I =1 be a partition of the unity associated to the decomposition (Ω ) , and by picking ψ = β α -1 ζ (w -w 1 ), in the second equation, we have the identity: For each ,

(2.14) Ω α -1 β ζ (w -w 1 ) 2 dx = Ω β α ∇(u -u 1 ) α -1 ∇ζ (w -w 1 ) dx.
Summing up, we get (2.15) a

I =1 α -1 β Ω ζ (w -w 1 ) 2 dx = a Ω β∇(u -u 1 ) ∇(w -w 1 ) dx ≥ 0.
Choosing the test function φ = u -u 1 in the first equation and using (2.15) and the positivity of α and β, we obtain:

∂ t (u -u 1 ), u -u 1 2 = -b β∇(u -u 1 ), ∇(u -u 1 ) 2 -λ D (u -u 1 ), (u -u 1 ) 2 -a β∇(u -u 1 ), ∇(w -w 1 ) 2 ≤ 0. (2.16) So that ∂ t (||u(t) -u 1 (t)|| 2 2 ) = 2 ∂ t (u -u 1 ), u -u 1 2 ≤ 0.
It follows that the function t → ||u(t) -u 1 (t)|| 2 2 is decreasing on R + . Since u(0) = u 1 (0), we get u = u 1 , in R + × Ω. It follows that w 1 and w are solutions of the same variational equation: for a.e. t > 0, r = w, respectively r = w 1 and

∂ t u, φ + b β∇u, ∇φ 2 + a β∇r, ∇φ 2 + λ D u, φ 2 = λ D f, φ 2 , ∀φ ∈ H 1 (Ω),
which admits a unique solution. Thus w = w 1 in R + × Ω and the problems (2.1), (2.12) are equivalent. In addition, from the elliptic regularity theory, we have that w ∈ L 2 (R + ; V ), thus for the discretization we will seek (u, w)

∈ L 2 (R + ; V ) × L 2 (R + ; V ).

Second gradient flow.

We introduce now a second gradient flow, namely the formulation H 1 -H -1 , which is the gradient flow associated to the energy functional J (u, u) in (2.7) and that we write as a second-order system : find (u, w) ∈ L 2 (R + ; V ) 2 such that:

(2.17)

           ∂ t u -a∆ α u + bu + w = 0, in R + × Ω, -∆ β w + λ D (u -f ) = 0, in R + × Ω, ∂u ∂n = ∂w ∂n = 0, on R + × ∂Ω, u(0, x) = f (x), in Ω.
Therefore, a weak solution of problem (2.17) is defined as a pair (u, w) ∈ L 2 (R + ; V ) 2 , u(0, x) = f (x) such that: for a.e t > 0

(2.18)

∂ t u, φ + a α∇u, ∇φ 2 + u, φ 2 + w, φ 2 = 0, ∀φ ∈ H 1 (Ω), β∇w, ∇ψ 2 -λ D u, ψ 2 = λ D f, ψ 2 , ∀ψ ∈ H 1 (Ω). Notice that w, φ 2 = ∆ -1 β (λ(u -f ), φ 2 = λ D (u -f ), φ -1,β
, therefore we obtain the uniqueness of the solution of problem (2.17) by choosing φ = u 1 -u 2 in the first line of the system and arguing as previously to deduce the desired result.

The discretization of this two PDEs system, based on mixed formulations, leads to similar discrete problems (with obvious changes), thus we only describe it for the first system.

3. Full discretization and adaptive strategy. We assume that the domain Ω is polygonal. We consider a regular family of triangulations T h made of elements which are triangles (or quadrilaterals) with a maximum size h, satisfying the usual admissibility assumptions, i.e., the intersection of two different elements is either empty, a vertex, or a whole edge. For h > 0, we introduce the following discrete space:

X h = v h ∈ C(Ω)|∀K ∈ T h , v h | K ∈ P 1 (K) ∩ V,
where P 1 (K) is the set of a polynomial function on K of total degree ≤ 1 and K is a triangle in T h .

For ∆t > 0, the discretized version of the problem (2.13) leads to finding a pair (u

h , w h ) ∈ X h × X h solution of (3.1)        u h -u old h ∆t , φ h 2 + b β∇u h , ∇φ h 2 +a β∇w h , ∇φ h 2 + λ D u h , φ h 2 = λ D f h , φ h 2 , ∀φ h ∈ X h , α∇u h , ∇ψ h 2 -w h , ψ h 2 = 0, ∀ψ h ∈ X h ,
where f h is a finite-element approximation of f associated with T h . Since X h ⊂ V , the well-posedness of the problem (3.1) follows as in the continuous case. Before we explicit the adaptive strategy, we present a sketch of the Γ-convergence analysis of such an approximation approach which shows how our diffusion model converges to the Mumford-Shah-H -s functional and therefore performs the inpainting by preserving the jump set and the curvature.

Γ-convergence approximations.

The Γ-convergence study of the adaptive strategy for a functional of second order was presented in [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF] for the optic flow estimation. The authors proved that this algorithm is equivalent to the adaptive one introduced by Chambolle-Dal Maso [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] and Chambolle-Bourdin [START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF] where a numerical discrete approximation of the Mumford-Shah energy to the segmentation problem was proposed. This approach, based on finite-element discretization and adaptive mesh strategy, turns out to be a well suited approximation with a family of discrete energies, in the Γ-convergence sense [START_REF] Braides | Gamma-Convergence for Beginners[END_REF] to the Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] (see [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF], [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF] for more details).

Let E denotes the energy functional associated with (2.3) (after applying ∆ -1 β to the equation)

E(u, v) = a 2 Ω α(x)|∇u| 2 dx + Φ(u, F (v)),
where

(3.2) Φ(u, F (v)) = b 2 u 2 2 + λ 0 2 u -F (v) 2 -1,β ,
which measures the data fitting in the norm • -1,β and acts as a filtering term. The first term in E(•) is a regularization part, and if the diffusion function can be "chosen" such that it equals 1 except on a given (compact) set, of null Lebesgue measure, where it equals zero then this part of the functional will constitutes an approximation, in the Γ-convergence sense, of the energy of such a set. In other words, in this case, this term will converges to H 1 (S u ).

To show this, we follow the proof in [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] which consists in introducing a family of discrete energies, that stem from a particular discretization of the regularizer. We briefly recall the method. For a fixed angle 0 < θ 0 ≤ 2π/3, a constant c ≥ 6, and for > 0, let T (Ω) = T (Ω; θ 0 ; c) be the set of all triangulations of Ω whose triangles K have the following characteristics:

(i) The length of each of the three edges of K is between and c. (ii) The three angles of K are greater than or equal to θ 0 . Let V (Ω) be the set of all continuous functions u : Ω -→ R such that u is affine on each triangle K of a triangulation T ∈ T (Ω). For a given u, we define T (u) ⊂ T (Ω) as the set of all triangulations adapted to the function u, i.e., such that u is piecewise affine on T. We consider a non-decreasing continuous function g : [0, +∞) -→ [0, +∞) such that:

lim t→0 g(t) t = 1, lim t→+∞ g(t) = g ∞ < +∞.
For any u ∈ L p (Ω), (p ≥ 1) and T ∈ T (Ω), the authors in [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] introduced the following minimization problem:

(3.3) G (u) = min T∈T (Ω) G (u, T), where G (u, T) =    K∈T |K ∩ Ω| 1 h K g(h K |∇u| 2 ), u ∈ V (Ω), T ∈ T (u), +∞ , otherwise.
When goes to zero and provided θ 0 is less than some Θ > 0, it was proven that the energy G Γ-converges to the Mumford-Shah functional:

G(u) =    Ω |∇u(x)| 2 dx + g ∞ H 1 (S u ), u ∈ L 2 (Ω) ∩ GSBV, +∞, u ∈ L 2 (Ω)\GSBV,
where is GSBV the generalized special function of bounded variation (see [START_REF] Ambrosio | Functions of Bounded Variation and Free Discontinuity Problems[END_REF]).

Remark 3.1. If Φ : X -→ [-∞, +∞] is continuous and (G ) Γ-converges to G then (Φ + G ) Γ-converges to Φ + G.
Here, we assume that the function β(•) is fixed. Then we have the following proposition: Proposition 3.2. Let f ∈ L 2 0 (Ω) and > 0 be a positive parameter. The sequence of functionals

(G (•) + Φ(•, F (v))) Γ-converges for -→ 0 in the topology of L 2 (Ω) to G(•) + Φ(•, F (v)).
Proof. We have

G (u ) + Φ(u , F (v)) = G (u ) + Φ(u, F (v)) + Φ(u , F (v)) -Φ(u, F (v))
From the result of the Γ-convergence of G to G, see [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF]Theorem 2] and the previous remark we may handle the first two terms. Notice that, for fixed β and v, the continuity of the term Φ(u , F (v)) of the functional in L 2 (Ω) follows from the continuity of ∆ -1 β , (i.e. stability of the elliptic problems and the Poincaré-Wirtinger inequality). Therefore, we have

G (u ) + Φ(u, F (v)) Γ-converges to G(u) + Φ(u, F (v)),
and standard convergence results of the finite element solution in the energy norm give:

Φ(u , F (v)) -Φ(u, F (v)) -→ 0 when → 0
In the theorem we note that the function α does not appear explicitly in G or G . To introduce such an α, let ψ be the Legendre-Fenchel transform of g. For a given triangulation T , it was proven in [START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF] that the minimization of G is equivalent to the minimization of the following function

G (u, v, T ) = K∈T |K ∩ Ω| v K |∇u| 2 + ψ(v K ) h K , over all u ∈ V (Ω) and v = (v K ) K∈T (Ω)
, piecewise constant on each K ∈ T . The new minimization problem over u and α is implemented by an alternate minimization scheme and we have: For a given u, the minimizer over each α is explicitly given by:

(3.4) α K = g (h K |∇u| 2 ).
The minimizer u satisfies

(3.5) div(α K ∇u) + ∇Φ(u, F (v)) = 0, in Ω,
where the ∇Φ(u, F (v)) denotes the Fréchet derivative of Φ(u, f ) over L 2 (Ω). Then, from equation (3.5), we can see now how the function α K appears in the minimization of the energy G and the link with our problem. The choice of α, as explained above, depends on the function g. An example of such a function is given by [START_REF] Chambolle | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF]

g(t) = 2g ∞ π arctan ( π t 2g ∞ )
which gives a smoothly decreasing sequence α (as the inverse of the gradient square). In what follows, we will use another way to decrease, in a geometric way, such an α. Moreover, we give a unified framework to perform the minimization in (3.3) and the choice of α based on appropriate a posteriori error indicators. We emphasize that the function β is not fixed in our approach and is also chosen adaptively to control the diffusion in the length term. Thus, writing

E(u, v) = a 2 Ω α(x)|∇u| 2 dx + b 2 u 2 2 + λ 0 2 Ω β(x)|∇w| 2 dx, where (3.6) -∆ β w = u -F (v), in Ω, ∂w ∂n = 0, on ∂Ω,
and

Ω β|∇w| 2 dx = Ω (F (v) -u) w dx.
Therefore, we may write equivalently

E(u, w; v) = a 2 Ω α(x)|∇u| 2 dx + b 2 u 2 2 + λ 0 2 Ω β(x)|∇w| 2 dx - λ 0 2 Ω (u -F (v)) w dx,
and for > 0, we define E (u, w; v) = E(u, w; v) + w 2 2 . Similarly to the discrete energy G(u ) in (3.3), we define the energy H (w ) which approximates, in the Γ-convergence sense, the term λ 0 2 Ω β(x)|∇w| 2 dx. Then, we have

H (w ) Γ-converges to H(w) = λ 0 2 Ω |∇w| 2 dx + λ 0 h ∞ H 1 (S w )
Since for any > 0, w solves (3.6) (with v and u in the right-hand side), we have from the finite element theory

λ 0 2 Ω (u -F (v))w dx converges to λ 0 2 Ω (u -F (v))wdx
Therefore, setting

Ψ(u, w; v) = λ 0 2 Ω λ D (u -F (v)) w dx,
we write

E(u , w ; v) = a 2 Ω α(x)|∇u | 2 dx + b 2 u 2 2 + λ 0 2 Ω β(x)|∇w | 2 dx+ Ψ(u, w; v) + (Ψ(u , w ; v) -Ψ(u, w; v)).
Arguing like previously we approximate E (u , w ; v) with G (u ) + H (w ) + Ψ(u, w; v) + (Ψ(u , w ; v) -Ψ(u, w; v)) and we obtain Proposition 3.3. Let f ∈ L 2 0 (Ω) and > 0 be a positive parameter. The sequence of functionals

(G (•) + H (•) + Ψ(•, •; v)) , Γ-converges for -→ 0 to (G(•) + H(•) + Ψ(•, •; v)) in L 2 (Ω)-topology.
3.2. The adaptive approach. The adaptive strategy is based on two ingredients [START_REF] Belhachmi | Control of the effects of regularization on variational optic flow computations[END_REF]: the usual mesh adaptation and a "functional" adaptation which consists in choosing locally the diffusion coefficients in order to "cut" high gradients of the computed solution. In fact, mesh adaptation in finite-elements method has become an efficient tool for the computation of a numerical solution of partial differential equations, in particular it allows us to localize the singular set of the solution. Such localization leads to a coarsening of the mesh in the homogeneous areas (regions where the error indicator is small) and to refine near the singularities. Besides near the singular set, the diffusion coefficients tends to smooth the solution so that decreasing their values in these regions allows us to preserve such features.

In image processing, it is well-known that the discontinuities (edges) are contained in regions where the brightness changes sharply and consequently where the gradient of the image is large. Then, for each element K ∈ T h , we use the following local quantities:

(3.7) η K = α 1 2 K ∇u h L 2 (K) and η K = β 1 2 K ∇w h L 2 (K) ,
as error indicators.

Remark 3.4. Rigorously, the error indicator should be a mean value of the projection of ∇ u h in ω K = {K ; |K ∩ K | = 0} that is all the elements sharing an edge with K [START_REF] Verfurth | A Review of A Posteriori Error Estimation and Adaptive Mesh-Refnement Techniques[END_REF]. This allows also to reduce the artifact in the refinement process. Nevertheless, our rough indicator η K is faster to compute and already gives satisfactory results, so that is why we use it to.

The error indicator η K plays the same role as η K but for the variable w so it captures the brightness change in the Laplacian ∆ α u which gives a supplementary information on the discontinuities of the "second kind". Thus, the two quantities in (3.7) act as an edge detector (see also [START_REF] Auroux | Application of the topological gradient method to color image restoration[END_REF][START_REF] Auroux | Contour detection and completion for inpainting and segmentation based on topological gradient and fast marching algorithms[END_REF][START_REF] Auroux | Image processing by topological asymptotic expansion[END_REF]) which makes it well suited to control and locally select the diffusion coefficient α, β using the following algorithm (see [START_REF] Belhachmi | Control of the effects of regularization on variational optic flow computations[END_REF]): Algorithm 3.1 Adaptivity steps 1. Choose initial guesses α 0 , β 0 and compute u 0,h solution of (2.1) with the initial grid T 0 h corresponding to the image.

Adaptive steps: For

k ≥ 1 • Compute u k,h on T k h with α = α k and β = β k fixed.
• Build an adapted mesh T k+1 h (in the sense of the finite element method, i.e., with respect to the parameter h) with the metric error indicator.

• Perform a local choice of α(x) and β(x) on T k+1 h to obtain new functions α k+1 and β k+1 .

If ||u

k+1,h -u k,h || ≤ or k ≥ k max , stop. Otherwise k = k + 1, go to Step 2.
During the adaptation, we use the following formula for each triangle K;

α k+1 K = max      α k K 1 + κ η K η ∞ -0.1 + , α thr     
, where α thr is a threshold parameter and κ is a coefficient chosen to control the rate of decrease in α, u + = max(u, 0). The same formula is used to update the function β(x) and where η is replaced by η . The mesh refinement is generated as follows: we cut the element K close to the edges of the image (i.e., where the error (η K ) K∈T h is larger than the mean value) into a finite number of smaller elements to decrease the error and to fit the geometry of the edges, while, the mesh is coarsened in the other regions. This process produces a tight localization of the singular set (see the next section). Then, we make an "optimal" choice of the functions α and β, following the error map furnished by the error indicators (η K ) K∈T h , in each element K. We decrease the values of the diffusion parameters in order to correctly approximate the edges. The values may be chosen as regular functions of the computed gradient (e.g.

1 1+|∇u h | 2 )
or in order to give a geometric sequence as in the previous algorithm. In all cases, the size of the singular set (its Hausdorff measure) is controlled by | {x ∈ Ω; α(x) = α m } |, a and b (where | • | denotes the Lebesgue measure). Note that the function β influences also this size (in a less explicit way) but plays a more important role in the filtering part of the data term. The overall result expresses a balance between the filtering and the segmentation effects which is the key point to obtain curvy and short reconstructed paths.

Numerical examples.

In this work, all the PDEs are solved with the open source finiteelement software FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF]. The goal is to reconstruct the missing information in the red parts, i.e., D, by the diffusion of the information from the intact part, i.e, Ω\D. In all examples, we set ∆t = 0.1, an initial guess u(0, x) = 0.5, ∀x ∈ D, λ 0 = 10 4 in the case of image restoration and λ 0 = 10 7 in the case of image inpainting.

4.1. Image restoration. We begin by testing our approach for an image denoising problem. In the left-hand plot of Fig. 4.1, we display the original (binary) image with a squared domain (120 × 120 pixels). In the middle, we display the noisy image f , obtained by adding a Gaussian noise, whereas the right-hand plot shows the restored one obtained by using our approach (2.1).

We initialized the algorithm with a large value of α = β = 50 and we performed 20 iterations of the adaptive algorithm. We plot in Fig. 4.2 the mesh, the function α and the error indicator at convergence. The latter indicates the regions where we have edges whereas α K plays in the present context a role similar to that of the z-field in the Ambrosio-Tortorelli approximation method for the Mumford-Shah energy. We can also see the "sparsification" effect on the mesh in the left-hand plot of Fig. 4.2, which emphasizes the low cost of the method.

We tested our model on a gray-scale image (200×200 pixels) in Fig. 4.3 where we displayed the original (noisy) image and the restored one. In Fig. 4.4, we vary the level of noise in the input images. In the top row, we superpose the three noisy images corrupted by Gaussian noise having zero mean and standard deviation 25 In the bottom row, we show the restored images. Large gap inpainting: We test the proposed model in interpolating large gaps along the inpainting domain in response to the values of the two weights a and b in (2.1). More precisely, we want to see if our approach fulfills the connectivity principle. In Fig. 4.5, we give the restored images using total variation [START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF], harmonic [START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF] and biharmonic models (our model for α = β = 1), respectively. The harmonic inpainting does not achieve any connectedness and produces a smooth solution u in D, blurring the edges. TV inpainting is unable to connect and preserve the edge of the strip. In the same figure, we display the restored images using the model (2.1) for different values of a and b. In contrast to harmonic, TV and biharmonic models, the proposed approach is able to connect large gaps and edges are well-captured. The results show that we have to favor the fourth-order term by choosing b small, because, matching edges across large distance is the main interest of higher-order PDEs. We notice here that this example illustrates the main difference between this combined model and the pure fourth-or second-order models proposed in [START_REF] Belhachmi | Weighted harmonic and complex Ginzburg-Landau equations for gray value image inpainting[END_REF][START_REF] Theljani | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF]. In fact, considering only fourthorder model in [START_REF] Theljani | A multiscale fourth-order model for the image inpainting and low-dimensional sets recovery[END_REF] allowed for reconstructing curved edges, and this can easy seen in the last image where the broken edges are connected but they suffer form a curvature which is not present in the original image. Whereas considering only the second-order model [START_REF] Belhachmi | Weighted harmonic and complex Ginzburg-Landau equations for gray value image inpainting[END_REF] allows for a similar result to the total variation, which means that the edges will be disconnected.

Curvature inpainting: In Fig. 4.6, we tested the adaptive algorithm for different values of the ratio b a in equation (2.1), in order to show the effects of each term (the fourth-and second-order one). We plot in Fig. 4.6 the restored images at the end of algorithm. It appears that if more weight is set on the second-order derivatives in (2.1), then the inpainted edge tends to be a straight line as expected (the length term is enforced). In fact, if we consider equation (2.1) without the fourth-order derivatives, the adaptive algorithm gives a solution u that converges to the one given by the Mumford-Shah model in the sense of the Γ-convergence [START_REF] Belhachmi | An adaptive approach for the segmentation and the TV-filtering in the optic flow estimation[END_REF][START_REF] Chambolle | Discrete approximation of the Mumford-Shah functional in dimension two[END_REF]. It is well know that the preferable edge curves in the Mumford-Shah model are those which have the shortest length because the penalization term acts on the length of the edge only. Therefore, promoting the second term allows the model to favors straight edges to commit the connectivity principle in perception [START_REF] Chan | Mathematical models for local non-texture inpainting[END_REF]. In Fig. 4.6, we can see that the curvature is well inpainted which proves that our approach, based on fourth order linear diffusion models, allows us to obtain a very interesting result that one might expect by solving some more complex models like the Euler's elastica [START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF], which is highly nonlinear and numerically difficult to solve.

Blood vessel inpainting: The experiment in Fig. 4.7 shows the efficiency of the proposed method in the reconstruction of damaged parts of two-dimensional image representing a narrowed blood vessel. The zones to be inpainted are marked with red color. We illustrate in this example three advantages of the model. First, the ability of producing curved edges in region 1, second the ability, in region 2, of connecting edges across large distance, and third, in region 3, the ability of separating three vessels which from a T -junction. In Fig. 4.8, we display the mesh at different adaptation iterations. The mesh is coarsed in the homogeneous regions which makes the overall numerical procedure fast.

Noisy images inpainting: We illustrate in Fig. 4.9 the efficiency of the approach in the restoration of noisy and damaged image. Here, we add to the intact parts of the image a Gaussian noise, with zero mean and a standard deviation equal to 20. We used our algorithm to perform restoration and inpainting simultaneously. We can see that edges are connected (see Fig. 4.9) and the curvature is well reconstruted (see Fig. 4.10) while the known regions of the image were restored in parallel.

Conclusion.

We have investigated a combined second-and fourth-order PDE model for the image restoration and inpainting problems. We proved existence of a H 1 -solution for the evolution equation as well as its corresponding stationary problem. Then, we proposed a simple discretization based on mixed finite element and we detailed an adaptive strategy for the choice of the regularization parameters. The adaptation strategy produces good solutions that one might obtain by using more complex higher order PDEs. We performed numerical experiments by denoising images that have been corrupted by Gaussian noise.

Figure 4 . 1 .

 41 Figure 4.1. Original, noisy and restored images, respectively.

Figure 4 . 2 .

 42 Figure 4.2. The diffusion function α, the mesh and the error indicator η at convergence.

Figure 4 . 3 .

 43 Figure 4.3. The noisy and the restored images.

Figure 4 . 4 .

 44 Figure 4.4. In the top row and from left to right: Input noisy image obtained by adding Gaussian noise with zero mean and standard deviation 10, 15 and 25, respectively. In the bottom: Restored images.

Figure 4 . 5 .

 45 Figure 4.5. Connectivity across large gaps: First image: the damaged. Second line: Total variation, harmonic and biharmonic models and our model for the ration b a = 10, respectively. Third line: our model for the ratio b a = 5, 1, 0.1, and our model when b = 0 (i.e., only fourth-order derivatives).

Figure 4 . 6 .

 46 Figure 4.6. Top row: Original and damaged images. Bottom row: Model (2.1) and adaptation: the ratio

Figure 4 . 7 .

 47 Figure 4.7. From left to right: Damaged image,d restored one at iteration 1, restored images at iteration 5 and 20, respectively.

Figure 4 . 8 .

 48 Figure 4.8. From left to right: Initial mesh, the adapted one at iteration 1, adapted meshes at iteration 5 and 20, respectively.

Figure 4 . 9 .

 49 Figure 4.9. Connectivity across large gaps for noisy image.

Figure 4 . 10 .

 410 Figure 4.10. The damaged & noisy and restored images, respectively.