Zakaria Belhachmi
email: zakaria.belhachmi@uha.fr

Moez Kallel
email: moez.kallel@ipeit.rnu.tn

AND Anis Theljani
email: thaljanianis@gmail.com

Keywords: Parial Differential Equations -Inverse problems -Image inpainting -Regularization techniques -Higher-order PDEs -Cauchy problems -, respectively. 4.3 Curvature inpainting, example 2 Image inpainting -Inverse problems -Regularization procedures -Adaptive finite elements Image inpainting -Inverse problems -Regularization procedures -Mixed finite elements AMS subject classifications. 65M32 -65M50 -65M22-94A08 -65N22 -35G15-35Q68

Image inpainting or disocclusion refers to the process of restoring a damaged image with missing information and it has many applications in different fields. Different mathematical approaches were suggested to deal with this problem. In particular, partial differential diffusion equations are extensively used. The underlying idea of PDE-based approaches is to fill-in damaged regions with available information from their surroundings. This is generally accomplished by exploiting a wide range of diffusion equations for which the available information in the boundary of these regions are usually used as Dirichlet boundary conditions.

The first purpose of this Thesis is to treat the case where this information is not available in a part of the boundary of the damaged region. We formulate the inpainting problem as a nonlinear boundary inverse problem for incomplete images. Then, we give a Nash-game formulation of this Cauchy problem and we present different numerical which show the efficiency of the proposed approach as an inpainting method. Typically, inpainting is an ill-posed inverse problem for it most of PDEs approaches are obtained from minimization of regularized energies, in the context of Tikhonov regularization. The second part of the thesis is devoted to the choice of regularization parameters in second-and fourth-order energy-based models with the aim of obtaining as far as possible fine features of the initial image, e.g., (corners, edges, . . . ) in the inpainted region. We introduce a family of regularized functionals with regularization parameters to be selected locally, adaptively and in a posteriori way allowing to change locally the initial model. We also draw connections between the proposed method and the Mumford-Shah functional. An important feature of the proposed method is that the investigated PDEs are easy to discretize and the overall adaptive approach is easy to implement numerically.

challenging problems, and motivated me in my own creative thinking. Their guidance helped me in all the time of research and writing of this thesis.

Second, I would like to express my deep thanks to my thesis reviewers Professor Didier AUROUX and Professor Simon MASNOU for having accepted to be members of the jury, for reading this manuscript carefully and for their support and helpful comments. I am also extremely grateful to Professor Faker BEN BELGACEM for his participation in my jury.

I also would like to thank all the people who contributed in some way to the work described in this thesis during my years at the Mathematical and Numerical Modeling in Engineering Sciences Laboratory (LAMSIN) and the Mathematics, Information Technology and Applications Laboratory (LMIA). Special thanks to all my colleagues and friends for creating the friendly, pleasant and fertile atmosphere to work.

Finally, I would like to express my deepest gratitude to my family for the love, support, encouragement and understanding over all these years of my research.

Introduction Inpainting as na inverse problem

Inpainting as forward smelborp

Résumé

Le traitement des images numériques voit ses débuts vers les années 1920 dans la transmission de données par câble mais ne connaît de vrai succès que dans les années 1960 avec le développement des ordinateurs. Au départ, les techniques de traitement d'images sont essentiellement des méthodes de restauration et de compression d'images. Ce domaine s'est développé avec les progrès de l'informatique, pour toucher de nombreux traitements, comme la segmentation, la restauration et la désocclusion, dans des domaines aussi variés (médical, télévision, imagerie satellite, multimédia, musées). Dans cette thèse, on s'intéresse essentiellement au problème de reconstruction ou désocclusion d'images connu sous le nom "inpainting". C'est une tache très importante et classique en analyse d'images qui consiste à retrouver des parties cachées ou endommagées d'une image. A l'instar du débruitage d'images (dont elle fait partie), elle occupe une place centrale tant dans l'acquisition de données (images) fiables que dans les domaines de restauration (films, peintures,...) ou retouches (réalité virtuelle). Typiquement, il s'agit d'un problème inverse. D'une façon générale, nous considérons une image u définie sur un domaine Ω ⊂ R n (n = 2, 3) et un sous domaine D ⊂ Ω où l'information est inconnue. Le but est de retrouver cette information en se basant seulement sur les données disponibles dans les parties connues de l'image, c-à-d, Ω\D. Des diverses approches ont été proposées pour résoudre ce problème; l'approche par modélisation stochastique, l'approche par ondelettes et l'approche par les équations aux dérivées partielles (EDP). Cette dernière est très populaire et profite des progrès récents en analyse des équations aux dérivées partielles au niveau théorique ainsi qu'au niveau de la résolution numérique. L'approche par EDP se base sur la diffusion des données fournies sur la frontière de D vers son intérieur en exploitant des diverses types d'équations de diffusion qui s'écrivent généralement de deux façons différentes. La première consiste à poser le problème localement sur le sous-domaine D comme suit : R(u) = 0, dans D, u = f, sur ∂D,

où R représente l'opérateur différentiel de diffusion. Cette approche s'avère convaincante si les informations présentes sur ∂D sont exactes. Néanmoins, elle n'est pas toujours la meilleure façon pour aborder le problème de désocclusion. En effet, outre qu'elle est endommagée, l'image peut être aussi bruitée dans plusieurs cas, et en diffusant les données présentes sur ∂D, ce bruit va nécessairement affecter la zone à reconstruire. Pour cela, dans plusieurs travaux le problème se ramène sur tout le domaine Ω comme suit :

-R(u) + λ D (uf ) = 0, dans Ω,

∂ n u = 0, sur ∂Ω, (2)
où λ D = λ 0 χ Ω\D et n est le vecteur unitaire et normal à ∂Ω. Le paramètre λ 0 prend généralement une grande valeur de telle sorte que la solution de (2) converge vers la solution de (1).

Dans le premier chapitre de la thèse, on s'intéresse à un cas moins étudié par les méthodes classiques, où le domaine de désocclusion D touche le bord de l'image ∂Ω. Un problème de Cauchy linéaire et elliptique est alors introduit par les equations suivantes :

       ∇•[k(|∇u| 2 )∇u] = 0, dans D, u = f, sur Γ c , k(|∇u| 2 )∂ n u = φ, sur Γ c , (3)
Dans ce cas, les données de type Dirichlet et Neumann sont seulement accessibles sur Γ c . Ceci est le cadre naturel d'un problème d'extrapolation de données, contrairement aux cas classiques où la désocclusion est traitée comme une interpolation spatiale de données. De ce point de vue, une modélisation par un problème de Cauchy qui consiste à déterminer une solution d'une équation aux dérivées partielles dans D avec des données aux bords accessibles seulement sur la partie Γ c , est la façon la plus naturelle pour retrouver l'information de l'image initiale dans D. La fonction k(•), qui représente le coefficient de diffusion, vérifie certaines conditions afin de garantir l'unicité du problème du Cauchy et aussi la solvabilité d'un problème direct associé à (3).

En général, un problème de Cauchy est connu d'être mal posé au sens d'Hadamard. L'existence d'une solution n'est pas toujours assurée pour des données quelconques qui résultent de mesures et des observations. En pratique, il y a toujours une erreur de mesure commise qui peuvent entacher les données fournies. Ceci rend la solution du problème de Cauchy, même si elle existe, instable, d'où la difficulté essentielle de ce problème. L'utilisation d'une méthode de régularisation vis à vis les petites perturbations est donc indispensable dans ce cas. Elle consiste à calculer une solution approchée du problème de Cauchy initial par la résolution d'un problème bien posé vérifiant les propriétés d'existence, d'unicité et de stabilité de la solution par rapport aux données initiales. On propose dans cette thèse une approche qui consiste à définir un problème de minimisation équivalent au problème de Cauchy (3). Ensuite, on utilise des algorithmes de minimisation, basées sur la théorie des jeux, pour se remédier au caractère mal-posé du problème du Cauchy. Cette approche est une extension de celles publiées récemment dans [START_REF] Habbal | Data completion problems solved as Nash games[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic cauchy problems[END_REF] sur la résolution du problème de complétion de données manquantes via une une formulation basée sur une stratégie de jeux de Nash. Le problème de Cauchy (3) est reformulé comme un jeux de Nash à deux joueurs qui vont résoudre en parallèle deux problèmes aux limites associés. Le premier suppose connu les données de type Dirichlet f , et utilise une condition de type Neumann préscrite sur la partie inaccessible du bord, comme une variable stratégique. Contrairement, le deuxième joueur considère les données φ sur la frontière accessible du bord et joue avec une condition de Dirichlet sur l'autre frontière. Chacun des deux joueurs minimise sa propre fonction coût sous la contrainte du problème direct qui lui est associé et en fixant l'information provenant de l'autre joueur. L'échange de l'information entre les deux joueurs se fait par un terme de couplage qui, en réalité, représente le terme de régularisation.

La deuxième partie de cette thèse concerne le contrôle des paramètres de régularisation. Dans un cadre variationnel, la régularisation d'un problème inverse exige l'ajout d'un ou plusieurs paramètres de régularisation. En général, le problème régularisé donne une meilleure approximation de la solution si les paramètres de régularisation sont choisis de manière appropriée.

Pour cela, l'élaboration de stratégies efficaces en traitement d'images, et même en d'autres problèmes inverses, a suscité d'intenses recherches et une variété d'approches ont été proposées. Le problème de sélection de paramètres de régularisation, traités souvent de façon empirique et par conséquent, ne tenant pas compte des différentes zones de l'image, est au cøeur de l'approche que nous avons développés dans cette deuxième partie. Dans notre approche, on construit une famille de problèmes d'optimisation bien posés où les paramètres sont choisis comme étant une fonction variable en espace. Ceci permet de prendre en compte les différents détails, à différents échelles dans l'image. L'apport de la méthode est de résoudre de façon satisfaisante et 'objective", le choix du paramètre de régularisation en se basant sur des indicateurs d'erreur et donc le caractère à "posteriori" de la méthode (i.e. indṕendant de la solution exacte, en génórale inconnue). En outre, elle fait appel à des techniques classiques d'adaptation de maillage, qui rendent peu coûteuses les calculs numériques. En plus, un des aspects attractif de cette méthode, en traitement d'images est la récupération et la détection de contours et de structures fines.

La première partie est consacrée au choix des paramètres de régularisation pour des EDP d'ordre deux. On considère le problème de désocclusion en minimisant l'énergie suivante :

Ω α(x)|∇u| 2 dx + Ω λ D (x)(u -f ) 2 dx. (4)
On propose une méthode adaptative pour le choix du paramètre de régularisation α. La valeur de α, qui varie selon les différentes échelles dans l'image, est choisie d'une façon local et adaptative. Le but est de créer une cartographie qui guide la diffusion pour capturer les singularités (arêtes) et avoir une image bien contrastée. On présente le cadre discret du modèle proposé et on définit des indicateurs d'erreur par résidu à posteriori qui permettent la sélection du paramètre par un algorithme adaptatif. Par la suite, on établit la relation, au sens de la Γ-convergence, entre la solution obtenue avec celle du modèle de Mumford

       ∂ t u + ∆(∆ α u) + λ D (u -f ) = 0, dans R + × Ω, u = ∆ α u = 0, sur R + × ∂Ω, u(0, x) = f, dans Ω, (5)
où ∆ α u = ∇(α(x)∇u). On effectue une étude théorique sur l'existence d'une solution pour ce modèle. Ensuite, on applique l' approche adaptative, développée dans le chapitre 2, pour le choix du coefficient de régularisation. On montre aussi que l'adaptation dans ce cas approche, au sens de la Γ-convergence, un modèle qui couple l'énergie de Mumford-Shah avec un second terme de fidélité mesurant la différence entre la solution et l'image initiale en norme H -1 . Cette approche est novatrice et les résultats numériques ont montré qu'elle permet de capturer des structures géométriques fines (arêtes et des coins)

Dans le chapitre 5, afin d'approcher l'énergie de Willmore, on introduit un nouveau modèle qui combine les dérivées partielles d'ordre quatre et celles d'ordre deux comme suit :

       ∂ t u + a∆ β (∆ α u) -b∆ β u + λ D (u -f ) = 0, dans R + × Ω, ∂ n u = ∂ n ∆ α u = 0, sur R + × ∂Ω, u(0, x) = f, dans Ω, (6)
où a, b > 0 sont deux constantes. Ce modèle peut être considéré comme une "approximation" plus simple du modèle de l'elastica d'Euler où la courbure et la longueur prises en compte respectivement par ∆ β (∆ α u) et ∆ β u.

Ensuite, on effectue l'étude théorique du modèle et on présente une discrétisation simple, basée sur la méthode des éléments finis mixtes. On introduit le cadre discret de l'algorithme adaptatif pour le choix des paramètres de régularisation dans ce modèle. Par la suite, une analyse de l'algorithme adaptatif montre que ce modèle d'ordre quatre approchent, au sens de la Γ-convergence, la solution d'un nouveau modèle qui couple l'énergie de Mumford-Shah avec un second terme de fidélité mesurant, en norme H -1 , la concordance de la solution avec l'image initiale. Les résultats numériques obtenus sont satisfaisants et ont montré la capacité de reconstruire des courbures.

"Without computers we will be stuck only proving theorems that have short proofs."

I. Kenneth (1932Kenneth ( -2013)).

Chapter 1 Introduction

Motivation

Our world is invaded by pictures (photographs, signals, films, photos...) and as it is often known "a picture is worth a thousand words". Nevertheless, these visual signals, which represent the most kind of information we deal with and interpret, usually have some imperfections during the transmission process due to captors, the atmosphere or noise, which make this information unreliable. Therefore, getting a good visual quality as well as inherent information which is easy to interpret is the main challenge for people who are interested in images. The first essential step is then to improve the visual quality of images in order to make the information more useful. This is exactly the role of imaging science which consists of operations that take as input an image and produces as output another image, a set of features or data (contours, objects, etc) or parameters related to the input one. "Image Processing" roots date back to the early 1960's with the emergence of computers and it has been developed in response to problems previously mentioned. It had challenged researchers for decades, even for simple visual tasks (finding an edge, denoising problems, etc). It is currently moving ahead and it has covered many other types of processing, such as image enhancement, image restoration, image segmentation, image matching, image compression, optic flow estimation, image recognition, etc, which are involved in varied fields closely related to our real world (satellite, medical imaging, military, robotic, Photoshop, etc). Image processing started out as a research field for computer scientists and electrical engineers and seems not to attract the interest of mathematicians. Gradually, it has become highly interdisciplinary and it has drawn much attention from the mathematical community. Nowadays, various mathematical approaches, in different frameworks, are available to process and to analyze images: Fourier transform methods, statistical and stochastic modeling [5,41,42], wavelets methods [26,37,[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], PDE methods [7,13,[START_REF] Masnou | Disocclusion a variational approach using level lines[END_REF][START_REF] Paragios | The Handbook of Mathematical Models in Computer Vision[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF].

This Ph.D. thesis mainly focuses on the use of Partial Differential Equations (PDEs) in the image inpainting problem, also called "disocclusion". It consists in recovering of missing or corrupted parts of a damaged image by an interpolation process using a great variety of mathematical approaches. Among them, Partial Differential Equations (PDEs) are widely used and they are proven to be efficient and very popular in the last two decades (see [15,29,31,39,[START_REF] Kallel | The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting[END_REF]). In fact, PDEs are well established theoretically and numerically the fact Introduction which provides useful theoretical framework and gives the possibility of using an efficient, stable and accurate algorithms to compute numerical solutions. In addition, geometric Partial Differential Equations and curvature driven flows provide powerful tools to deform or reconstruct curves, surfaces and geometric objects which represent the most crucial components of the image.

This thesis is divided into two parts. The first, which is treated in Chapter 2, is related to formulating the inpainting problem as a Cauchy problem for linear or nonlinear models. In fact, the basic idea of PDE's method for image inpainting that have been proposed in the literature is to fill-in damaged regions with available information from their surroundings by interpolating image data across all boundaries of the damaged domain via a diffusion equation. The Cauchy problem, which is-known as an ill-posed one, will be then introduced when this information is not available near a part of the boundary of the damaged region. To avoid ill-posedness, we use a game theory framework and we rephrase the Cauchy problem as a game with two Nash-game players. The second part is the subject of Chapter 3, Chapter 4 and Chapter 5 and it is dedicated to the study of an adaptive control and regularization methods for second-and higher-order PDEs in image inpainting. This strategy is based on the use of simple but effective PDEs and regularization techniques. We construct a family of simple regularized functionals and we select the regularization parameters in an adaptive way which allows tractable choice for controlling the diffusion process in order to best fit the geometric specificities of the image. The adaptation process is completely a posteriori in the sense that it uses only the computed solution of the PDE instead of any prior guess on it, which allows us to select and modify dynamically the initial model.

Before discussing the material presented in this thesis, we present some fundamental concepts about image processing which is the key word of this work.

Basic concepts in image processing

A digital image is an approximation of a real scene from the real world and it is generated through image capture devices called digitizers such as a digital camera. The digital image is described in a 2-dimensional discrete space and is defined as a finite set of points, called "Pixels" (Picture Elements) which represent its smaller components. They are spaced on a regular grid and every pixel is described by its position (i, j) on the grid and a numerical value u(i, j) which represents the intensity value. The image is called binary when the intensity of each pixel takes only the value 0 (black) or the value 1 (white). Whereas, it is called a gray-scale image if the intensities are ranging between 0 (black) and 255 (white). In the RGB system, color images are described by three channels that represent the intensities of the three primary colors "Red Green Blue". Each color channel is ranged, as gray-values, from 0 to 255. By mixing the three primary colors, any color can be obtained. For example, TVs and computer monitors use RGB to create the colorful images we see on the screen.

The mathematical analysis in the framework of PDEs is usually done in the continuous setting. A continuous representation of a digital image is defined on two dimensional (in general rectangular) domain Ω as a function u : Ω -→ [0, 255] for gray-value image. The color images are analogically defined as a vector valued function u : Ω -→ [0, 255] 3 , whereas a binary image is regarded as u : Ω -→ {0, 1}.

Introduction

Elements of an inverse problem

In a wide range of image processing tasks, the challenge is to recover a desired "original" image from an observed noisy and/or damaged one. The latter usually suffers from undesirable imperfections like noise, loss, artificial artifacts and even defaults in the acquisition devices which make the investigated equations under-determined or numerical instable. Mathematically, it leads either to an operator not one-to-one (inpainting) or to solution extremely sensitive to errors that are always present in the noisy (measured) data: small changes in the measurements f can result in large changes in the solution u. Hence, image inpainting is known as an ill-posed inverse problem. In order to explain image processing problems from inverse problem point of view, we first need to understand the basic concepts of inverse problems and the ill-or well-posedness.

"Inverse problem" is very extensive subject and its applications can be found in diverse fields (geophysics, astronomy, imaging, etc). An inverse problem can be simply defined as a map between two Banach or Hilbert functional spaces X and Y as follows:

f = T u, u ∈ X and f ∈ Y. (1.1)

Introduction

The direct (forward) problem consists in determining the output f given the input u and the transform operator T . While the inverse problem is to find the input u, assuming that both T and the output of the system f are given. A meaningful concept of well-posedeness of the latter problem was defined by Hadamard [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equations[END_REF].

Definition 1.1. Let X and Y be a normed spaces, T : X -→ Y a (linear or nonlinear) mapping. The equation T u = f is called well-posed if the following three conditions hold: i) Existence: For every f ∈ Y there is (at least one) u ∈ X such that T u = f .

ii) Uniqueness: For every f ∈ Y there is at most one u ∈ X with T u = f .

iii) Stability: The solution u depends continuously on f ; that is, for every sequence (

u n ) ⊂ X with T u n → T u (n → ∞), it follows that u n → u (n → ∞).
If (at least) one of the previous conditions fails, we say that the problem (1.1) is ill-posed in the sense of Hadamard.

In image processing problems, the ill-posedness is mostly related to the instability issue. In fact, the right-hand side f ∈ Y is never known exactly, but only up to an error δ > 0, which is responsible for instability, and problem (1.1) is transformed then to a perturbed one as follows:

f = T u + δ. (1.2)
Various digital processors can be applied to the digitized image f to generate a new digital image u. We summarize in the next the most important treatments that use the Partial Differential Equations.

Image segmentation

Image segmentation consists in partitioning an input image into its constituent parts or objects. Mumford and Shah [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] have proposed a variational framework for image segmentation by minimizing the following energy:

J (u, Γ) = γ 2 Ω\Γ |∇u| 2 dx + λ 0 2 Ω (u -f ) 2 dx + |Γ|, (1.3)
which depends on two variables of different kinds; Γ is a closed subset and u is a function belonging to the Sobolev space H 1 (Ω\Γ). The first term in (1.3) is a regularization term and acts in the homogeneous regions of the image and the second one is a data fidelity term. The solution image u is intended to be partitioned in smooth regions Ω i , separated with sharp boundaries, i.e, Γ. In practice, a difficultly in minimizing (1.3) comes from the different natures of the two unknowns. For this reason, another formulation involves identifying the curve Γ with the jump discontinuities S u . Hence, the energy (1.3) depends now only on u and it is transformed then into the following functional:

J (u) = γ 2 Ω\Γ |∇u| 2 dx + λ 0 2 Ω (u -f ) 2 dx + |S u |, (1.4)
Introduction over the space of special functions of bounded variations SBV (Ω). Several approaches have been suggested to solve this minimization problem, for instance, the elliptic approximation [2], approximation by nonlocal terms [20], approximation by finite-difference schemes [START_REF] Gobbino | Finite difference approximation of the Mumford-Shah functional[END_REF]. The level set method was greatly used in segmentation in order to directly approach the minimization of (1.4). It is an efficient algorithm particularly dedicated to the minimization of energies with contours and curves variables. The curve is merely captured as the zero levelset of a Lipschitz continuous scalar function is negative inside and positive outside. Several works in image segmentation are related to this strategy, for instance, piecewise constant level set [START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF], multiphase level set segmentation [START_REF] Kallel | A Nash-game approach to joint image restoration and segmentation[END_REF][START_REF] Lie | Piecewise constant level set methods and image segmentation[END_REF]. The main issue in the level set methods is to prescribe the normal velocity for the level set.

Active contours [19], were extensively used in segmentation [6,[START_REF] Zhang | Active contours with selective local or global segmentation: A new formulation and level set method[END_REF]. Their main idea is to evolve a parametric curve Γ into the object to be detected, see Fig: 1.3. The evolution of the curve is governed by the following equation:

∂Γ(r; t) ∂t = v(r; t), Γ(r; 0) = Γ 0 (r), (1.5)
where Γ 0 is an initial curve and v is the velocity of the curve evolution. Kass et al. [START_REF] Kass | Snakes: Active contour models[END_REF] Figure 1.3: The evolution of the parametric curve Γ into an object.

introduced the snakes model as a particular case of classical active contours approaches. It consists in minimizing the following energy:

J(Γ) = a b |Γ (r)| 2 dr + a b |Γ (r)| 2 dr + λ a b g(∇u(Γ(r)) dr, (1.6)
where the function g(.) plays the role of an edge detector. The first and second terms in the above energy control the smoothness of the curve while the third attracts the contours towards the edge of the object in the image u. The minimizer curve of the energy (1.6) corresponds to the sought contour. However, active contours approaches, including snakes model, heavily depend on the initial curve Γ 0 , and need special handing to changes in the topology of the evolving curves. These shortcomings were improved by Caselles et al. [24] by considering the geodesic active-contour model.

Introduction

Image restoration

Image restoration is a fundamental task in image processing. It consists in recovering an original image u from an observed one f which is degraded and contaminated by noise and usually modeled by problem (1.2). The purpose is to improve the image quality in order to get more reliable information which can be used in a wide range applications (medical imaging, astronomical imaging, to forensic science, etc). An equivalent mathematical formulation consists in approximating u solution of (1.2) by the minimizer of the following energy:

||T u -f || 2 Y . (1.7)
One of the major problems of (1.7) is its ill-posedness and minimizing (1.7) alone is very sensitive to data changes and yields in most cases highly oscillating minimizers. So there is a growing need for carefully dedicated theoretical and numerical methods for its treatment.

Regularization There is no way to avoid ill-posedness unless additional information about the solution is available in order to enlarge the solution space. Therein regularization techniques by building properties (regularity, curvature, gradient, etc) into the model and by expressing of constraints are required. The problem turned then to a new regularized and well-posed one which generally expressed as follows:

min u∈Y αR(u) + λ 0 ||T u -f || 2 Y , (1.8)
where α > 0 and λ 0 are two regularization parameters. The first part R(u) is the regularization function which represents the smoothing effect of the model (diffusion term) and has to be chosen according to desired properties of the solution. The minimization problem (1.8) is transformed by means of Euler-Lagrange equation into solving a Partial differential equation having the form of (1.12).

Another possibility to approach the issue of restoration is to work directly on the equations, without thinking of any energy. The image is viewed as a smoothed solution of Partial Differential Equation formally written as follows:

     ∂ t u -R(x, u, Du, • • • , D k u) = 0, in (0; T ) × Ω, u(x, 0) = f, in Ω, ∂ n u = 0, on (0; T ) × ∂Ω, (1.9)
where R is a k th -order differential operator. The earliest and the simplest equation used in the literature of image restoration is heat equation, i.e., R(u) = ∆u. It has a strong smoothing effect and it smooths noise, but also blurs important features such as edges which make them harder to identify. Then, nonlinear models were proposed in order to overcome the strong smoothness and to slow down diffusion at the edges. These models are based on the idea of applying diffusion, depending on the local properties of the image . They consist in incorporating the magnitude of the gradient in the diffusion process which is interpreted as a measure of edge likelihood. A prototype nonlinear isotropic diffusion model was introduced Introduction by Perona and Malik [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] and it corresponds to the following diffusion operator:

R(u) = ∇ • [g(|∇u| 2 )∇u],
where g nondecreasing diffusivity function g of |∇u| 2 such that g(0) = 1, g(s) > 0 and lim ∞ g(s) = 0. J. Weickert in [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF] extended the same technique to anisotropic processes. He suggested a diffusivity tensor D, depending on |∇u σ | 2 , instead of a scalar diffusivity g. The diffusion becomes variable in space and is driven according to the directional information contained in the image structure. The proposed anisotropic model corresponds to

R(u) = ∇ • [D(∇u σ )∇u].
Other alternatives extend the idea of Perona and Malik to higher order equations called "low curvature image simplifiers" and have the general form: R = ∇ • (g(m)∇∆u). The function m is chosen such that it measures some information about the image u. For example, You and Kaveh in [START_REF] You | Fourth-order Partial Differential Equations for Noise Removal[END_REF] proposed to use m = |∇u| which indicates the set jumps of u, Tumblin and Turk used the curvature in [START_REF] Tumblin | LCIS: A boundary hierarchy for detail-preserving contrast reduction[END_REF], Bertozzi and Greer [17] proposed the use m = |∆u| which locates the set jumps of the gradient of u, ∇u. to fill-in an incomplete/damaged region D ⊂ Ω based upon the image information available outside D (i.e., in Ω\D) in such a manner, that viewers can not detect the restored parts.

Image inpainting

State-of-the-art of image inpainting

Various mathematical approaches were proposed to solve the inpainting problem. They can be divided into three categories according to the idea used to mathematically model the problem. The first group, called texture-oriented methods, is concerned with the synthesis and the reconstruction of texture. The second and the third approaches are called geometryoriented methods and they are devoted to the recovering of the geometric part of images.

Texture image inpainting

Texture inpainting methods were triggered in the pioneering works of Efros and Leung [38] and Wei and Levoy in [START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF][START_REF] Wei | Texture synthesis over arbitrary manifold surfaces[END_REF] about texture synthesis. The overall idea of this approach is to synthesize a new image by producing texture from a given sample of a texture. The strategy has been extensively used for inpainting and the main algorithm can be summarized as follows: a damaged pixel x at the boundary of D is filled with the value of the known pixel y in Ω \ D, such that y being the center of the square patch P y which verifies:

P y = arg min Pz ,z ∈Ω \D d(P x , P z ),
where d(P x , P z ) is the Sum of Squared Differences (SSD):

d(P x , P z ) = i,j |P x (i, j) -P z (i, j)| 2 , (1.10)
and the indices (i, j) describe pixels in the patch. Therefore, the missing gap is filled-in recursively pixel by pixel. The main shortcoming of this algorithm is its computational cost. Criminisi et al. [34] improved the previous work by considering Patch-based inpainting algorithm. The model focuses on a patch-based filling approach as opposed to pixel-based ones as in Efros and Leung [38] which makes the method considerably faster. Their approach is also able to correctly inpaint straight edges, that are not in large scale, by promoting damaged pixels located near edges to be inpainted prior others pixels. However, it fails to reproduce large-scaled geometric information, such long edges.

Introduction

E(ϕ) = D Br |u 0 (x + y) -u 0 (ϕ(x) + y)| 2 dy dx,
where B r = B r (0) is the ball of center 0 and radius r and it represents the patch. The unknown variable is the correspondence map ϕ. The integral of the quadratic distance measures the visual closeness of two patches. This model has been the subject of further theoretical analysis and improvement by Aujol et al. [8] by proposing extensions and proving the existence of a solution.P. Arias, V. Caselles, G. Facciolo [4]proposed more easily computable variational approach for the problem of non-local image inpainting based on non-convex models. The idea consists in minimizing a non-local means type energy with dynamic weights. The latter are determined by minimizing an entropy term and they measure the similarity between patches centered in the inpainting domain D and in its complement.

Level lines approaches

The spirit of these approaches is the pioneering work of Nitzberg et al. in [START_REF] Nitzberg | Filtering, Segmentation, and Depth[END_REF]. It takes its origin from the ideas of Gestalt psychology and Kanizsa [START_REF] Kanizsa | Grammaire du Voir[END_REF] about the ability of the visual system to complete partially occluded edges following a principle of good continuation, called amodal completion principle. The method consists in identifying points, called T -junctions, where visible edges intersect the occluding objects and performing good continuation between each pairs of them likely belong to the same edge. The aim is to join them by an optimal curve with respect to appropriate criterion. Based on Nitzberg and [START_REF] Nitzberg | Filtering, Segmentation, and Depth[END_REF] work, Masnou and Morel [START_REF] Masnou | Filtrage et désocclusion d'images par méthodes d'ensembles de niveau[END_REF][START_REF] Masnou | Disocclusion a variational approach using level lines[END_REF][START_REF] Masnou | Level-lines based disocculusion[END_REF] considered another representation of the image in term of the level set:

X r u = {x ∈ Ω; u(x) ≥ r},
and its correspondent boundaries Γ r := {x ∈ Ω | u(x) = r}, called level lines. The family of X r u gives a complete representation of the image u and the level lines represent its edges, which transform the problem to level lines completion one. In this case, T -junctions are defined as points where ∂D intersects the level lines and two T -junctions T i and T j are compatible if both are associated with the same level set and having the same orientations θ i and θ j of ∇u. The underlying idea of the proposed level line approach is to search for an Introduction optimal curve Γ called a completion curve that joins two compatible T -junctions T 1 and T 2 , among the minimizers of the following Euler elastica energy:

Γ (α + βκ p )dH 1 + (θ 1 , N 1 ) + (θ 2 , N 2 ),
where a, b are positive constants, p ≥ 1, κ is the curvature of Γ, and the last two terms denote the angles between θ i and N i , the normal to γ at T i (for i = 1, 2). Then by considering all the where F r denotes the family of completion curves associated with the level sets X r . For the energy minimization for p = 1, the authors used a particular practical implementation based on dynamic programming to compute the completion curves. The results obtained in this work were completed with a theoretical issues related to the energy (1.21) in [START_REF] Masnou | Disocclusion a variational approach using level lines[END_REF] for p = 1 and in [1,[START_REF] Masnou | On a variational theory of image amodal completion[END_REF] for p > 1.

Partial Differential Equations: Formulation and overview

In the last few years, a number of PDE-based inpainting models have appeared. Some models are written down directly, based on some mathematical understandings on the properties of diffusion PDEs (the heat equation, the Cahn-Hilliard equation, etc.) or in analogy with linear and nonlinear filtering [START_REF] Oliveira | Fast digital image inpainting[END_REF], by the construction of gradient flows [40,[START_REF] Haehnle | Mumford-Shah-Euler flow with sphere constraint and applications to color image inpainting[END_REF] or by active contours modeling [39,[START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF]. Others are formulated as an adequate optimization problem having the form (1.8). We better explain in the next section PDE-based approaches and we give an overview on them.

The underlying idea of PDE-based methods is to interpolate image data across all boundaries of the missing area D via a process modeled by the following k th -order evolutionary PDE:

∂ t u -R(x, u, Du, • • • , D k u) = 0, in D, u = f, on ∂D, (1.11)

Introduction

where R is a k th differential operator. In the noisy case, various approaches state the problem in all the image domain Ω as follows:

∂ t u -R(x, u, Du, • • • , D k u) + λ D (u -f ) = 0, in Ω, ∂ n u = 0, on ∂Ω, (1.12)
where λ D = λ 0 χ D and equation (1.12) is equivalent to (1.11) when λ 0 is large. The choice of the operator R is inspired by phenomenological modeling like Navier-Stokes and fluid dynamics in [13,14], or inspired from phase field models like Cahn-Hilliard equation [15,16,18]. In other instance, it comes from the idea of propagating local information with smoothness constraints, in analogy with physical phenomena like heat propagation in physical structures by a diffusion process. These diffusion PDEs take their origin from the following general continuity equation:

∂ t u = ∇ • j, (1.13)
where j is expressed by Fick's law as follows:

j = -D∇u,
and D is a diffusion tensor intended to be positive definite symmetric matrix which encodes the properties of the diffusion. In case where j and ∇u are parallel, equation (1.13) is called isotropic and we may replace the diffusion tensor by a positive scalar-valued diffusivity. In general, i. e., where j and ∇u are not parallel, the equation is called anisotropic. Various investigations are focusing on the choice of the diffusion tensor D: a positive diffusion tensor D(x), i. e., depends only on x and not on the image u, leads to linear, isotropic or anisotropic PDEs which are the simplest PDE methods and mostly used for image smoothing. However, the strong smoothing character makes them less exploited in image inpainting because they produce blurring and smooth edges. Then, other alternatives are focusing on adaptive nonlinear diffusion methods and they are especially oriented towards understanding the mathematical properties of the diffusion and to drive it by choosing D depending on specific and local properties of the image u.

Overview on PDE methods

A large number of relevant PDE-based approaches can be found in the literature. This area is still very active and it is technically difficult to cover all the great approaches in a single work . Next, we mention some contributions that, in our opinion, have a major impact on this area. We also emphasize the difference between second-and higher-order models in inpainting.

Second-order PDEs

Second-order PDEs have the advantage of being well-established theoretically and easy to solve numerically. However, they do not usually perform well in inpainting chiefly when the size of the hole in the image is important. In fact, the major drawbacks are either the disconnection of level lines over large distances (Connectivity Principle), either their smooth propagation into the damaged domain or the inability of reproducing features of higher order Introduction (curvature, corners,...) due to the lack of information (second-order). Let us now start with the presentation of some existing second-order inpainting models.

Smooth inpainting via Harmonic extension [31] The aim is the Harmonic extension of the known image f on Ω\D by a function u which solves the following PDE: ∆u = 0, in D and u| ∂D = f, (1.14) or to combine, for noisy images, inpainting and denoising in a single equation by focusing a function u ∈ H 1 (Ω) which solves the following PDE:

-∆u + λ D (u -f ) = 0, in Ω, ∂ n u = 0, on ∂Ω. (1.15)
Equations (3.2) and (3.1) are isotropic linear diffusion PDEs and they are good candidates for smooth images. However, the strong smoothness property makes them unable to restore occluded edges which represent a crucial object in image processing problems.

Topological gradient method Auroux and Masmoudi in [9,10] used techniques that have been developed for the inverse conductivity problem in the context of crack detection. The missing contours in the domain D were identified as cracks σ which minimize the energy:

j(σ) = 1 2 u d -u n 2 L 2 (Ω) ,
where u d and u n solve:

         ∆u d = 0, in D \ σ, u d = f, in ∂D, ∂ n u d = 0, on σ, u d = f, in Ω \ D, and
         ∆u n = 0, in D \ σ, ∂ n u n = ∂ n f, in ∂D, ∂ n u n = 0, on σ, u n = f, in Ω \ D, (1.16)
respectively. Edges are localized using the topological asymptotic analysis and defined by the most negative points of the topological gradient. The method provides a promising numerical results and the computing time was very short comparing to another approaches. The only shortcoming of this work is the getting of disconnected contours in the missing part chiefly when the image is noisy. As the solution is harmonic in each subdomain of D, a discontinuous contour would lead to some blurred and smooth image. The authors considered a thresholded version of the topological gradient in order to manage the problem of identifying connected contours. However, this strategy is not very relevant because a bad choice of the threshold parameter does not give sharp edges and does not distinguish between them and noisy points.

In order to overcome this limitations, minimal paths and fast marching techniques were used in [35].

Total variation and related works Total variation (TV) method was originally developed by Rudin, Osher, and Fatemi [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] for image denoising and it was exploited in inpainting Introduction by Chan and Shen [31]. It consists in regularization term, proven to be equivalent to a penalization of the length of level lines in D, and minimizing the following regularized energy:

F (u) = Ω |∇u|dx + Ω λ D (u -f ) 2 dx, (1.17)
over the space of functions of bounded variation BV (Ω). The model allows the recovering of sharp geometric structures and promote discontinuities. However, it is unable to joint level lines and edges across large domain D because of the penalization term of the length. Another shortcoming of this method is that it does not always keep the direction of the isolines across the boundary of the damaged parts. This work was improved by Chan and Shen proposed in [30] by incorporating, in its Euler-Lagrange equation, a diffusivity function g depending on the curvature of level lines and then solve the following Curvature Driven Diffusion (CDD) equation:

   -∇• g(|κ|) |∇u| ∇u + λ D (u -f ) = 0, in Ω, ∂ n u = 0, on ∂Ω, (1.18)
where κ = ∇• ∇u |∇u| is the scalar curvature and g : R + -→ R + is a continuous function.

The new information about the geometry makes possible to extend edges over large distance and to fulfill the connectivity principle. In other alternative [32,[START_REF] Strong | Adaptive total variation minimization image restoration[END_REF][START_REF] Strong | Relation of Regularization Parameter and Scale in Total Variation Based Image Denoising[END_REF], a weighted total variation (TV) functional Ω α(x)|∇u|dx was proposed, leading to adaptive choice of α.

Mumford-Shah

The Mumford-Shah model appeared in the works of Tsai, Yezzi, Willsky [START_REF] Tsai | Curve evolution implementation of the Mumford-Shah functional for image segmentation, denoising, interpolation, and magnification[END_REF] and it was originally designed for the segmentation and denoising problems. Afterward, it was considered in [39] for image inpainting purpose and it consists in minimizing the following energy functional:

γ 2 Ω\Γ |∇u| 2 dx + 1 2 Ω λ D (u -f ) 2 dx + H 1 (Γ), (1.19)
where Γ denotes the collection of edges, H 1 denotes the one dimensional Hausdorff measure, which generalizes the length notion for regular curves. The shortcomings of this model are parallel to those of the TV inpainting one. In fact, it favors straight edges and cannot connect isolines across large distances due to the penalization term of the length. As a second-order model, it is then unable to reconstruct curved objects.

Higher-order PDEs

The shortcomings of second-order methods gave rise to a new class of higher-order diffusion models which in general perform better. In fact, whenever the image is contaminated by noise, they damp the oscillations and high frequencies (noise) faster than second-order based diffusion models. In addition, they are richer than second order ones and have the benefits of connecting isolines across missing parts while preserving the curvature (of third-order). Moreover, the use of boundary conditions for both the solution u(x) and its derivatives gives Introduction a supplementary information on the isoline directions and allows matching edges across large distances.

Several contributions are focusing on higher-order PDEs by incorporating higher-order derivatives or a sophisticated combination of higher-and second-order ones.

Bertalmio approach & fluid interpretation The pioneering works within PDE-based approaches is that of Bertalmio et al. [14] where the authors proposed to propagate both the gradient direction (geometry) and the gray-scale values of the image inside the region to be filled by solving the PDE :

∂ t u -∇ ⊥ • ∇∆u = 0, in D, u = f, on ∂D, (1.20)
where ∇ ⊥ denotes the perpendicular gradient (-∂ y , ∂ x ), and ∆ is the Laplace operator. In [13], Bertozzi et al. exploited ideas from classical computational fluid dynamics in the image inpainting problem. Their approach is directly based on the Navier-Stokes equations for fluid dynamics and is given by the following equation:

w = ∆u, ∂ t w + v • ∇w = σ∇.(g(|∇w|)∇w),
where the velocity field v = ∇ ⊥ u and the function σ allows for anisotropic diffusion of the smoothness w. In analogy with fluid dynamics, the image intensity is token as a "stream function" for two-dimensional incompressible flow and its Laplacian plays the role of the vorticity of the fluid.

Euler's Elastica Model T. Chan et al. proposed in [29] a slightly different model for the Euler-elastica functional. They considered the following functional:

u -f L 2 (Ω) + K (a + bκ 2 ) dH 1 (x), (1.21)
where K is a closed regular subset of Ω and κ = ∇ • (∇u/|∇u|) is the curvature of level sets γ r := {x ∈ K | u(x) = r}. It is a higher-order variational model where the regularization term combines the total variation, sensitive to the length of the isolines, and the square of the curvature, which favors curves lines than straight ones. Minimizing the Euler-elastica leads to a highly nonlinear PDE and therefore its numerical solution is a non trivial task and was the subject of many investigations [1,33,39,[START_REF] Masnou | Filtrage et désocclusion d'images par méthodes d'ensembles de niveau[END_REF].

In [21,33], a fast and efficient numerical algorithms were proposed to solve minimization problem of Euler's elastica energy. They consist in reformulating the minimization problem as a constrained one by an operator splitting method as follows:

min u,p,n Ω (a + b(∇ • n) 2 |p|dx + η 2 Ω (u -f ) 2 dx
γ 2 Ω\Γ |∇u| 2 dx + 1 2 Ω λ D (u -f ) 2 dx + Γ (α + βκ 2 )ds, (1.22)
where α and β are non-negative regularization parameters. Its numerical realization was based on the Γ-convergence approximations originally conjectured of Ambrosio and Tortorelli [2,3], and De Giorgi [43]. The minimization leads to a fourth-order PDE, facing some difficult issues in the numerical computation, such as cost, convergence, stability, etc.

Combined first and second order derivatives Recently, the authors in [START_REF] Papafitsoros | A combined first and second order variational approach for image reconstruction[END_REF] have introduced a new high-order nonlinear model called T V -T V 2 image inpainting and it consists in minimizing the following energy:

||u -T f || L 2 (Ω) + αT V (u) + βT V 2 (u), (1.23)
where α and β are non-negative regularization parameters chosen empirically, T V (u) and T V 2 (u) are the total variations of u and ∇u, respectively. It is higher-order extension of the well-known Ruden-Osher-Fatemi functional [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] (total variation minimization). The underlying idea is to extend the works in [27] (see also, e.g., [START_REF] Hinterberger | Variational methods on the space of functions of bounded Hessian for convexification and denoising[END_REF][START_REF] Lysaker | Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[END_REF][START_REF] Lysaker | Iterative image restoration combining total variation minimization and a second-order functional[END_REF][START_REF] Scherzer | Denoising with higher order derivatives of bounded variation and an application to parameter estimation[END_REF]) where energy functionals were proposed and they are straightforward convex combinations of first-and second-derivatives having the following general form:

u -f L 2 (Ω) + Ω G 1 (∇u) dx + Ω G 2 (∇ 2 u) dx, (1.24)
where G 1 (•) and G 2 (•) are given functions.

Cahn-Hilliard inpainting Bertozzi, Esedoglu and Gillette proposed in [16] the Cahn-Hilliard equation for binary image inpainting. It is a semilinear fourth-order PDE originally introduced in material sciences by John W. Cahn and John E. Hilliard in [23] and it describes the evolution of an interface separating two stable states. The equation is the following:

Introduction

Organization of the Thesis

This thesis is organized into five chapters and a conclusion.

Chapter 1 is a general introduction to the problem treated in the thesis. Chapter 2, concerns the work appeared in [START_REF] Kallel | The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting[END_REF]. It is dedicated to modeling the inpainting problem as a suitable Cauchy one, where the image information in the incomplete/damaged region D is completed by the Cauchy problem solution. According to the type of the image (smooth or containing edges), a linear and nonlinear models will be investigated and will be solved as a Nash-game.

Chapter 3 deals the work presented in [12] which is dedicated to the use of efficient strategies to select locally and adaptively the regularization parameters for second-order variational models in image inpainting problems in order to account for the geometric specificities of the image.

In Chapter 4, we present the works in [12]. We consider simple and effective fourth-order variational model for solving image inpainting problem. Then, we perform a tractable control of the regularization parameter on posterior way combined with mesh adaptation techniques.

In Chapter 5, we propose a fourth-order models in image inpainting and restoration problems. The model combines second-and fourth-order derivatives, in analogy with Eulerelastica energy. Then, we introduce a multi-scale approach, rendered by a variable diffusion coefficients in the same spirit of the previous chapter.

".. although nature begins with the cause and ends with the experience we must follow the opposite course, namely begin with the experience and by means of it end with the cause."

Leonardo da Vinci (1452-1519).

Chapter 2

The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting

The purpose of this chapter is the formulation of the inpainting problems as a nonlinear Cauchy problem and the use of the game theory framework to solve it. It is based on the article [START_REF] Kallel | The Cauchy problem for a nonlinear elliptic equation: Nash-game approach and application to image inpainting[END_REF] which is summarized below.

Summary

We consider an image f defined on a domain Ω in R d (d = 2, 3), usually a rectangular domain with piecewise smooth boundary ∂Ω. Let D ⊂ Ω be a damaged subregion in Ω where the information is not available. Local PDEs inpainting methods consider the diffusion of the image information available in the boundary ∂D into the interior of D. However, two typical cases are considered in inpainting depending on the position of damaged area D in the image domain Ω. The first corresponds to the case where D is totally inside Ω (i.e., ∂D ∩ ∂Ω = ∅). In this case, inpainting can be achieved by solving the Dirichlet problem (1.11) or the problem (1.12). The second is more complicated and it corresponds to the case where ∂D ∩ ∂Ω = ∅. In such situation, we consider the following decomposition of ∂D:

Γ c Γ i = ∂D, Γ i ⊂ ∂Ω, Γc Γi = ∅, Γc ∂Ω = ∅. (2.1)
Therefore, solving inpainting problem by solving problem (1.11), which is posed locally in D, involves the use of Dirichlet Boundary conditions on all the boundary ∂D, which are not accessible on Γ i . In addition, the choice of the homogeneous Neumann condition, which comes from restoration problems, in (1.12) it not reliable and this condition is not fulfilled in some situations. Therefore, there is need for specialized inpainting methods and unusual boundary conditions that are more relevant for the case where ∂D ∩ ∂Ω = ∅. Disposing only boundary conditions on one side of ∂D is the natural setting to formulate inpainting as

Inpainting as a Cauchy problem

an extrapolation problem which will be achieved by solving the following inverse nonlinear boundary value problem:

       ∇•[k(|∇u| 2 )∇u] = 0, in D, u = f, on Γ c , k(|∇u| 2 )∇u • n = φ, on Γ c , (2.2)
where the diffusivity function k(•) fulfills some assumptions. n is the unit outward normal to ∂D, f ∈ H Inverse problem (2.2) is ill-posed in the Hadamard sense [START_REF] Hadamard | Lectures on Cauchy's Problem in Linear Partial Differential Equations[END_REF]. In fact, existence of solutions for arbitrary Cauchy data f and φ is not guaranteed and depends on the compatibility of f and φ. They are said to be compatible if the Cauchy problem (2.2) has a solution. We also say that f and φ are consistent for (2.2). We proved in the nonlinear case that there exists always a compatible data. More precisely, for fixed f ∈ H 1 2 (Γ c ), we prove that the set of data φ ∈ (H

1 2 00 (Γ c )) which are compatible with f is dense in (H 1 2
00 (Γ c )) . As regards the uniqueness, we prove that, if the Cauchy data f and φ are compatible, then the Cauchy problem (2.2) has a unique solution in H 1 (D). Nevertheless, Ill-posedeness of the Cauchy problem (2.2) is also related to the instability of the solution, even it exists, with respect to small perturbation of the Cauchy data. For that reason, classical numerical methods are usually inappropriate because they are unstable. Hence, there is a growing need for carefully stabilized and dedicated computational methods for the numerical treatment of the Cauchy problem (2.2). In this work, we give a Nash-game formulation of problem (2.2). The work follows the one in introduced in [START_REF] Habbal | Data completion problems solved as Nash games[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic cauchy problems[END_REF] for solving the Cauchy problem for a linear operator, i.e., k(•) depends only on x ∈ D.

Loosely speaking, the Cauchy problem is formulated as a two-player Nash-game with the following two costs:

     J 1 (η, τ ) = 1 2 (k(|∇u 1 (η)| 2 )∇u 1 (η) • n -φ|| 2 (H 1 2 00 (Γc)) + 1 2 ||u 1 (η) -u 2 (τ ) 2 H 1 2 (Γ i ) , J 2 (η, τ ) = 1 2 ||u 2 (τ ) -f || 2 H 1 2 (Γc) + 1 2 ||u 1 (η) -u 2 (τ )|| 2 H 1 2 (Γ i ) , (2.3)
where (η, τ ) ∈ (H

1 2 00 (Γ i )) × H 1 2 (Γ i ).
The first player u 1 (η) is given the known Dirichlet data f and uses the Neumann condition prescribed over the inaccessible Γ i part of the boundary ∂D as strategy variable . The second player u 2 (τ ) is given the known Neumann data φ, and plays with the Dirichlet condition prescribed over the inaccessible boundary. The two players solve in parallel two associated direct boundary value problems (SP 1) and (SP 2), respectively. The goal is to achieve an optimum which represents an equilibrium situation to all players, in the sense that none of the players has an interest to optimize its cost function. This optimum is called "Nash equilibrium" and is fundamental concept in game theory. For the particular case where k(•) = 1, the authors in [START_REF] Habbal | Data completion problems solved as Nash games[END_REF][START_REF] Habbal | Neumann-Dirichlet Nash strategies for the solution of elliptic cauchy problems[END_REF] proved that there always exists a unique Nash equilibrium, which turns out to be the reconstructed data when the Cauchy problem has a solution. They also proved that the completion process by the Nash-game approach is stable with respect to noisy data. For the nonlinear case, i.e., k(•) depends on u, and when the Cauchy problem (2.2) has a solution u, we prove that the Inpainting as a Cauchy problem

pair (η c , τ c ) = (k(|∇u| 2 )∇u • n| Γ i , u| Γ i )
is unique "Pareto optimal Nash equilibrium" for the two-players game involving the costs functionals J 1 and J 2 .

To compute the Nash equilibrium, we use a classical alternating minimization algorithm with relaxation [START_REF] Uryas | On relaxation algorithms in computation of noncooperative equilibria[END_REF], also referred to as the inertial Nash equilibration process. The minimization claimed the computation of the gradients of J 1 and J 2 by means of the adjoint state method and involved the resolution of the direct problems (SP 1) and (SP 2). From classical calculus of variations, we prove the existence and uniqueness of solution in H 1 (D) of the direct problem (SP 1) and the proof for (SP 2) is similar. We also linearized it and we proved the convergence of the sequence of solutions of the linearized problems to the solution of the nonlinear direct problem. Finally, we present numerical results carried out on different images and they are obtained using the Nash-game algorithm. We also present numerical results which illustrate the stability, at least numerically, of our method with respect to noisy data. To demonstrate the efficiency of the proposed method from data completion viewpoint, we make a comparison with the Kozlov algorithm [START_REF] Kozlov | An iterative method for solving the Cauchy problem for elliptic equations[END_REF] which is an alternative method for solving Cauchy problems. In particular, many variational models have appeared in the literature. These models give rise to partial differential equations for which Dirichlet boundary conditions are usually used. The basic idea of the algorithms that have been proposed in the literature is to fill-in damaged regions with available information from their surroundings. The aim of this work is to treat the case where this information is not available in a part of the boundary of the damaged region. We formulate the image inpainting problem as a nonlinear Cauchy problem. Then, we give a Nash-game formulation of this Cauchy problem and we present different numerical experiments using the finite-element method for solving the image inpainting problem.

1. Introduction. The Cauchy problem for elliptic equations, which is generally an ill-posed problem, remains a challenge for the analysis of inverse problems. The question on the ill-posedness was first studied by Hadamard [21]. He showed that the Cauchy problem associated with the Laplace operator is ill-posed. There are two issues related to ill-posedness. First, existence of solutions for arbitrary Cauchy data is not guaranteed. The Cauchy data must be compatible (or consistent), i.e., the corresponding Cauchy problem has a solution. Second, even if such a solution exists, it is unstable with respect to small perturbation of the data. Since the work of Lavrentev [32,33], Maz'ya and Havin [38], and Tarkhanov [40] there has been a renewed interest in this classical ill-posed problem, see e.g. the recent works [1,5,7,8,19,22]. For this ill-posed problem, the classical numerical methods are usually inappropriate because they are unstable. So there is a growing need for carefully stabilized and dedicated computational methods for its numerical treatment. Recently, in [18,19] the authors formulated the Cauchy problem for a linear elliptic operator as a Nash game. They proved that there always exists a unique Nash equilibrium, which turns out to be the reconstructed data when the Cauchy problem has a solution. They also proved that the completion process by the Nash-game approach is stable with respect to noisy data.

In the last decade, analytical and numerical approaches for nonlinear elliptic Cauchy problems were also developed. The iterative method used by Kozlov et al.

854

Moez Kallel, Maher Moakher and Anis Theljani [30] was extended in [4] for a nonlinear equation arising from a model in glaciology. Leitao et al. proposed in [14,31] a level-set approach for solving a nonlinear elliptic Cauchy problem. They formulated the Cauchy problem as a minimization problem with Tikhonov regularization. Such nonlinear problems can be linearized by the Cole-Hopf transformation [12]. We also mention that Ly [35] have developed a new method for solving Cauchy problems for the p-Laplace operator.

In this paper, we consider a nonlinear Cauchy problem arising in image inpainting which consists in reconstructing lost or deteriorated parts of an image. We refer the reader interested in the image inpainting problem and its applications to [6,10,11,9,15,37]. Different techniques can be applied to solve this problem, in particular, partial differential equations (PDEs) are widely used and are proven to be efficient.

Let Ω ⊂ R d (d = 2, 3) denote the entire image domain, the problem is to fillin image information in the incomplete/damaged region D ⊂ Ω based upon the image information available outside D (i.e., in Ω\D). We assume that ∂D is sufficiently smooth and composed of two connected components Γ c and Γ i . When the information is available near all the boundary ∂D = Γ i ∪ Γ c , it can be used as Dirichlet boundary conditions for the partial differential equation that propagates the information inside D such as in [10,11]. The aim of this work is to treat the case where this information is not available in the part Γ i . The inpainting problem is formulated as a Cauchy problem which is linear in the case of smooth images and nonlinear for images containing edges. In both cases, we consider the following inverse boundary-value problem for the unknown image function u:

(1)

     ∇•[k(|∇u| 2 )∇u] = 0, in D, u = f, on Γ c , k(|∇u| 2 )∇u • n = φ, on Γ c ,
where n is the unit outward normal to ∂D, f ∈ H (2)

                 (a) α 1 ≥ k(r) ≥ α 0 > 0, ∀r ∈ R + ; (b) k (r) ≤ 0, ∀r ∈ R + ; (c) k(r) + 2rk (r) ≥ α 2 > 0, ∀r ∈ R + ; (d) α 3 > 0, (B(ω)X, X) ≥ α 3 |X| 2 , ∀X ∈ R d , where B(ω) = k(|ω| 2 )I d + 2k (|ω| 2 )ω ⊗ ω, ω ∈ R d and I d is the identity in R d .
We remark that the linear case corresponds to k(•) ≡ 1 which fulfills all conditions listed above.

We note that (a), (b), (c) and (d) are the natural sufficient conditions to guarantee the solvability of the direct nonlinear boundary-value problem:

(3)

     ∇•[k(|∇u| 2 )∇u] = 0, in D, u = f, on Γ c , k(|∇u| 2 )∇u • n = ϕ, on Γ i , where f ∈ H 1 2 (Γ c ) and ϕ ∈ (H 1 2
00 (Γ i )) are given functions.

Inverse Problems and Imaging

Volume 9, No. 3 (2015), 853-874

The remainder of this paper is organized as follows: In Section 2, the existence and uniqueness of solution in H 1 (D) of the direct problem were proved. We construct a sequence of solutions of the linearized problems which converges in H 1 -norm to the solution of the nonlinear problem. In Section 3, a density result for compatible data of the Cauchy problem and uniqueness of its solution have been proven. Section 4 gives the Nash-game formulation of the Cauchy problem and presents a new algorithm for its numerical solution. In Section 5, we formulate the image inpainting problem as a Cauchy problem. We illustrate the efficiency and robustness of the proposed method by treating a wide range of situations.

2. The direct problem. Nonlinear problems similar to (3) have been studied in several works [2,26,27,28,29]. In [23] Hazanov introduced an inverse coefficient problem for such operator and studied the direct problem with homogeneous Dirichlet condition. He also linearized the problem and proved the convergence of the sequence of solutions of the linearized problems to the solution of the nonlinear one. We will study in this section the direct nonlinear boundary-value problem (3). We then give the linearized problem and prove the convergence, in H 1 -norm, of its solution to the solution of the nonlinear direct problem.

Let

u * ∈ H 1 (D) such that u * | Γ c = f
. Then, it is easy to see that problem (3) can be written in the following form (with unknown function w):

(4)      ∇•[k(|∇(w + u * )| 2 )∇(w + u * )] = 0, in D, w = 0, on Γ c , k(|∇(w + u * )| 2 )∇(w + u * ) • n = ϕ, on Γ i .
A weak form of problem ( 4) can be stated as follows:

(5)

   Find w ∈ W = {w ∈ H 1 (D); w| Γ c = 0} such that D k(|∇(w + u * )| 2 )∇(w + u * ) • ∇v dx = Γ i ϕv ds, ∀v ∈ W.
To find the solution w of the nonlinear problem (5), we start from an initial guess, w 0 , and we construct the sequence (w n ) n>0 of solutions of the linearized problems: (6) Find

w n ∈ W = {w ∈ H 1 (D); w| Γ c = 0} such that A(w n-1 ; w n , v) = Γ i ϕv ds -(w n-1 ; v), ∀v ∈ W, where (7)
     A(w n-1 ; u, v) = D k(|∇(w n-1 + u * )| 2 )∇u • ∇v dx, ∀u, v ∈ W, (w n-1 ; v) = D k(|∇(w n-1 + u * )| 2 )∇u * • ∇v dx, ∀v ∈ W .
The proof of convergence of the sequence of solutions of the linearized problems (6) can be easily adapted from the proof given by Hasanov [23] for the homogeneous case, i.e., u = 0 on Γ c . The only difference is that one has to consider at each iteration the linear form (w n-1 ; v). While conditions (a), (b) and (c) given in (2) imply uniqueness of solution in the homogeneous case, we have added condition (d) to guarantee uniqueness in our case. For the proof of convergence, we need the following Lemma. Lemma 2.1. Let J be the functional on W defined by

(8) J(v) = 1 2 D |∇(v+u * )| 2 0 k(r) dr dx.
Then, for all w n-1 and w n in W we have:

1 2 A(w n-1 ; w n , w n ) + (w n-1 ; w n ) -1 2 A(w n-1 ; w n-1 , w n-1 ) -(w n-1 ; w n-1 ) -J(w n ) + J(w n-1 ) ≥ 0. Proof. Let (9) Q(t) = t 0 k(r) dr.
Then, it is easy to verify that the function Q is concave. In fact, we have:

Q (t) = k (t) ≤ 0, ∀t > 0.
From the property of concave differentiable functions, we have:

Q (t 1 )(t 2 -t 1 ) -Q(t 2 ) + Q(t 1 ) ≥ 0, ∀t 1 , t 2 .
By taking

t 1 = |∇(w n-1 + u * )| 2 and t 2 = |∇(w n + u * )| 2
, we obtain the proof. (ii) Let be (w n ) n ⊂ H 1 (D) the sequence of solutions of the linearized problems (6). Then, the sequence

(w n + u * ) n converges, in H 1 -norm as n → ∞, to the weak solution u ∈ H 1 (D) of the nonlinear problem (3).
Proof. We start by proving (i). Problem ( 4) is the Euler-Lagrange equation for the minimization problem (10)

min w∈W L(w),
where L(w) = J(w) -

Γ i
ϕw ds and the space W is defined in (6). In fact,

L (w), h = d dt {L(w + th)} t=0 = D k(|∇(w + u * )| 2 )∇(w + u * ) • ∇h dx - Γ i ϕh ds, ∀h ∈ W.
Therefore, the functional L(w) is Gâteaux differentiable and hence by using Green's formula we obtain the equivalence between the two problems. Now, we will prove that minimization problem (10) has a unique solution w ∈ W . First, we prove that L(w) is continuous. Let h ∈ W , then we have

L(w + h) -L(w) = 1 2 D Q(|∇(w + u * + h)| 2 ) -Q(|∇(w + u * )| 2 ) dx - Γ i ϕh ds,
where Q is defined in (9). From the Mean-Value Theorem applied for the function Q we conclude that there exists 0 < δ < 1 such that

Q(|∇(w + u * + h)| 2 ) -Q(|∇(w + u * )| 2 ) = Q (R) |∇(w + u * + h)| 2 -|∇(w + u * )| 2 , where R = (1 -δ)|∇(w + u * )| 2 + δ|∇(w + u * + h)| 2 .
Inverse Problems and Imaging Volume 9, No. 3 (2015), 853-874

The Cauchy problem for a nonlinear elliptic equation 857 Therefore,

L(w + h) -L(w) = 1 2 D Q (R) |∇(w + u * + h)| 2 -|∇(w + u * )| 2 dx - Γ i ϕh ds = 1 2 D k(R) |∇(w + u * + h)| 2 -|∇(w + u * )| 2 dx - Γ i ϕh ds = 1 2 D k(R) -2∇(w + u * ) • ∇h + |∇h| 2 dx - Γ i ϕh ds
Due to condition (a), Poincaré's inequality and the continuity of the trace operator from H 1 (D) to L 2 (Γ i ), for all h ∈ W we have

|L(w + h) -L(w)| ≤ c 1 ∇(w + u * ) L 2 (D) + c 2 ϕ L 2 (Γ i ) + c 3 h H 1 (D) h H 1 (D) ,
where c 1 , c 2 and c 3 are positive constants. This inequality shows that the functional L is continuous in W . We now prove that the functional L(w) is strongly convex, which means that it is coercive. We have

L (w; v, h) = d dt D k(|∇(w + u * + th)| 2 )∇(w + u * + th) • ∇v dx t=0 = D [k(|∇(w + u * )| 2 ) + 2k (|∇(w + u * )| 2 )|∇(w + u * )| 2 ]∇v • ∇h dx.
Using condition (c) and applying Poincaré's inequality, we obtain

L (w; h, h) ≥ α 2 h H 1 (D) , ∀h ∈ W,
which implies the strong positivity of the second Gâteaux derivative. Then, we have L(w) is a strongly convex function (hence coercive) and the minimization problem (10) admits a unique solution w ∈ W . Now, it is easy to see that u = w + u * in H 1 (D) is a weak solution to the direct problem (3) which clearly depends only on f = u * | Γ c . To prove uniqueness, let u 1 and u 2 be two weak solutions of (3) and set U = u 1u 2 and u t = u 1 + t(u 2u 1 ). Then we have

k(|∇u 1 | 2 )∇u 1 -k(|∇u 2 | 2 )∇u 2 = 1 0 d dt k(|∇u t | 2 )∇u t dt = M (u 1 , u 2 )∇U,
where

M (u 1 , u 2 ) = 1 0 B(∇u t ) dt,
and the matrix-valued function B(•) is given in (d). Therefore, U is a solution of the following problem (11)

     ∇•[M (u 1 , u 2 )∇U ] = 0, in D, U = 0, on Γ c , M (u 1 , u 2 )∇U • n = 0, on Γ i .
For Γ c with positive Hausdorff measure, it is well known that (11) For (ii), let w n be the solution of the linearized problem (6). Using the coercivity of the bilinear form A(w n-1 ; •, •) and condition (a) we get

α 0 w n -w n-1 2 ≤ A(w n-1 ; w n -w n-1 , w n -w n-1 ) = A(w n-1 ; w n , w n -w n-1 ) -A(w n-1 ; w n-1 , w n -w n-1 ) = A(w n-1 ; w n , w n ) -A(w n-1 ; w n , w n-1 ) -A(w n-1 ; w n-1 , w n ) + A(w n-1 ; w n-1 , w n-1 ) = A(w n-1 ; w n-1 , w n-1 ) -A(w n-1 ; w n , w n-1 ) -A(w n-1 ; w n-1 , w n ) + Γ i ϕw n ds -(w n-1 ; w n ) = A(w n-1 ; w n-1 , w n-1 ) -A(w n-1 ; w n , w n-1 -w n ) -A(w n-1 ; w n , w n ) -A(w n-1 ; w n-1 , w n ) + Γ i ϕw n ds -(w n-1 ; w n ) = A(w n-1 ; w n-1 , w n-1 ) -2 Γ i ϕw n-1 ds + 2 (w n-1 ; w n-1 ) -A(w n-1 ; w n , w n ) -2 (w n-1 ; w n ) + 2 Γ i ϕw n ds.
Using the inequality in Lemma 2.1, we obtain

(12) α 0 w n -w n-1 2 ≤ 2(L(w n-1 ) -L(w n )).
Then, we proved that sequence (L(w n )) n is monotone decreasing and bounded below, so it converges. From inequality (12) we deduce that (w n ) n is a Cauchy sequence. Therefore, it converges to some w in H 1 (D). Passing to the limit in

D k(|∇(w n-1 + u * )| 2 )∇(w n + u * ) • ∇v dx - Γ i ϕv ds = 0, we obtain D k(|∇(w + u * )| 2 )∇(w + u * ) • ∇v dx - Γ i ϕv ds = 0,
for each v ∈ W . Thus, w is the unique weak solution of the nonlinear problem (4). It follows that

u n = w n + u * → u = w + u * in H 1 -norm as n → ∞.
3. Existence and uniqueness of the Cauchy solution. The existence of a solution of the Cauchy problem (1) depends on the compatibility of the Cauchy data f and φ. They are said to be compatible if the Cauchy problem (1) has a solution. We also say that f and φ are consistent for (1).

Proposition 1. For f ∈ H 1 2 (Γ c ), let M = φ ∈ (H 1 2
00 (Γ c )) ; f and φ are compatible .

Then M is dense in (H

Proof. Let u * ∈ H 1 (D) such that u * | Γ c = f , the Cauchy problem (1) is then reduced to find w ∈ H 1 (D) solution of (13)      ∇•[k(|∇(w + u * )| 2 )∇(w + u * )] = 0, in D, w = 0, on Γ c , k(|∇(w + u * )| 2 )∇(w + u * ) • n = φ, on Γ c .
To prove that M is dense is equivalent to prove that

M = φ ∈ (H 1 2
00 (Γ c )) ; φ and 0 are compatible is dense in (H 1 2 00 (Γ c )) . Otherwise, the Hahn-Banach theorem would guarantee the existence of a functional g ∈ H 1 2 00 (Γ c ) such that g ≡ 0 and φ, g = 0, ∀φ ∈ M .

Given an arbitrary test function, θ ∈ C ∞ 0 (Γ i ), the following mixed boundary-value problem:

(14)      ∇•[k(|∇(w + u * )| 2 )∇(w + u * )] = 0, in D, w = 0, on Γ c , k(|∇(w + u * )| 2 )∇(w + u * ) • n = θ, on Γ i ,
has a unique solution w ∈ H 1 (D). Now, let v be the unique solution of the following linear mixed boundary-value problem:

(15)      ∇•[k(|∇(w + u * )| 2 )∇v] = 0, in D, v = g, on Γ c , k(|∇(w + u * )| 2 )∇v • n = 0, on Γ i ,
where w is the solution of (14). Multiplying (14) by the function v and (15) by the function w and integrating over D, we get:

D ∇•[k(|∇(w + u * )| 2 )∇(w + u * )]v dx - D ∇•[k(|∇(w + u * )| 2 )∇v]w dx = 0.
Using Green's formula we obtain:

- D k(|∇(w + u * )| 2 )∇u * • ∇v dx + Γ i θv ds + Γ c k(|∇(w + u * )| 2 )∇(w + u * ) • ng ds = 0. We have k(|∇(w + u * )| 2 )∇(w + u * ) • n| Γ c ∈ M , therefore: Γ c k(|∇(w + u * )| 2 )∇(w + u * ) • ng ds = 0.
Taking u * solution of: ( 16) [36, Example 5.1]). It follows that: So v| Γ i = 0, and then v is also solution of:

∇u * = 0, in D, u * = f, on Γ c , which is unique (in the variational sense) in H 1 (D) for f ∈ H 1 2 (Γ c ) (see
Γ i θv ds = 0, ∀ θ ∈ C ∞ 0 (Γ i ).
(17)      ∇•[k(|∇(w + u * )| 2 )∇v] = 0, in D, v = 0, on Γ i , k(|∇(w + u * )| 2 )∇v • n = 0, on Γ i .
Using the unique continuation property for linear elliptic equations [25], we conclude that v ≡ 0 in D, contradicting the choice of g. Proof. Let u 1 and u 2 be two solutions of the Cauchy problem (1) and let U = u 1 -u 2 . Then, as we have done for the uniqueness of solution of the direct problem, we conclude that U is solution of the following linear elliptic Cauchy problem:

(18)          ∇•[M (u 1 , u 2 )∇U ] = 0, in D, U = 0, on Γ c , M (u 1 , u 2 )∇U • n = 0, on Γ c .
Using the condition (d) and the unique continuation property for elliptic equations [25], we deduce that U ≡ 0 in D, and hence u 1 ≡ u 2 .

4. Nash game formulation of the Cauchy problem. The mathematical theory of games was invented by John von Neumann and Oskar Morgenstern (1944). It is formal study of conflict and cooperation designed for the case of multi-objective and multidisciplinary optimization. It substitutes the notion of optimum, irrelevant when more then one criteria is under consideration, by the notion of Nash equilibrium which is a fundamental concept in the theory of games. Applications of game theory started with economic problems. Nowadays, it has been extensively studied and applied to a wide range of disciplines. For important techniques and results in game theory and its connections to PDEs the reader is directed to [17,20]. We refer to [34,42] for a general introduction and proof of convergence of computational methods for Nash equilibria, and to [3,13] for a study of alternating algorithms which are closely linked to our present approach.

In [18,19], the authors introduced a game theory based algorithm was introduced for solving the Cauchy problem for a linear operator and convergence results were proved. Here, we extend it to solve our nonlinear Cauchy problem. The Cauchy problem is formulated as a two-player game. The first player is given the known Dirichlet data f and uses the Neumann condition prescribed over the inaccessible Γ i part of the boundary ∂D as strategy variable. The second player is given the known Neumann data φ, and plays with the Dirichlet condition prescribed over the inaccessible boundary. Following the work in [18], for all (η, τ ) ∈ (H

1 2 00 (Γ i )) × H 1 2 (Γ i ), we define      J 1 (η, τ ) = 1 2 (k(|∇u 1 (η)| 2 )∇u 1 (η) • n -φ 2 (H 1 2 00 (Γ c )) + 1 2 u 1 (η) -u 2 (τ ) 2 H 1 2 (Γ i ) , J 2 (η, τ ) = 1 2 u 2 (τ ) -f 2 H 1 2 (Γ c ) + 1 2 u 1 (η) -u 2 (τ ) 2 H 1 2 (Γ i )
,

where u 1 (η) and u 2 (τ ) are the solutions to:

(SP 1 )      ∇•[k(|∇u 1 | 2 )∇u 1 ] = 0, in D, u 1 = f, on Γ c , k(|∇u 1 | 2 )∇u 1 • n = η, on Γ i , and
(SP 2 )      ∇•[k(|∇u 2 | 2 )∇u 2 ] = 0, in D, u 2 = τ, on Γ i , k(|∇u 1 | 2 )∇u 2 • n = φ, on Γ c .
The two players solve in parallel the associated boundary-value problems (SP 1 ) and (SP 2 ). We seek to find a couple (η c ; τ c ), called "Pareto optimal Nash equilibrium" in the game theory vocabulary.

Definition 4.1. A pair (η N , τ N ) ∈ (H 1 2 00 (Γ i )) × H 1 2 (Γ i ) is called:
• Nash equilibrium for the two-player game involving the costs functionals J 1 and J 2 if:

J 1 (η N , τ N ) ≤ J 1 (η, τ N ), ∀ η ∈ (H 1 2 00 (Γ i )) , J 2 (η N , τ N ) ≤ J 1 (η N , τ ), ∀ τ ∈ H 1 2 (Γ i ).
• Pareto optimal Nash equilibrium for the two-player game involving the cost functionals J 1 and J 2 if there does not exist another Nash equilibrium (η n , τ n )

∈ (H 1 2 00 (Γ i )) × H 1 2 (Γ i ) such that: J 1 (η n , τ n ) ≤ J 1 (η N , τ N ) and J 2 (η n , τ n ) ≤ J 1 (η N , τ N ).
Proposition 3. When the Cauchy problem (1) has a solution u, the pair

(η c , τ c ) = (k(|∇u| 2 )∇u • n| Γ i , u| Γ i )
is the unique Pareto optimal Nash equilibrium for the twoplayer game involving the costs functionals J 1 and J 2 .

Proof. First, we prove that (η c , τ c ) is a Nash equilibrium. In fact, we have:

J 1 (η c , τ c ) = J 2 (η c , τ c ) = 0.
From the positivity of J 1 and J 2 , we obtain:

J 1 (η c , τ c ) ≤ J 1 (η, τ c ), ∀ η ∈ (H 1 2 00 (Γ i )) and J 2 (η c , τ c ) ≤ J 1 (η c , τ ), ∀ τ ∈ H 1 2 (Γ i ). Then (η c , τ c ) is a Nash equilibrium. Let (η * , τ * ) ∈ (H 1 2 00 (Γ i )) × H 1 2
(Γ i ) be another Nash equilibrium, then we have:

J 1 (η c , τ c ) < J 1 (η * , τ * ) and J 2 (η c , τ c ) < J 1 (η * , τ * ).
Otherwise, we obtain J 1 (η c , τ c ) = J 1 (η * , τ * ) = 0 and J 2 (η c , τ c ) = J 2 (η * , τ * ) = 0, from the uniqueness of solution of the the Cauchy problem, we obtain (η c , τ c ) = (η * , τ * ). So (η c , τ c ) is the unique Pareto optimal Nash equilibrium. To compute the Nash equilibrium, we use the classical alternating minimization algorithm with relaxation [42], also referred to as the inertial Nash equilibration process, by means of the following iterative process: Algorithm 1:

• Set k = 0 and choose an initial guess S (0) = (η (0) ; τ (0) ).

• Step 1: Compute η (k) solution of min η J 1 (η, τ (k) ). • Step 2: Compute (in parallel) τ (k) solution of min τ J 2 (η (k) , τ ). • Step 3: Set S (k+1) = (η (k+1) ; τ (k+1) ) = t(η (k) ; τ (k) ) + (1 -t)(η (k) ; τ (k) ), 0 < t < 1. • If S (k+1) -S (k) ≤ , stop. Otherwise k = k + 1, go to Step 1.

Remark 2.

i) For the case k(r) ≡ 1, we have a Cauchy problem for the Laplace operator. Nash game formulation of this problem was studied in [18,19]. The costs functionals are elliptic and convex. The authors proved that there exists always a unique Nash equilibrium, which turns out to be the reconstructed data when the Cauchy problem has a solution. They also proved that the completion algorithm is stable with respect to noise. The numerical procedure described in Algorithm 1 is proven to be convergent in the linear case [42]. ii) For the nonlinear case, i.e., k(r) ≡ 1, the partial functionals η → J 1 (η, τ ) and τ → J 2 (η, τ ) are coercive but not convex. The convergence of the numerical iterative procedure for the nonconvex case is more complicated and can be found in [13]. In fact, the functional J 1 +J 2 must clearly satisfy the conditions in [13] which hold in our case except the Kurdyka-Lojasiewicz inequality.

The gradients can be efficiently computed by means of the adjoint state method (wherever it is involved, the inner product in the space (H 1 2 00 (Γ i )) is represented using a classical integral notation). They are given in the following proposition whose proof is given in Appendix A.

Proposition 4. We have the following two partial derivatives:

∂J 1 ∂η (η, τ )h = - Γ i λ 1 h ds, ∀h ∈ (H 1 2 00 (Γ i )) ,
where λ 1 ∈ W solves the adjoint problem

D B(∇u 1 )∇σ • ∇λ 1 dx = Γ c (k(|∇u 1 | 2 )∇u 1 • n -φ)(B(∇u 1 )∇σ • n) ds + Γ i (u 1 -u 2 )σ ds, ∀σ ∈ W. and ∂J 2 ∂τ (η, τ )ξ = Γ i (B(∇u 2 )∇λ 2 • n + u 2 -u 1 )ξ ds, ∀ξ ∈ H 1 2 (Γ i ),
where λ 2 solves the adjoint problem where the matrix-valued function B(•) is given in (d), u 1 and u 2 are the solutions to (SP 1 ) and (SP 2 ).

     ∇ • (B(∇u 2 )∇λ 2 ) = 0 in Ω, λ 2 = 0 on Γ i , (B(∇u 2 )∇λ 2 ) • n = f -u 2 on Γ c ,
5. Numerical results. In this work, all the PDEs are numerically solved using the free finite-element software FreeFem++ [24]. We have used a structured mesh adapted to the image where every node of the mesh corresponds to a pixel in the image.

5.1. Application to image inpainting. We consider an image f defined on a domain Ω, usually a rectangular domain with piecewise smooth boundary ∂Ω. Let D ⊂ Ω be a damaged subregion in Ω where the information is not available. A PDE inpainting model consists in filling-in image information in the region D of incomplete information based upon the image information available outside D, i.e., in Ω\D.

In the last few years, a number of PDE-based inpainting models have appeared. Some models are written down directly, based on some mathematical understandings on the properties of diffusion PDEs (the heat equation, the Cahn-Hilliard equation, etc.). Others are formulated as a minimization for an appropriately regular functional. In [10], Chan and Shen exploited the total variation based image denoising model of Rudin, Osher, and Fatemi [39] for the inpainting purpose. They proposed in [11] the Curvature Driven Diffusion (CDD) approach to the problem of image inpainting. The diffusion coefficient depends on the curvature of the isophotes. Another inpainting approach is based on the complex Ginzburg-Landau equation [16].

In our work, we consider the following function

k α (r) = 1 √ r + 2 + α, ∀r ∈ R + ,
where α > 0 and > 0 are small parameters. The function k α (•) is a regularized version of the diffusion function used in total variation (TV) [10] which consists of joining the same level lines on both sides of D and minimizing its lengths. In the TV model, i.e., without the regularizing parameter α, the solution must be in W 1,1 (D). Existence and uniqueness of the solution of the mixed direct problem were studied in [41]. Unfortunately, until now, nothing can be said about Cauchy problems in W 1,1 (D) space. In order to work in the Hilbert space H 1 (D), in which the Cauchy problem has been studied theoretically and numerically, we have added the constant α in the diffusion function k α (•). This choice is also compatible with the natural conditions (2) that guarantee the existence and uniqueness of solutions of the direct and Cauchy problems. We distinguish two typical cases of image inpainting problems. The first, is when the region D is totally inside Ω (i.e., ∂D ∩ ∂Ω = ∅), see Fig. 1(a). The information available near the boundary ∂D can be used as Dirichlet boundary conditions for the partial differential equation that propagates the information inside D. We then obtain the following problem:

(19) ∇•[k α (|∇u| 2 )∇u] = 0, in D, u = f, on ∂D.
The second, is when ∂D ∩ ∂Ω = ∅, see Fig. boundary into two parts Γ i and Γ c such that

(20) Γ c Γ i = ∂D, Γ i ⊂ ∂Ω, Γc Γi = ∅, Γc ∂Ω = ∅,
where Γc and Γi are the interiors of Γ c and Γ i , respectively. Here, Dirichlet boundary conditions are unknown on Γ i and thus classical methods are not appropriate.

In this work, we introduce a new model that is more relevant for this case in which we ignore the condition on Γ i and we use both Dirichlet and Neumann boundary conditions on Γ c . The later condition can be computed from information available in Ω\D. The main originality of this approach lies on the use of the second Neumann condition on Γ c . Therefore, the inpainting problem can be formulated as the inverse nonlinear boundary-value problem (1)

with k = k α .
In what follows, we present the numerical results obtained using the algorithm described in the previous section. To evaluate the effectiveness of the proposed algorithm, a numerical study is carried out on different images. In all numerical tests, we have chosen images where the missing region D is such that ∂D ∩ ∂Ω = ∅, see Fig. 1(b). In all examples, the damaged regions are marked with red. We also present numerical results which illustrate the stability of our method with respect to noisy data.

We recall that the MSE (Mean Square Error) and PSNR (Peak-Signal-to-Noise Ratio) are used to estimate the quality of an image I 2 with respect to a reference image I 1 , and are defined by the expressions:

M SE = 1 mn m-1 i=0 n-1 j=0 |I 1 (i, j) -I 2 (i, j)| 2 and P SN R = 10 • log 10 M AX I 2 M SE .
Here, m and n are respectively the numbers of lines and columns of the the matrix representing the image, and M AX I 2 is the maximum possible pixel value of the image I 2 .

5.1.1. Linear case. We consider the case k(r) ≡ 1, which was studied in the work proposed in [18,19]. We start by applying this for the image inpainting problem where the missing information part does not contain edges. The example in Fig. 3 deals with the application of this model to color images with a 1% Gaussian noise. We recall that color image contains three channels (u R , u G , u B ) representing red, green and blue intensities. The treatment of color images is a natural extension of the treatment of scalar (i.e., gray-scale) images, for which we solve a Cauchy problem for each channel independently. For one channel, we present in Fig. 2 the L 2 -error of Nash strategies, S (k) -S (k-1) and L 2 -error of missing Dirichlet data τ k -u| Γ i as functions of overall Nash iterations. We remark that both errors behave well and decrease with the number of iterations.

Nonlinear case.

Because the linear model has a strong regularization effect, it can not be used to inpaint images with edges since it can not distinguish between edges (discontinuities) and noise. To take care of edges in non-smooth images, we use the nonlinear model. We give some numerical results where we have chosen an image of a broken line. We present two cases; a rectangular (Fig. 4) and an L-shaped (Fig. 7) damaged domain D. We have tested different values of α in our algorithm. In the first case, we illustrate in Fig. 5 the results obtained using our algorithm for α ∈ 1, 10 -2 , 10 -6 . In Fig. 6, we present the PSNR and MSE of the reconstructed image as functions of different values of α. It can be seen that we obtain a better solution by choosing α small enough. first restore the known portion of the noisy image and then use the restored image as an input in the inpainting step. The results using these two approaches are similar. In Fig. 8 we present the reconstructed image for 1%, 3% and 5% Gaussian noise using the first approach. The PSNR is 31.06, 30.2 and 27.16, respectively. The MSE is 7.9 × 10 -4 , 9.6 × 10 -4 , and 1.9 × 10 -3 , respectively.

Data completion over Γ

i .
To demonstrate the efficiency of the proposed method from data completion viewpoint, we make a comparison with the Kozlov algorithm [30] which is an alternative method for solving the Cauchy problem. It was introduced for solving ill-posed problem and was recently used for solving a similar nonlinear Cauchy problem in glaciology [4]. The basic idea is to reduce the ill-posed problem to a sequence of well-posed mixed boundary-value problems and consists of the following steps:

     ∇•[k α (|∇v 0 | 2 )∇v 0 ] = 0, in D, v 0 = µ 0 , on Γ i , k α (|∇v 0 | 2 )∇v 0 • n = φ, on Γ c , to obtain ψ 0 = k α (|∇v 0 | 2 )∇v 0 • n.
Step 2. For an iteration j > 1, solve the following well-posed mixed boundaryvalue problem:

(22)      ∇•[k α (|∇u j | 2 )∇u j ] = 0, in D, u j = f, on Γ c , k α (|∇u j | 2 )∇u j • n = ψ j-1 , on Γ i .
Set µ j = u j | Γ i and solve the following well-posed mixed boundary-value problem:

(23)      ∇•[k α (|∇v j | 2 )∇v j ] = 0, in D, v j = µ j , on Γ i , k α (|∇v j | 2 )∇v j • n = φ, on Γ c .
Step 3. Set ψ j = k α (|∇v j | 2 )∇v j • n and repeat Step 2 until a prescribed stopping criterion is satisfied.

More detailed results concerning the convergence of this algorithm in the nonlinear case and the stopping criterion can be found in [4]. We illustrate in Fig. 9 the missing data τ N and η N at convergence of the Algorithm 1 dedicated to the computation of the Nash equilibrium for different values of α. In Fig. 10 we give the results obtained using Kozlov algorithm. For both algorithms, the numerical Dirichlet solution remains good (see the left-hand plots of Fig. 9 and Fig. 10) for α = 10 -8 . On the other hand, concerning the numerical Neumann solution (see the right-hand plots of Fig. 9 and Fig. 10), it can be seen that we have two peaks which correspond to edges in the image.

We can say that the Nash game algorithm gives better approximation specially for the numerical Neumann solution. Indeed, the height of the peak in the Neumann solution using to game method is larger then the peak in the Neumann solution using Kozlov algorithm.

To see the sensitivity to the noise, Dirichlet and Neumann numerical solutions using Nash game and Kozlov algorithms are shown in Fig. 11 This work can be extended in many different ways. For instance, one might consider using our approach for other inpainting models such as the one based on Ginzberg-Landau equation [16]. One can also theoretically study the convergence of the numerical method used for solving the Nash-game formulation.

Appendix A. Proof of Proposition 4. To evaluate the gradient of J 1 and J 2 , we must calculate the Fréchet derivatives of the function Θ 1 : η -→ u 1 (η) which is solution of (SP 1 ) and Θ 2 : τ -→ u 2 (τ ) which is solution of (SP 2 ). For (h, ξ) ∈ (H

1 2 00 (Γ i )) ×H 1 2 (Γ i )
, we denote by u 1 (h) and u 2 (ξ) the derivatives dΘ 1 dη ψ and dΘ 2 dτ ξ,

respectively. Then we have:

Lemma A.1. [4]
The Fréchet derivatives u 1 (h) and u 2 (ξ) are solutions of

(24)      ∇ • (B(∇u 1 )∇u 1 ) = 0 in Ω, u 1 = 0 on Γ c , (B(∇u 1 )∇u 1 ) • n = h on Γ i , and (25)
     ∇ • (B(∇u 2 )∇u 2 ) = 0 in Ω, u 2 = ξ on Γ i , (B(∇u 2 )∇u 2 ) • n = 0 on Γ c ,

respectively.

Inverse Problems and Imaging Volume 9, No. 3 (2015), 853-874 Proof. (of Proposition 4) For (η, τ ) ∈ (H

1 2 00 (Γ i )) × H 1 2 (Γ i ), we define the following Lagrangian L(η, τ, u 1 , u 2 , λ 1 , λ 2 ) = 1 2 (k(|∇u 1 (η)| 2 )∇u 1 (η) • n -φ 2 (H 1 2 00 (Γ c )) + 1 2 u 2 (τ ) -f 2 H 1 2 (Γ c ) + 1 2 u 1 (η) -u 2 (τ ) 2 H 1 2 (Γ i ) + Ω k(|∇u 1 | 2 )∇u 1 • ∇λ 1 dx - Γ i ηλ 1 ds + Ω k(|∇u 2 | 2 )∇u 2 • ∇λ 2 dx - Γ c φλ 2 ds, where (u 1 , u 2 , λ 1 , λ 2 ) ∈ H 1 (D) × H 1 (D) × W × W 1 and W 1 = {u ∈ H 1 (D); u| Γ i = 0}.
We have

∂J 1 ∂η (η, τ )h = Γ c (k(|∇u 1 | 2 )∇u 1 • n -φ)(B(∇u 1 )∇u 1 (h)) • n ds (26) + Γ i (u 1 -u 2 )u 1 (h) ds ∀h ∈ (H
+ Γ i (u 1 -u 2 )σ ds + Ω (B(∇u 1 )∇σ) • ∇λ 1 dx = 0, ∀ σ ∈ H 1 (D).
Using the weak formulation associated to u 1 (h), ( 26) and (27)

for σ = u 1 (h), we obtain ∂J 1 ∂η (η, τ )h = - Γ i λ 1 h ds, ∀h ∈ (H 1 2 00 (Γ i )) ,
where λ 1 solves the adjoint problem

D B(∇u 1 )∇σ • ∇λ 1 dx = Γ c (k(|∇u 1 | 2 )∇u 1 • n -φ)(B(∇u 1 )∇σ • n) ds + Γ i (u 1 -u 2 )σ ds, ∀σ ∈ W = {u ∈ H 1 (D); u| Γ c = 0}.
We now compute the other derivative. We have

(28) ∂J 2 ∂τ (η, τ )ξ = Γ c (u 2 -f )u 2 ds - Γ i (u 1 -τ )ξ ds, ∀ξ ∈ H 1 2 (Γ i ),
and Therefore, we get ( 30)

∂L ∂u 2 σ = Γ c (u 2 -f )σ ds + Ω (B(∇u 2 )∇σ) • ∇λ 2 dx = 0, ∀ σ ∈ W 1 . (29
Ω ∇ • (B(∇u 2 )∇λ 2 )σ dx + Γ c (B(∇u 2 )∇λ 2 ) • n -(f -u 2 ))σ ds.
From (30) we deduce the following problem (31)

     ∇.(B(∇u 2 )∇λ 2 ) = 0, in Ω, λ 2 = 0, on Γ i , (B(∇u 2 )∇λ 2 ) • n = f -u 2 , on Γ c .
The weak formulation associated to (31) is

Ω (B(∇u 2 )∇λ 2 • ∇v) dx = Γ c (f -u 2 )v ds + Γ i (B(∇u 2 )∇λ 2 • n)v ds, ∀v ∈ H 1 (D). From Lemma A.1, we have u 2 | Γ i = ξ. Taking v = u 2 in the above, we obtain Ω (B(∇u 2 )∇λ 2 • ∇u 2 ) dx = Γ c (f -u 2 )u 2 ds + Γ i (B(∇u 2 )∇λ 2 • n)ξ ds.
Using the weak formulation associated to u 2 (ξ) and the boundary condition λ 2 | Γ i = 0, we then obtain

Ω (B(∇u 2 )∇λ 2 • ∇u 2 ) dx = 0.

Inverse Problems and Imaging

Volume 9, No. 3 (2015), 853-874

Therefore,

Γ c (u 2 -f )u 2 ds = Γ i (B(∇u 2 )∇λ 2 • n)ξ ds, and ∂J 2 ∂τ (η, τ )ξ = Γ i (B(∇u 2 )∇λ 2 • n + u 2 -u 1 )ξ ds, ∀ξ ∈ H 1 2 (Γ i ),
where λ 2 solves the adjoint problem (31).

"... partial differential equations are the basis of all physical theorems."

Bernhard Riemann (1826-1866).

Chapter 3

Weighted harmonic and complex Ginzburg-Landau equations for gray value image inpainting

The aim of this chapter is to provide some practical and efficient strategies for the choice of the regularization parameters in image inpainting and restoration. It summarizes the work presented in [12].

Summary

Most of variational approaches to image inpainting and restoration problems assume that image is to be reconstructed by considering the following minimization problem:

αR(u) + 1 2 Ω λ D (u -f ) 2 dx, (3.1)
where α and λ D = λ 0 χ Ω\D are regularization parameters. In general, the optimization problem (3.1) gives a better representation of the solution if the regularization parameters α and λ 0 are suitably chosen. However, various methods use uniform parameters chosen, mostly, in an empirical way. This choice may not be reliable in all practical situations. For that reason, the development of efficient strategies for parameters regularization choice in image processing was a vast field to discuss and different methods have been proposed in various frameworks. Images are generally composed of multiple objects at different scales. This calls forth desirable "localized" values of α(x) and λ 0 (x) chosen adaptively depending on the different features and scales in the image. In this work, we propose the following linear equation:

-∇.(α(x)∇u α ) + λ D (u α -f ) = 0, in Ω, ∂ n u α = 0, on ∂Ω. (3.2)
Then, we consider a posteriori method for the choice of these regularization parameters. The selection is performed at the discrete level in the framework of the finite element method with the aim to obtain in the restored region the fine features of the initial image. First, we use a metric mesh-adaptation allowing both the refinement and the coarsening the grid in order to

Weighted Harmonic and Ginzburg Landau equations fit the geometry of the solution. Thus, the number of degrees of freedom in the homogeneous areas decreases which makes the method considerably fast. Second, we select adaptively and locally the diffusion coefficient α by decreasing it near edges in order to improve the visual quality of image and by keeping it large homogeneous regions. The selection is performed with the help of suitable and posteriori error indicators which are known to be equivalent to the H 1 -norm of the finite element error, however, they also encode some priory information about large gradients (edges) of the solution. Eventually, we perform a mathematical analysis and we draw connections between the proposed strategy and the Mumford-Shah functional using Γ-convergence tools and based on the works presented in [11,25,28]).

In image inpainting, there is a strong connection between noise removal in the available part of the input image and the inpainting process. To overcome this sensitivity to noise, we then modify the adaptive strategy in the non damaged regions in order to improve the fitting to the data term which allows us to handle simultaneously the inpainting task and the denoising of the available part of the input image. The two steps are coherently integrated by controlling of both parameters α and λ D . Afterward, we extend the adaptive choice of α from the linear diffusion to the complex-Ginzburg-Landau equation [START_REF] Grossauer | Using the complex Ginzburg-Landau equation for digital inpainting in 2D and 3D. Scale Space Methods in Computer Vision[END_REF].

Finally, we present several numerical simulations to test the efficiency of the proposed approach. We make a comparison with the nonlinear anisotropic diffusion method [START_REF] Weickert | Theoretical foundations of anisotropic diffusion in image processing[END_REF], mostly used in restoration problems. The obtained results prove that, although the model is linear, the tactical choice for spatially varying regularization parameters allows us to fit the geometric specificities of the image and to recover sharp edges. Abstract. We consider two second-order variational models in the image inpainting problems.

The aim is to obtain in the restored region some fine features of the initial image, e.g. corners, edges, . . . . The first model is a linear weighted harmonic method well suited for binary images and the second one is its extension to the complex Ginzburg-Landau equation for the inpainting of multi-gray level images. The approach that we introduce consists of constructing a family of regularized functionals and to select locally and adaptively the regularization parameters in order to capture fine geometric features of the image. The parameters selection is performed, at the discrete level, with a posteriori error indicators in the framework of the finite element method. We perform the mathematical analysis of the proposed models and show that they allows us to reconstruct accurately the edges and the corners. Finally, in order to make some comparisons with well established models, we consider the nonlinear anisotropic diffusion and we present several numerical simulations to test the efficiency of the proposed approach.

Introduction

Image inpainting (or disocclusion) refers to restoring a damaged image with missing information. This type of image processing is very important and has many applications in various fields (painted canvas, movies restoration, augmented reality, . . . ). In fact, many images are often scratched and damaged, and the goal in the inpainting problems is to restore deteriorated or missing parts, so that a viewer cannot distinguish them from the rest. Various mathematical and heuristic techniques were considered to address this problem, such as statistical methods [23], mathematical programing and computational geometry methods [31], we refer to the article [11] and the references therein where an exhaustive review is given for this problem and for the various approaches developed to solve it. In this article we will be concerned by the Partial Differential Equations (PDE) approach which belongs to the class of the widely used methods ( [6,12,19,20]). Let Ω ⊂ R d (d = 2, 3), denotes the entire domain of a given image f , the basic idea in the PDE approach, is to fill-in the damaged region D ⊂ Ω, where the pixels of f are altered or lost, by an interpolation from the available part in Ω\D. Usually, the PDE-based models are obtained from the mathematical knowledge of the properties of some differential operators, and aim to fulfill some a priori expectations and assumptions on the final solution. The diffusion operators are the mostly used to this end (e.g. the heat equation, the Cahn-Hilliard equation,... [10,12,15,20,26]). Usually such 2000 Mathematics Subject Classification. 65M32, 65M50, 65M22, 94A08. where the image f is given in Ω\D and n is a Gaussian noise. R(u) denotes the regularizing term, mostly a semi-norm of a functional space fixed a priori to enforce some expectations on the solution (e.g. a Sobolev space H s , Bounded Variations functions space BV, . . . ) and u is the image to be reconstructed. The unconstrained formulation of (1) reads:

(2)

αR(u) + 1 2 Ω λ D (u -f ) 2 dx,
where α is a regularization parameter and λ D = λ 0 χ Ω\D for λ 0 0, a penalization factor, and χ Ω\D is the indicator function of the sub-domain Ω \ D. These two parameters α and λ 0 are chosen in order to balance the regularization term R(u) and the data fitting term.

Various methods use uniform parameters α and λ, chosen in general empirically or within the regularization theory, e.g. with Morozov's criterium when the magnitude of the noise is given. In many applications, this choice is not reliable and may produces the loss of some relevant features of the image such as the edges (see Fig: 1). Therefore, based on the importance of the scale-space representation of the Figure 1. Harmonic inpainting (T. Chan and J. Shen [19]) image, spatially varying choices of the parameter α were proposed in the literature. We mention as an example the variant of the total variation (TV) functional, considered by D. Strong and T. Chan [33] which results in a multi-scale strategy with a uniform α updated at each scale [4]. Others strategies to choose such parameters are also developed within the statistical approach or using some a priori PDE [29] for the denoising problem. Note that the topological gradient method leads implicitly to such a choice by allowing the modification of the diffusion coefficients [5,6].

We consider in this article a novel approach which consists of an adaptive method for the choice of such spatially varying regularization parameters. The method is well-suited for images with few textures and was successfully applied to the segmentation problems [8]. Loosely speaking, we start with a simple model (e.g. linear diffusion with a variable coefficient), then iteratively, an adaptive selection of the parameters based on some local information on the gradient magnitude is performed. The gradient information are available at the discrete level from the computed solution, thus the process is completely an a posteriori method without any reference to the continuous solution of (2). This amount to change dynamically the reconstruction model in order to capture accurately the fine geometric structures of the image. This approach was introduced by Hecht and Belhachmi in [9] for the optic flow estimation problem, where it was demonstrated to have several attractive features such as: the efficiency (e.g. the cost of computations, best representation of the solution,. . . ) as well as a good edge-preserving property. Moreover, it was proven in [8] that it allows one to approximate, in the Γ-convergence sense [14], the Mumford-Shah functional (see [8,17,18]) although formally the continuous model remains linear (with respect to the principale variable).

The article is organized as follows: In Section 2, we introduce a weighted regularizing functional to obtain the suitable modified version of the harmonic model and we establish its properties. In Section 3, we introduce the discrete framework of the method and we make a selective diffusion, controlled by a suitable error indicators. Using ( [8,17,18])), we perform the Γ-convergence analysis of the method. We also modify the adaptive strategy in the non damaged regions in order to improve the fitting to the data term which allows us to handle simultaneously the inpainting task and the denoising of the available part of the input image. In section 4, we extend such an approach from second order linear diffusion to the complex-Ginzburg-Landau energy which is known, at least numerically, to enhance the contrast in inpainting problems and is well suited for multi-gray level images [3,24]. We present several numerical simulations to show the performances of the method for the considered models. We also make some comparisons with the nonlinear anisotropic diffusion method which belongs to the well established techniques in the image inpainting [35].

Weighted harmonic inpainting

We assume that the domain Ω is partitioned into a disjoint finite number of subdomains Ω i , i = 1, . . . , I. and we consider a function α which is scalar, piecewise constant in Ω and such that α = α i in Ω i , i = 1, . . . , I.

We denote by α m = min 1≤i≤I α i , α M = max 1≤i≤I α i , and we assume that α m > 0. We consider the following linear equation:

(3) -∇.(α(x)∇u α ) + λ D (u α -f ) = 0, in Ω, ∂ n u α = 0, on ∂Ω.
Remark 1. It should be emphasized here that the parameter λ 0 is intended to be large enough to penalize the constraint u α = f in Ω \ D and ( 3) is equivalent to the the following transmission problem:

(4)

           -∇.(α(x)∇u α ) + λ 0 (u α -f ) = 0, in Ω\D, ∇.(α(x)∇u α ) = 0, in D, [u α ] = 0, on ∂D, [α∇u α • n] = 0, on ∂D, ∂ n u α = 0, on ∂Ω,
where [•] denotes the jump across ∂D.

We define the subspace V = {u ∈ H 1 (Ω); D u dx = 0}. Therefore, under the previous assumptions on the function α, we have:

Proposition 1. Let f ∈ L 2 (Ω)
, then the problem (3) admits a unique weak solution u α in V .

Proof. Equation ( 3) is the optimality condition of the following minimization problem: (5) min One may check directly that F α is convex and weakly lower semi-continuous in H 1 (Ω). For u ∈ V , we have:

v∈V {F α (v) = Ω α(x)|∇v| 2 dx + Ω λ D (v -f ) 2 dx}.
F α (u) ≥ α m Ω\D |∇u| 2 dx + λ 0 Ω\D (u -f ) 2 dx + α m D |∇u| 2 dx.
Using the previous inequality and applying the Poincaré-Wirtinger inequality in D, we get:

F α (u) ≥ c u 2 H 1 (Ω)
, where the constant c is dependent on α m , λ 0 and the geometry of D. which implies that F α is coercive. Thus, the functional F α admits a minimizer u α ∈ V . The uniqueness is guaranteed by the strict convexity of F α .

The weak formulation of (3) reads: (6) find u α ∈ V, such that:

a α (u α , v) = l(v), ∀v ∈ V, where (7)
     a α (u, v) = Ω α(x)∇u • ∇vdx + Ω λ D uvdx, l(v) = Ω λ D f vdx.
The equivalence of the problems ( 6) and ( 5) follows by standard arguments. Note that if Ω is Lipschitz-continuous, f ∈ L 2 (Ω) and λ D ∈ L +∞ (Ω), the following regularity result holds [7, Proposition 2.5] Proposition 2. There exists a constant c only depending on the geometry of Ω, such that a weak solution u α of problem (6) belongs to H s+1 (Ω), for all real numbers s < s 0 , where s 0 is given by

s 0 = min 1 2 , c| log(1 - α m α M )| .
Remark 2. This result reminds us that even non-smooth the solution of ( 6) is H 1 (Ω) and therefore admits no jump inside Ω. Nevertheless, our approach consists in decreasing the diffusion coefficient α in high gradient zones (formally to zero) encouraging possible jumps in these areas.

2.1. Discrete problem and adaptivity. 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,
and the following notations: for

u h , v h ∈ X h : (8)
     a α,h (u h , v h ) = Ω α h (x)∇u h • ∇v h dx + Ω λ D u h v h dx, l h (v) = Ω λ D f h v h dx,
where f h is a finite element approximation of f associated with T h . The discrete problem leads to:

(9) find u α,h ∈ X h , such that: a α,h (u α,h , v h ) = l h (v h ), ∀v h ∈ X h .
Proposition 3. There exists a unique solution u α,h in X h of the discrete problem (9).

Proof. The proof can be carried out by applying the Lax-Milgram Lemma. Furthermore, we have following finite element error:

u α -u α,h V ≈ O(h).
Remark 3. We do not impose any compatibility of the mesh with the"partition" D ∪ Ω \ D. We are given a regular mesh over Ω similarly to the fictitious domain methods.

Adaptive local choice of α.

For an element K ∈ T h , we denote by E K the set of its edges not contained in the boundary ∂Ω. The union of all E K , K ∈ T h is denoted by E h . With each edge e ∈ E h , we associate a unit vector n e normal to e and we denote by [φ] e the jump of the piecewise continuous function φ across e in the direction n e . For each K ∈ T h , we denote by h K the diameter of K and we denote by h e the length of e, e ∈ E K and f h a finite element approximation of f . We define the residual error indicator as follows: for each element K ∈ T h , we set:

η K = α -1 2 K h K ||λ 1 2 D (u α,h -f h )+α h ∆u α,h || L 2 (K) + 1 2 e∈E K α -t 2 e h 1 2 e ||[α∇u α,h •n e ] e || L 2 (e) ,
where α e = max(α K 1 , α K 2 ), K 1 and K 2 being the two elements adjacent to e.

On the triangulation T h , we compute the solution u α,h of problem ( 9) and the corresponding error indicator which is well know to be equivalent to the H 1 -norm of the finite element error (see [8] for details) and allows mostly mesh adaptation. η K gives the error distribution of the computation of u α,h , and includes information about edges in the following term:

(10)

1 2 e∈E K α -1 2 e h 1 2
e ||[α∇u α,h • n e ] e || L 2 (e) .

In fact, the edges in the image are characterized by the brightness changes (large gradients). Therefore the quantity (10) acts as a measure locating regions of edges and will be used next to control the parameter α.

Remark 4. The gradient can represent the change in gray level and his magnitude provides information about the strength of the edge. Since all error indicators are (mainly) equivalent [34], we may change the error indicator η K by the following local energy:

(11) η K = α 1 2 K h 1 2 K ||∇u α,h || L 2 (K)
, which might be well suited in the adaptation steps and behaves like the residual error indicator.

Adaptive strategy. We control the diffusion process by following the adaptive algorithm: Given the initial grid T 0 h in Ω, we: (1) Compute u α 0 ,h solution of the problem (3) on T 0 h with a large constant α = α 0 .

(2) We build an adapted isotropic mesh T 1 h (in the sense of the finite element method, i.e. with respect to the parameter h) with the metric error indicator ( [8]).

(3) We perform an automatic local choice of α(x) on T 1 h to obtain a new function α 1 (x) in D.

α k+1 K = max      α k K 1 + κ * η K ||η|| ∞ -0.1 + , α thr      ,
where α trh is a threshold parameter and κ is a coefficient chosen to control the rate of decrease of α, (u + ) = max(u, 0). Here η is the piecewise-constant function such that η| K = η k , ∀K ∈ T 1 h . The formulae (12) means that in the regions of high gradients, one decreases the values of α. Actually, if the error indicator deviates more then 10% from its mean value, then there is a large error which indicate that the element contains a part of the singular set of u. Therefore, decreasing α (nearly as a geometric sequence with the iteration number) produces an edge location.

The adaptive algorithm consists of two steps. First, given α, we solve a linear equation ( 3) and build an adapted isotropic mesh T h . The adapted mesh is obtained by coarsening the initial grid in the homogeneous regions and by refinement to obtain smaller elements 'close' to the jump set of u. Second, we update the value of α in every element K of T h in accordance with the formulae (12).

Γ-convergence analysis of the adaptive algorithm.

A Γ-convergence study of the adaptive strategy is performed in [8] for the optical flow problems. Analyzing this strategy, the authors proved that it is equivalent to the adaptive algorithm introduced for denoising, by Chambolle-Dal Maso in [18] and Chambolle-Bourdin in [17] where a similar method for the numerical discrete approximation of the Mumford-Shah energy was proposed. They proved that this method, based on finite element discretization and adaptive mesh strategy, is a good approximation in the Γ-convergence sense [14] of the Mumford-Shah energy. We briefly recall the results and the numerical approximation of this method. For a fixed angle 0 < θ 0 ≤ π/3, a constant c ≥ 6, and for > 0, we set T (Ω) = T (Ω; θ 0 ; c) be the set of all triangulations of Ω whose triangles K have the following characteristics:

i) The length of all three edges of K is between and c.

ii) The three angles of K are greater than or equal to θ 0 . Let V (Ω) the set of all continuous functions u : Ω -→ R such that u is affine on any triangle K of a triangulation T ∈ T (Ω) and for a given u, T (u) ⊂ T (Ω) is the set of all triangulations adapted to the function u, i.e., such that u is piecewise affine on T. They introduce 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 [18] introduced the following minimization problem:

(13) 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 (Ω), +∞, Otherwise.
For going to zero and provided θ 0 is less than some Θ > 0, they proved that the energy G Γ-converges to the Mumford-Shah functional:

G(u) = Ω |∇u| 2 dx + g ∞ H 1 (S u ), u ∈ L 2 (Ω) ∩ GSBV (Ω), +∞, u ∈ L 2 (Ω)\GSBV (Ω),
where H 1 is the 1-dimensional Hausdorff measure and GSBV (Ω) is the generalized special function of bounded variation (see [1]). It follows from the Γ-convergence to G [18, Theorem 2]:

Theorem 2.1. Let (u ) >0 be a family of functions such that u ∈ V (Ω) for all > 0 and

(14) sup >0 G (u ) + ||u || L 2 (Ω) < +∞.
Then there exists u ∈ GSBV (Ω) and a subsequence (u j ) j converging to u, a.e. in Ω, such that:

G(u) ≤ lim inf G j (u j ),
and, if for each , u is a solution of:

(15) min v G (v) + Ω λ D |v -f | 2 dx,
then the limit u solves:

(16) min v G(v) + Ω λ D |v -f | 2 dx,
and u j converges strongly to u.

From convex analysis, we can write:

g(t) = min v∈[0,1] tv + ψ(v),
where ψ is the Legendre-Fenchel transform of g. The minimum is achieved for v = g (t) and therefore, for a given triangulation T , the minimization of G is then equivalent to minimize the following functional:

G (u, v, T ) = K∈T |K ∩ Ω| 1 h 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 . For a fixed u, the minimizer over each v is explicitly given by:

(17) v K = g (h K |∇u| 2 ),
and the optimal u for fixed v solves an elliptic equation. The adaptive algorithm minimizes G with v = α.

Remark 5. .

• Given a function α, the computation of the minimization problem with respect to u is simple and very fast because after each adaptation step, one solves a linear problem with the number of nodes of the mesh which decreases.

• For computational reasons, we have chosen the formulae (12) to update the diffusion function α (which is similar to the choice g(t) = λ min(t, µ) for given constants λ and µ). This allows for stable computations, other possible choices to update α are possible (e.g. in [17], the authors considered a smooth function g(t) = M arctan( t M ), for M > 0 which leads according to the formulae ( 16) to the diffusion function

α K = 1 1 + (h K |∇u|) 2 .
• It may be noticed that the parameter α acts as a "phase field function"

and plays a role similar to the z-field in the Tortorelli-approximation [2] method for the Mumford-Shah energy. However, the edges obtained in our case seems sharper and their "thickness" is controlled by the refinement strategy. This behavior of α is shown in Fig: 9 in the numerical examples.

2.4. The inpainting with the image restoration. Numerical evidences show that there is a strong connection between the noise removal in the available part of the initial image and the inpainting process. To overcome this sensitivity to the noise, a natural idea is to perform simultaneously the denoising in Ω \ D and the fill-in in D. The two steps are coherently integrated in the method. Thus, we replace now the previous constant λ 0 with a spatially varying function λ(x) and we will select locally it values in order to decide wether we should encourage the penalization (by increasing λ) or not. As mentioned in the remark (1), the residual error indicator splits as for the transmission problem:

(18)                   
If K ∩ Ω\D = ∅, we have :

η K = α -1 2 K h K ||λ 1 2 K (u α,h -f h ) + α h ∆u α,h || L 2 (K) + 1 2 e∈E K α -1 2 e h 1 2
e ||[α e ∇u α,h • n e ]|| L 2 (e) , and if K ∩ Ω \ D = ∅, we have :

η K = α -1 2 K h K ||α K ∆u α,h || L 2 (K) + 1 2 e∈E K α -1 2 e h 1 2
e ||[α e ∇u α,h • n e ]|| L 2 (e) .

For an element K ⊂ Ω\D, the error indicator contains a supplementary term, i.e.

(

) E K = α -1 2 K h K ||λ 1 2 K (u α,h -f h ) + α h ∆u α,h || L 2 (K) 19
, where λ K is the constant value of λ in the element K. In this term the parameters are competing and finding a balance is not obvious. Since our purpose is to make a noise filtering we choose to keep α fixed and to adjust only λ. In fact, starting with a large value of α 0 in the previous algorithm will smooth the input image f in Ω\D at the first iterations which produces some undesirable blurring in this part of the domain. Thus, in Ω \ D, we keep such α constant and increase λ to enhance the fidelity term. In D, the process is unchanged. This calls for a slight modification of the previous algorithm. We start with a constant value λ 0 . If K ∩ Ω\D = ∅, then we update λ as follows:

(20)

λ k+1 K = min{λ k K * 1 + κ * E K ||E|| ∞ -0.1 + , λ 0 },
The explanation of this formula is identical to that of α previously, that is if the error indicator is 10% larger than the mean value, the fidelity constraint should be enforced by increasing λ at that location.

Adaptive strategy in the complex Ginzburg-Landau equation

We extend such an approach from the linear diffusion to the complex-Ginzburg-Landau energy. This model was originally developed by Ginzburg & Landau in [27] to describe phase separation and it is given by: ( 21)

-∆u + W (u) 2 = 0,
where u : Ω -→ C, is a small positive parameter and W (u) = (1 -|u| 2 ) 2 . It is the Euler-Lagrange equation associated to the minimization of the following energy:

(22) 1 2 Ω |∇u| 2 dx + Ω W (u) 2 2 dx.
For digital image inpainting purposes, this equation was developed by H. Grossauer and O. Scherzer in [24]. The key advantage of this model is that its solutions are known to produce effects like vortices and shockwaves of the phase when -→ 0 and the solution reveals high contrast in the inpainting domain, which makes it well suited for multiple gray level images. The real Ginzburg-Landau equation ( 21) is appropriate only for two-scale images, while the minima of the potential W are attained in the sphere |u| = 1. For grayscale images, we follow the same methodology of Grossauer and Scherzer in [24]. We rescale the intensity of the input image f (x) to the interval [-1, 1]. Then f is identified with the real part of the complex valued function f re : Ω -→ C. We then define:

(23)    f = f re + if im , where: f re = f 0 (the initial damaged image), f im = 1 -f 2 0 .
By this choice, the complex valued solution u will also have an absolute value equal to 1 but our inpainting f re may contain any value in the interval [-1, 1]. The aim is to minimize the following Ginzburg-Landau energy:

(24)

F (u) = Ω α(x) 2 |∇u| 2 dx + Ω W (u) 2 2 dx + 1 2 Ω λ D (u -f ) 2 , over V = {u ∈ H 1 (Ω, C); D u dx = 0}
and where H 1 (Ω, C) is the space of complex functions equipped with the norm:

||u|| 2 1 = Ω uudx + Ω ∇u • ∇udx.
For the sake of clarity, we omit the dependence for the minimizers of F . Then, u α satisfies the following Euler-Lagrange equation:

(25) -∇.[α(x)∇u α ] + 1 2 u α (|u| 2 α -1) + λ D (u α -f ) = 0, in Ω, α(x)∂ n u α = 0, on ∂Ω.
It is readily checked (see [13] for details) Proposition 4. The functional (24) admits a minimizer u α in V . Moreover, u α is a weak solution of ( 25) and |u α | ≤ 1. Evolution equation and discretization. We consider the associated evolution problem:

(26) ∂u α ∂t -∇.[α(x)∇u α ] + W (u α ) 2 + λ D (u α -f ) = 0, in R + × Ω,
with homogeneous Neumann boundary conditions and the initial time condition u α (t = 0, x) = f (x) ∀x ∈ Ω. We assume without loss of generality that Ω f = 0 and f ≤ 1. The weak solution of ( 26) solves: find u ∈ L 2 (0, T ; V ) such that

Ω ∂u α ∂t φdx+ Ω α(x)∇u α ∇φdx+ Ω W (u α ) 2 φdx+ Ω λ D (u α -f )φdx = 0 ∀φ ∈ H 1 (Ω, C)
Time discretization. We use the linearly implicit Euler scheme: given u n α , find

u n+1 α ∈ H 1 (Ω, C) such that Ω u n+1 α -u n α δt φdx + Ω α h ∇u n+1 α ∇φdx + 1 2 Ω (|u n α | 2 -1)u n+1 α φdx + Ω λ D (u n+1 α -f )φdx = 0 ∀φ ∈ H 1 (Ω, C).

Comparison with the nonlinear diffusion methods

For the sake of completeness, we will make a comparison with the nonlinear diffusion method [35] that we recall now. The earliest model considered is the so-called nonlinear isotropic diffusion by Perona and Malik [30]. The diffusion coefficient was a nondecreasing function g of |∇u| 2 with g(0) = 1, g(s) > 0 and lim s-→+∞ g(s) = 0. The method was extended to the anisotropic case by J. Weickert in [35] who replaced the scalar diffusion function g with a diffusion tensor D depending on |∇u σ | 2 , where u σ is a smoothed version of u (convolution with a smoothing kernel). The diffusion is adjusted according to the directional information contained in the image structure. In our case, the anisotropic model may be written as follows:

(27)

     ∂u ∂t -∇.[D(∇u σ )∇u] + λ D (u -f ) = 0, in R + × Ω, u(x, 0) = f, in Ω, < D(∇u σ )∇u • n >= 0, on R + × ∂Ω.
Following [35], let {v 1 , ..., v d }, be an orthonormal basis of R d (d = 2, 3) such that v 1 ||∇u σ . The matrix D is symmetric, positive definite and {v 1 , ..., v d } represent its eigenvectors with corresponding eigenvalues {Λ 1 , ..., Λ d }. Then we have:

D(∇u σ ) = (v 1 |...|v d ) diag(Λ 1 , ..., Λ d )(v 1 |...|v d ) T = (v 1 |...|v d ) diag(Λ 1 , ..., Λ d )(v 1 |...|v d ) T . ( 28
)
These eigenvalues are chosen to be functions of |∇u σ | in order to obtain a diffusion tensor adapted to the local structure of the image (i.e., homogeneous area or edges). Let g ∈ C ∞ ([0, ∞), (0, 1]) be a Lipschitz continuous scalar function which is represented in [0, +∞) by a convergent power series as follows:

g(s) = ∞ k=0 a k s k ,
and we consider the tensor product J 0 (∇u σ ) = ∇u σ ⊗∇u σ . For the two dimensional case, the matrix:

(29) D(∇u σ ) = g( |J 0 (∇u σ )|) = ∞ k=0 a k ( |J 0 (∇u σ )|) k ,
defines a diffusion tensor with an orthonormal system of eigenvectors (v 1 , v 2 ) where:

v 1 = v 1 1 v 2 1 = ∇u σ |∇u σ | and v 2 = v 2 1 -v 1 1 .
The choice that prevents the diffusion over the edge lead to the matrix D:

D =   | | v 1 v 2 | |   Λ 1 0 0 Λ 2   | | v 1 v 2 | |   T .
Different classes of anisotropic models and diffusivity functions "g" may be used (see [21,22,29,35]). In this article we choose the following function:

(30) g(s) = 1 + s 2 /R 2 ,
where R and are a contrast and a resolution parameters, respectively. The variational formulation of problem ( 27) reads:

(31)    Find u ∈ C((0, T ); V ), such that: Ω ∂u ∂t v dx + a(u; u, v) = l(v), ∀v ∈ V, where (32)
     a(w; u, v) = Ω D(∇w σ )∇u • ∇vdx + Ω λ D uvdx, l(v) = Ω λ D f vdx.
Iterative scheme: To find the solution u of the nonlinear problem (27), we start from an initial guess u 0 and we use the explicit Euler scheme with ∂u ∂t replaced by u n+1u n δt (where δt is the time step). We then obtain the following semi-implicit problem:

(33)    Given u n , find u n+1 ∈ V, such that: Ω u n+1 -u n δt v dx + a(u n , u n+1 , v) = l(v) ∀v ∈ V.
At each iteration, the bilinear form a(u n , u n+1 , v) is symmetric and positive definite and the problem ( 33) is well defined. The proof of the next proposition is a straightforward application of the analysis in [16] (see also Weickert [35]).

Proposition 5.

(i) Let f ∈ L 2 (Ω); then there is a unique solution u(t, x) of ( 27), u ∈ C([0, T ]; V ) ∩ L 2 (0, T ; V ), V stands for the dual of V . (ii) Let (u n ) n , denotes the sequence of solutions of the linearized problems (33).

Then, the sequence (u n ) n converges, in C([0, T ]; L 2 (Ω)) as n → +∞, to the solution u of problem (27).

Numerical examples

All the PDEs considered in this section are solved using the finite element open source software FreeFem++ [25]. In all the examples, the damaged regions are marked with red color rectangles. All the examples are for 2-D images. Examples 1. In the first example, the task is to fill broken edges in a white strip. We display in Fig. 2 the evolution of the meshes for iterations 1, 10 and 20. one can see that the meshes are progressively sparse. The first mesh is T 0 which is a regular grid where every pixel corresponds to a node. The harmonic inpainting without adaptation in Fig. 3 does not achieve any connectedness and produces a smooth solution u in D, blurring the edges. At the contrary, with the adaptive algorithm the edges of the strip are well captured. Examples 2. In the second example, we have chosen a 220 × 250 gray-scale image containing some edges and jumps. The damaged regions are numbered (see Fig. 4). We show in Fig. 5 the results obtained using the total variation (left) and the harmonic models without adaptation (middle), and with the adaptation, using the error indicator η K (right). Note that the total variation is approximated here with + |∇u| 2 ( = 0.001). In Fig. 6, we display the inpainted images using the weighted harmonic model where we adapt with the error indicator η k (left), the Ginzburg-Landau with adaptation (middle) and the anisotropic model (right). We have performed the adaptation with the two residual error indicators of Section 2. In both cases, we initialized the algorithm by a large value of α = α 0 = 50 and we performed 20 iterations for the error indicator η K and 40 iterations for η K . We give in Fig. 9 the two error indicators η K and η K at the final iteration of the algorithm. We note that numerically, the error indicator η K gives better results for these examples. In addition, we emphasize that the adaptive method, both the weighted harmonic and the the Ginzburg-Landau equation, gives visually comparable results to those obtained using the anisotropic nonlinear model. The dissimilarities are only seen by zooming (so at few pixels scale). We present in Fig. 10 and Fig. 11 the evolution of the mesh for some iterations (1,2, 10 and 20) where we used η K as an error indicator. The number of elements decreases very quickly at the first iterations and produces sparse solutions requiring few degrees of freedom in the homogeneous zones. This is shown by the curve in Fig. 12 where we presented, in a semi-log scale, the degrees of freedom (right) and the L 2 -error ER = ||u k αu exact || between the restored and the exact image (left) as a function of the number of iterations.

The numerical experiment in Fig. 13 shows the efficiency of the proposed method for textured image. The aim is to remove the foreground text in the input image. We initialized the algorithm with a large value of α = 50 and we performed 10 iterations of the adaptive algorithm. The text have been successfully removed and the restored image is close to the original one. We display the difference between the original and the restored images at iterations 1 and 10 respectively and we give in in Fig. 14 a zoom caption (the resolution is degraded in the initial image however the "blurring" produced by the adaptive method -α > 0-is discernible at this scale). Examples 3. We present in Fig. 15 the result for the simultaneous image inpainting and the denoising. The input image f is corrupted by a Gaussian noise in the region Ω\D. We started the computation with α = 50 in the entire image domain adaptation allows us to approximate and to capture the corner, and reveals high contrast, which is the key advantage of this model compared to the other ones in this case. We display in Fig 18 the evolution of the mesh at iterations 0, 7 and 30 for the complex Ginzburg-Landau equation. Examples 5. The latest example in Fig 19 deals with the reconstruction of the curvature. the inpainted edge in the missing part tends to be a straight line and it is clear that our adaptive algorithm behaves like the Mumford-Shah model (see [17] for the denoising treatment). This behavior is expected because the preferable edge with the Mumford-Shah model are those which have the shortest length due to the penalization term on the length. This example is an extreme case for second order PDEs methods which fail to capture the curvature contrary to the high order PDEs.

Conclusion and perspectives

In this paper, we have considered an adaptive approach for image inpainting based on a local selection of the different parameters in the models, and on mesh adaptation techniques. We started with the formulation of a linear variational model, and detailed its numerical implementation based on the finite element discretization which approximate in the sense of the Γ-convergence the Mumford-Shah functional. We extended the model to the Ginzburg-Landau equation, more suited for multi-gray level images. In order to make some comparisons, we presented the nonlinear diffusion model which is a standard and high quality method in image processing. The numerical experiments for the various examples presented here demonstrate the efficiency of the adaptive method and tends to confirm that finding fine structures in the reconstructed images is a matter of the diffusion more than a non-linearity in the source term in the PDE. We may say that the multi-scale strategy, based on a rigorous adaptive selection of the diffusion "rate" and location, leads to comparable results that one might expect from the nonlinear PDE considered in this field while presenting an evident advantage from the numerical point of view. Finally, the combination of this approach and the complex Ginzburg-Landau model yields very encouraging results.

The adaptive approach of the article may be applied to other problems in image analysis. Since the anisotropic diffusion remains one of the best methods, a first tentative to improve the adaptive approach is to derive an anisotropic version, which means considering α as a matrix (this is an ongoing work). A second step will be to extend it to the fourth-order PDEs which are more suitable for the inpainting by preserving high curvatures ( [28,32]). Finally applying the adaptive approach to the vectorial setting (e.g. color images) would be a challenging problem and will constitutes a breakthrough in the field.

Multiscale fourth-order models

We analyze the proposed model in the framework of calculus of variations and with Γconvergence tools and we establish its connections with a Mumford-Shah model. We prove that the adaptivity steps applied on the previous fourth-order PDE allows us to approximate a Mumford-Shah like energy coupled with with H -1 -fidelity term. In addition, the model has the advantage of being easy to solve numerically and the computation is fast.

Finally, we recall the Cahn-Hilliard in order to make a comparison with our approach. We present several numerical examples, in image inpanting and restoration problems, which show the good quality in the recovery of low dimensional sets (edges, corners) and curvature in the inpainted zone. We consider a fourth-order variational model, to solve the image inpainting problem, with the emphasis on the recovery of low-dimensional sets (edges, corners) and the curvature of the edges. The model permits also to perform simultaneously the restoration (filtering) of the initial image where this one is available. The multiscale character of the model follows from an adaptive selection of the diffusion parameters which allows us to optimise the regularization effects in the neighborhoods of the small features that we aim to preserve. In addition, since the model is based on the high-order derivatives, it favors naturally the accurate capture of the curvature of the edges, hence to balance the tasks of obtaining long curved edges and the obtention of short edges, tip points and corners. We analyze the method in the framework of the calculus of variations and the Γ-convergence to show that it leads to a convergent algorithm. In particular, we obtain a simple discrete numerical method based on a standard mixed-finite elements with well established approximation properties. We compare the method to the Cahn-Hilliard model for the inpainting and we present several numerical examples to show its performances.

Keywords: Image inpainting; Inverse problems; Regularization procedures; Mixed finite elements.

AMS Subject Classification: 65M32, 65M50, 65M22, 94A08, 65N22, 35G15, 35Q68 1

Introduction

The digital image inpainting started with the works of engineers and computer scientists in the mid-nineties of last century and it refers to the restoration of a scratched or damaged image. In the image processing, this task is very important and has many applications in various fields (painted canvas and movies restoration, augmented reality, . . . ). Let Ω ⊂ R 2 denotes the entire image domain and let D ⊂ Ω be a damaged part of Ω. The basic idea in the image inpainting is to fill-in D, in such a manner that a viewer can not detect the restored part from the rest of the image. Different techniques have been applied to solve such a problem, e.g. some heuristic "copy-paste" procedures 31 , statistical approaches 5,34 , we refer to the review article 13 and the references therein for an exhaustive overview on this problem. Among the various methods, the Partial Differential Equations (PDEs) are now widely used and are proven to be efficient in this fields 12,14,17,26,27,32,36,38,39 . The underlying idea in the PDE-based methods is to interpolate with a differential operator the data given outside of D. A large number of differential operator have been considered, in particular, several second-order diffusion equations have been proposed in the literature 17,27,32 . Such operators have the advantage of being wellestablished theoretically and to lead to various simple numerical methods. However, they may suffer from major drawbacks such as the disconnection of level lines over large distances (violating the connectivity principle), the smooth propagation of the edges into the damaged domain (blurring effect), or the failure in reproducing some high-order features (curvatures, . . . ). Note also that these methods do not work when the damaged region D touches the boundary of Ω. In such a situation, the authors in 36 proposed an inpainting method based on the data completion with a nonlinear Cauchy problem.

The shortcomings of the second-order partial differential equations gave rise to a class of high-order diffusion models which in general perform better, particularly for restoring the curvatures and for matching edges across large distances. In fact, the supplementary information on the level lines directions which may be provided from the use of boundary conditions both on the solution and its gradient enhance the quality of the reconstruction. In addition, in the homogeneous areas of Ω the high-order of the operator damps the noise faster than any second order one.

Overview on higher-order PDE models

Bertalmio et al. 14 pioneered a two dimensional fourth-order PDE which consists in propagating both the gradient direction (geometry) and the gray-scale values of the image inside the region to be filled-in by solving the following equation:

∂ t u -∇ ⊥ • ∇∆u = 0, in D, u = f, on ∂D, (1.1)
where ∇ ⊥ denotes the perpendicular gradient (-∂ y , ∂ x ) and f is the initial damaged image. This model was the subject of other improvements in 12 based on the Navier-Stokes equations. IOther authots 39,41,42 considered the Euler's elastica functional

Multiscale fourth-order model for image inpainting 3

to minimize the following energy:

∥u -f ∥ L 2 (Ω) + ∫ K (a + bκ 2 ) dH 1 (x), (1.2)
where H 1 denotes the Hausdorff measure, K is a closed regular subset of Ω and κ = ∇•(∇u/|∇u|) is the curvature of the level sets γ r := {x ∈ K | u(x) = r} and a, b being two positive constants. In this high-order variational model, the regularization term (Euler's elastica) combines the total variation, sensitive to the length of the isolines, and the square of the curvature, which favors curved lines rather than the straight ones. Minimizing (1.2) leads to a highly nonlinear PDE and therefore its numerical solution is usually a non trivial task. It was the subject of many numerical investigations 2,30,32,37 . In the same spirit, Esedoglu and Shen proposed in 32 the Mumford-Shah-Euler image inpainting model which is a high-order correction of the Mumford-Shah model, incorporating the Willmore energy (integral of the square of the curvature). They propose a numerical scheme based on the Γ-convergence approximation of Ambrosio and Tortorelli 3,4 . Another high-order approach was proposed in 16 based on the Cahn-Hilliard equation for binary images inpainting. This semi-linear fourth-order PDE originally introduced in material sciences by John W. Cahn and John E. Hilliard in 22 describes the evolution of an interface separating two stable states. A variant named T V -H -1 model was proposed in 21 as a generalization to the gray-scale images. In 28,29 , the authors studied the existence of global attractors to some generalized Cahn-Hilliard equations and the effect of the non-linear potential on the model.

The multiscale approach

It ressorts from this overview that most high-order approaches to the inpainting problem lead to the minimization of an energy of type (1.2). Such models are generally highly nonlinear and require some sophisticated discretizations to be solved numerically. In another hand, the easiest way to obtain a fourth-order PDE is to minimize ∫ D |∆u| 2 dx 18,23,43 which leads to solve the isotropic fourth-order (stationary) diffusion equation:

∆ 2 u = 0, in D, u = f + another boundary condition on ∂D. (1.3)
Of course, because of its strong smoothing effect, this model cannot capture some relevant features of the image like the corners and the edges which belong to the singular set of the image f . In this article, we consider the following equation:

       ∂ t u + ∆(∆ α u) + λ D (u -f ) = 0, in R + × Ω, u = ∆ α u = 0, on R + × ∂Ω, u(0, x) = f, in Ω, (1.4)
where ∆ α u = div(α(x)∇u) and λ D = λ 0 χ Ω\D for λ 0 ≫ 0 and χ Ω\D is the indicator function of the sub-domain Ω \ D. The values of the diffusion function α, which encodes different scales in the image, are dynamically and locally chosen in order to control the "amount" of smoothing of the operator. Note that the homogeneous boundary condition for u is not a restrictive condition and general boundary data can be handled by a lifting operator.

As the model remains linear with respect to u, it seems not well suited for capturing fine geometric structures of an image at a first glance. However, we prove that this approach allows us to restore such relevant features. The reason for this, is that the adaptive process which is the key of the method is in fact a nonlinear process, where the values of α are modified during the resolution. Together with a mesh adaption technique, we construct a nonlinear discrete approximation to make the method sensitive to the low-dimensional sets contained in the singular set of u. In addition, the process turns out to be Γ-convergent. Loosely speaking, the method allows us to adjust the model by modifying locally the diffusion coefficients instead of using a priori some nonlinear diffusion law or some nonlinear potential like in the Cahn-Hilliard based equations. This yields -and it is not a small advantage-a simple numerical method which consists of solving only linear problems on u and updating the diffusion parameters from the informations given by some error indicators. In addition, we emphasize that a strong point in this approach is that it conciliates a local character of the parameter selection which enhance the accuracy and the global character given by the PDEs framework which ensures the coherence of the reconstruction for large features. Thus, it constitutes, as showed by the numerical simulations a balance between the performances of the local methods, mainly based on "copy-paste" procedures and well suited for textured images, and the performances of global methods (typically PDEs).

Organization of the paper

The remainder of this article is organized as follows: In Section 2, we prove by standard variational techniques and a fixed point theorem the existence of H 1solutions for the image inpainting problem. We obtain our equation as the gradient flow of an energy based on an H -1 -fidelity term. In Section 3, we describe in details the adaptive strategy. In particular, we show that it is a two-step approach where, a mesh-adaptation based on the metric error indicator is used first to fit the geometry of the computed solution, and second a residual type error indicator is used to locally select the value of α. We perform the Γ-convergence analysis of this process and we show that the solution u generated by the adaptive strategy approximates a solution of a new model which combines the Mumford-Shah functional and the H -1 -fidelity term. In Section 4, we recall the Cahn-Hilliard equation for image inpainting that we will use for the comparison with our model. Finally, in Section 5 we implement our approach and treat several numerical examples to test its efficiency and robustness. and measurable. In fact, we recall the following regularity result for the operator -div(α∇(•)): 9 Proposition 2.1. There exists a constant c depending only on the geometry of Ω, such that a weak solution u of the associated Dirichlet problem belongs to H s+1 (Ω), for all real numbers s < s 0 , where s 0 is given by:

s 0 = min { 1 2 , c| log(1 - α m α M )| } .
Following the methodology used in 21 , we now prove the existence of the H 1 -weak solution of the stationary problem by means of Schauder's fixed-point theorem.

Proposition 2.2 (Schauder's fixed-point theorem 33 ). Let X be a Banach space and K ⊂ X be a compact and convex set of X. Assume that the operator T : K -→ K is continuous. Then T admits a fixed point.

For R > 0, to be chosen later, we define the convex set:

V = {u ∈ H 1 0 (Ω); ∥u 2 ∥ ≤ R},
and we consider the following minimization problem:

min u∈H 1 0 (Ω) J (u, v), (2.6)
where

J (u, v) =    ∫ Ω α(x) 2 |∇u| 2 dx + 1 2λ 0 ∥λ 0 u -λ D f -Λ D v∥ 2 -1 , if u ∈ V, +∞ , otherwise, (2.7)
with f and v ∈ L 2 (Ω) and where we set Λ D = λ 0λ D for brevity.

Proposition 2.3. Let v ∈ L 2 (Ω), the functional J (•, v) admits a unique minimizer u ∈ V .
Proof. The functional J (•, v) is strictly convex. In fact, let u 1 and u 2 be two functions in H 1 0 (Ω) such that u 1 ̸ = u 2 and t ∈]0, 1[, we have:

tJ (u 1 , v) + (1 -t)J (u 2 , v) -J (tu 1 + (1 -t)u 2 , v) = t(1 -t) 2 [∫ Ω α(x)|∇u 1 | 2 dx + ∫ Ω α(x)|∇u 2 | 2 dx -2 ∫ Ω α(x)∇u 1 • ∇u 2 dx ] + t(1 -t) 2λ 0 [ ∥λ 0 u 1 -λ D f -Λ D v∥ 2 -1 + ∥λ 0 u 2 -λ D f -Λ D v∥ 2 -1 ] - t(1 -t) λ 0 [∫ Ω ∇∆ -1 (λ 0 u 1 -λ D f -Λ D v) • ∇∆ -1 (λ 0 u 2 -λ D f -Λ D v)dx ] = t(1 -t) 2 ∫ Ω α(x)|∇(u 1 -u 2 )| 2 dx + t(1 -t) 2λ 0 ∥λ 0 (u 1 -u 2 )∥ 2 -1 > 0. Since L 2 (Ω) → H -1 (Ω), we get ∥v k ∥ -1 ≤ C∥v k ∥ L 2 (Ω) ≤ CR and also ∥λ D f ∥ -1 ≤
C ′ for some given constants C, C ′ > 0. Accordingly, we obtain the following estimate:

J (u k , v k ) ≤ C ′ + CR 2 ,
and then (u k ) k≥0 is uniformly bounded in H 1 (Ω). Thus, we can consider a convergent subsequence

u k j ⇀ j→∞ u ∈ H 1 (Ω) and u k j -→ j→∞ u in L 2 (Ω). Hence, the unique (weak) solution T (v k ) = u k of: { -∆ α u k -∆ -1 (λ D (f -u k ) + Λ D (v k -u k )) = 0, in Ω, u k = ∆ -1 (λ D (f -u k ) + Λ D (v k -u k )) = 0, on ∂Ω,
weakly converges to the unique weak solution u of (2.8). From the uniqueness of the solution in V , 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. In addition, this solution satisfies (2.3), or equivalently, is a stationary solution of the problem (1.4).

Remark 2.1. The solution u of equation (2.2) verifies the estimate:

α m ∥∆u∥ 2 + λ 0 2 ∥u∥ 2 ≤ λ 0 2 ∥f ∥ 2 .
Therefore, we may choose R = ∥f ∥ 2 .

The evolution equation

Now, observe that the previous existence proof of a stationary solution applies to show that the following unbounded operator (in L 2 ):

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

⟨A(u), u⟩ 2 ≥ I ∑ ℓ=1 α m ⟨∆u, ∆u⟩ L 2 (Ω ℓ ) + ⟨λ D u, u⟩ 2 ≥ 0,
which means that it is monotone. Thus it follows from the theory of maximal monotone operators 20 , that the evolution problem (1.4) admits a unique solution u ∈ L 2 (0, T ;

H 1 0 (Ω) ∩ H 1 (0, T ; H -1 (Ω)) ∩ C([0, T ] ; L 2 (Ω)).
As in the stationary case, problem (1.4) can be splitted into two second-order equations by introducing an auxiliary function w such that:

           ∂ t u -∆w + λ D (u -f ) = 0, in R + × Ω, -∆ α u = w, in R + × Ω, u = w = 0, on R + × ∂Ω, u(0, x) = f (x), in Ω.
(2.10)

Multiscale fourth-order model for image inpainting 9

The weak formulation reads then: Find a pair (u, w) ∈ L 2 (0, T ; H 1 0 (Ω)∩H 1 (0, T ; H -1 (Ω)))×L 2 (0, T ; H 1 0 (Ω)), u(0, x) = f (x), such that:

{ ⟨∂ t u, ϕ⟩ 2 + ⟨∇w, ∇ϕ⟩ 2 + ⟨λ D u, ϕ⟩ 2 = ⟨λ D f, ϕ⟩ 2 , ∀ϕ ∈ H 1 0 (Ω), ⟨α∇u, ∇ψ⟩ 2 -⟨w, ψ⟩ 2 = 0, ∀ψ ∈ H 1 0 (Ω).
(2.11)

Let u the solution of problem (1.4). Then, it is easy to verify that the pair (u, -∆ α u), is a weak solution of (2.11). We consider another solution (u 1 , w 1 ) ∈ H 1 0 (Ω)×H 1 0 (Ω) of the system (2.11), we then have:

{ ⟨∂ t (u -u 1 ), ϕ⟩ 2 + ⟨∇(w -w 1 ), ∇ϕ⟩ 2 + ⟨λ D (u -u 1 ), ϕ⟩ 2 = 0, ∀ϕ ∈ H 1 0 (Ω), ⟨α∇(u -u 1 ), ∇ψ⟩ 2 -⟨(w -w 1 ), ψ⟩ 2 = 0, ∀ψ ∈ H 1 0 (Ω).
Let (ζ ℓ ) I ℓ=1 be a partition of unity associated to the decomposition (Ω ℓ ) I 1 , and picking ψ = α -1 ℓ ζ ℓ (ww 1 ), in the second equation, we have the identity:

∫ Ω ℓ α -1 ℓ ζ ℓ (w -w 1 ) 2 dx = ∫ Ω ℓ α ℓ ∇(u -u 1 ) α -1 ℓ ∇ζ ℓ (w -w 1 ) dx, ∀l = 1, • • • , I.
(2.12) Summing up, we get (after integrating by parts twice the right-hand side):

I ∑ ℓ=1 α -1 ℓ ∫ Ω ℓ ζ ℓ (w -w 1 ) 2 dx = ∫ Ω ∇(u -u 1 ) ∇(w -w 1 ) dx ≥ 0.
(2.13) By choosing the test function ϕ = uu 1 in the first equation and using (2.13) and the positivity of λ D , we obtain:

⟨∂ t (u -u 1 ), u -u 1 ⟩ 2 = -⟨∇(u -u 1 ), ∇(w -w 1 )⟩ 2 -⟨λ D (u -u 1 ), (u -u 1 )⟩ 2 ≤ 0. 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 which implies that w = w 1 . Thus, problems (1.4) and (2.11) are equivalent.

Semi-discrete time problem

For the discretization of the time derivative ∂ t u, we use the forward Euler scheme.

Let

uu old ∆t be an approximation of ∂ t u, where ∆t is a time step, u old and u are the solutions at time t old and t = t old + ∆t, respectively. Therefore, time discretization together with splitting scheme (2.10) leads to the following time-stepping problem:

   ⟨ u -u old ∆t , ϕ⟩ 2 + ⟨∇w, ∇ϕ⟩ 2 + ⟨λ D u, ϕ⟩ 2 = ⟨λ D f, ϕ⟩ 2 , ∀ϕ ∈ H 1 0 (Ω), ⟨α∇u, ∇ψ⟩ 2 -⟨w, ψ⟩ 2 = 0, ∀ψ ∈ H 1 0 (Ω).
(2.14)

Multiscale fourth-order model for image inpainting 11

3.1. Γ-convergence analysis and the adaptive algorithm

In this section, we will prove that the approach considered in the article allows us to approximate, in the Γ-convergence sense, a new model that couples a Mumford-Shah functional with an H -1 -term, which we will call MS-H -1 .

Γ-convergence analysis A Γ-convergence study of this adaptive strategy was presented in 11 for optic flow estimation. The authors proved that this algorithm is equivalent to the adaptive one introduced by Chambolle-Dal Maso 25 and Chambolle-Bourdin 24 where a numerical discrete approximation of the Mumford-Shah energy was proposed. This approach, based on finite-element discretization and adaptive mesh strategy, is a good approximation, in the Γ-convergence sense 19 of the Mumford-Shah energy 40 (see 25 for more details). We briefly recall the results and the numerical approximation of this 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 25 introduced the following minimization problem:

G ϵ (u) = min T∈Tϵ(Ω) Gϵ (u, T), (3.2)
where

Gϵ (u, T) =    ∑ K∈T |K ∩ Ω| 1 h K g(h K |∇u| 2 ), u ∈ V ϵ (Ω), T ∈ T ϵ (Ω), +∞, 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 1 ).

Remark 3.1. If F : X -→ [-∞, +∞] is continuous and (G ϵ ) ϵ Γ-converges to G then (F + G ϵ ) ϵ Γ-converges to F + G.
From the result of the Γ-convergence of G ϵ to G, see 25 , the continuity of the second term of the functional in L 2 (Ω) (which follows from the continuity of ∆ -1 , i.e., the stability in the elliptic problems) and Remark 3.1, we have: Proposition 3.1. Let f and v in L 2 (Ω) be two given functions and ϵ > 0 be a positive parameter. Therefore, the sequence of functionals

G ϵ (u ϵ ) + 1 2λ 0 ∥λ 0 u ϵ -λ D f -(λ 0 -λ D )v∥ 2 -1 , Γ-converges for ϵ -→ 0 in the topology of L 2 (Ω) to G(u) + 1 2λ0 ∥λ 0 u -λ D f -(λ 0 -λ D )v∥ 2 -1 .
In the theorem we note that the variable α do 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 25 that the minimization of G ϵ is equivalent to the minimization of the following functional:

G ′ ϵ (u, v, T ϵ ) = ∑ K∈T ϵ |K ∩ Ω| 1 h 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 ϵ . For a fixed u, it follows from standard convex analysis that the minimizer over each v is explicitly given by:

v K = g ′ (h K |∇u| 2 ). (3.3)
Choosing α = v leads to the adaptive strategy that we describe now.. Remark 3.2. The analysis presented here is carried out with the Neumann boundary conditions on u, which is the framework used in 11,24,25 for a denoising problem or optic flow estimation. The application to the Dirichlet case requires some (tedious but non-essential) modifications and the result still holds.

Adaptive procedure

For each element K ∈ T h , the following local discrete energy:

η K = α 1 2 K h 1 2 K ∥∇u h ∥ L 2 (K) , (3.4)
contains some information on the error distribution of the computed solution u h . In fact, the discontinuities (edges) are contained in regions where the brightness changes sharply and consequently where this error indicator is large. Moreover, it may be proven that the gradient of u h captures this change in brightness and its magnitude provides an information about the "strength" of the edges (see 10 ). Thus, the quantity (3.4) acts as an edge detector and locates such regions. Furthermore, this local error indicator is in some sense equivalent to the energy norm mostly used as an edge detection in the topological gradient based-methods 6,7,8 . This particularity makes it well suited to control and locally select the diffusion coefficient α using the following algorithm:

Algorithm

(1) Start with the initial grid T 0 h corresponding to the image. (2) Adaptive steps:

• Compute u 0 h on T 0 h with a large constant α = α 0 by solving (3.1). • Build an adapted mesh T 1 h (in the sense of the finite element method, i.e., with respect to the parameter h) with a metric error indicator (to solve (3.2)).

• In the triangles where η K is large (with respect to its mean value), we perform a local choice of α(x) on T 1 h to obtain a new function α 1 (x) following (3.3).

(3) Go to steps 1. and 2. and compute u 1 h on T 1 h .

We emphasize that we have to solve only linear problems on u and then to update α. Actually, in place of introducing a function g like in the theorem, we use the following formula for the update of α in each triangle K:

α k+1 K = max      α k K 1 + τ * (( η K ∥η∥ ∞ ) -0.1 ) + , α trh      , (3.5)
where α trh is a threshold parameter and τ is a coefficient which controls the rate of decrease in α, u + = max(u, 0). Here η is the piecewise-constant function such that η| K = η k , ∀K ∈ T 1 h .

Remark 3.3.

There are other several possibles choices of functions g to compute α.

The formula (3.5) corresponds to min(t, g ∞ ) and is the one we implemented for the numerical computations here. Loosely speaking, it may be understood as follows: in regions of high gradients, it decreases the values of α when the error indicator deviates more then 10% from its mean value. α decreases nearly as a geometric sequence with the iteration number, until a given threshold is attained.

Let us give more details on the implementation of this algorithm. First, we build an adapted mesh T 1 h as follows: close to the jump sets of u h , the error is large, we then cut the element K into a finite number of smaller elements to decrease such an error and to fit the edges, while, far from these jump sets, there is no restriction on how to choose the triangles and the initial grid is coarsened. The adopted meshes have small number of degrees of freedom in the homogeneous regions which makes the method considerably fast. Second, we perform an "optimal" choice of the function α following the map furnished by the error indicator (η K ) K∈T h , on each element K and in accordance with (3.5) in order to correctly approximate edges.

To conclude this section let us try to interpret in term of a PDE the Γ-limit G. We may say that the sequence (u ϵ , α ϵ ) constructed by the adaptive algorithm is Γ-convergent to a Mumford-Shah-H -1 functional (MS-H -1 ) and that it solves

∆(∆ α ϵ u ϵ ) + λ D (u ϵ -f ) = 0, in Ω,
after eliminating w ϵ from the mixed formulation (3.1). Of course there are others Γ-convergent approximation functionals to the MS-H -1 , e.g. The Cahn-Hilliard model.

Cahn-Hilliard equation

For the sake of completeness, we will make a comparison with the Cahn-Hilliard model 15,16,17 that we now recall. The Cahn-Hilliard equation 22 originally refers to John W. Cahn and John E. Hilliard 22 and was introduced to phenomenologically describe phase separation, i.e., the change from one state to another. It is a fourthorder semi-linear PDE and is obtained as the H -1 -gradient flow of the following Ginzburg-Landau energy:

1 2 ∫ Ω |∇u| 2 dx + 1 ϵ 2 W (u), (4.1)
where phase separation is modeled by the smooth free energy, e.g., double-well potentials W such that:

W (u) = (1 -u 2 )u 2 , (4.2)
or a non-smooth free energy as follows:

W (u) = ψ(u) + I [0,1] (u), (4.3)
where ψ(u) = 1 2 u(1u), and

I [0,1] (u) := { 1 2 u(1 -u), 0 ≤ u ≤ 1, +∞, otherwise.
The analogy between the Cahn-Hilliard model in material sciences and Cahn-Hilliard model in image inpainting is that the two state phases are considered as two homogeneous regions and the interface is considered as an edge. In image inpainting, it was exploited by Bertozzi et al. in 16 and afterwards in 15,17,21 by considering the following equation:

   ∂ t u + ∆(ϵ∆u -1 ϵ W ′ (u)) + λ D (u -f ) = 0, in R + × Ω, u = f and ϵ∆u -1 ϵ W ′ (u) = 0, on R + × ∂Ω, u(0, x) = f, in Ω, (4.4)
which was obtained by incorporating the data fidelity term λ D (fu). Note that the classical Chan Hilliard equation (4.4) is only appropriate for two-scale (binary)

Multiscale fourth-order model for image inpainting 15

images inpainting due to the double well potential W which vanishes on only the values 0 and 1.

Remark 4.1. For the time discretization of Cahn-Hilliard equation, we use a convexity splitting scheme (see 17,21 ). The idea is to divide the energy functional into two parts; a convex part plus a concave one. The convex part is then treated implicitly while the concave part is treated explicitly.

Numerical experiments

In this work, all the PDEs are solved with the open source software FreeFem++ 35 . In all examples, the damaged/missed regions are delimited by the red contour. We give some examples for the application of our proposed approach to image inpainting. 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 λ 0 = 10 7 , and ∆t = 0.1 and an initial guess u(0, x) = 0.5, ∀x ∈ D. Curvature inpainting. We present the numerical results for the adaptive inpainting approach in order to illustrate the effectiveness of the proposed algorithm for the inpainting of curvature. In Fig. 1, we present the reconstruction of quarter of a circle. We display the evolution of the restored image for iterations 1, 5 and 20. In the first iteration (α is constant), we solved a biharmonic equation which gives a curved, but a very smooth (blurred), edge in D. By iterations, we can see the efficiency of the adaptation process in the damaged region where the edge was inpainted sharply by simultaneously keeping its curvature.

Other examples are presented in Fig. 2 and Fig. 4. The curvature in Fig. 2 is well inpainted which proves that our approach, based on fourth-order linear diffusion model, allows us to obtain a very interesting result that one might expect by solving comlpex PDEs like the Euler's elastica 42 , which is highly nonlinear and numerically difficult to solve. We give in Fig. 5 a zoom caption in the damaged region 2 which proves that the missing part is well restored and is very close to the original one.

"Real world"-image inpainting. The experiments in Fig. 6 show the efficiency of the proposed method in a real image inpainting. The portions of unknown pixels are 45%, 55% and 75%, respectively. From these experiments, we can see that the proposed model can successfully recover the inpainting domain even when up to 75% of pixels are unknown.

Comparison with the Cahn-Hilliard model. In Fig. 7, we have chosen the same image presented by Bertozzi, Esedoglu and Gillette (see 16 ). We give the result obtained by solving the Cahn-Hiliard equation and the one obtained using our approach. We display the evolution of the restored image at iteration 1, 5 and 10 which show that the edges are progressively and sharply approximated and the four corners are very accurately matched. Notice that the image for Cahn-Hilliard equation is computed in a two-steps process. In the first step, the authors solved their equation with a large value of ϵ, e.g., ϵ = 0.1, until the numerical scheme is close to a steady state. In this step, the level lines are continued into the missing domain. In a second step, they used the previous result as an initial time condition u 0 for a smaller ϵ (e.g., ϵ = 0.01) in order to sharpen the contours. This is an adaptive choice for ϵ, however, subjected to a hand tuning and being uniform in the entire domain.

Broken three bars. In Fig. 8, we show the inpainting result for a synthetic image which represents three broken bars forming a T -junction. We give the restored images using total variation, harmonic and biharmonic models, respectively. As excepted, total variation inpainting model is unable to connect the edges between the three bars because it does not fulfill the connectivity principle. Harmonic and biharmonic inpainting models produce a smooth solution u in D, blurring the edges. However, the proposed approach is able to sharply connect the edges between three bars while enforcing the curvature.

Conclusion

In this article, we have investigated a multiscale approach to a diffusion operator of fourth-order for the image inpainting. The model is linear with respect to the variable u representing the image to be reconstructed and is depending on a spatially varying diffusion function. We introduced an adaptive approximation procedure based on some a posteriori error indicators, which allows us to select the diffusion coefficient using the available data and the informations on the computed solution. We analyzed this approach from the variational point of view and we established its connections with a Mumford-Shah-like energy, in the sense of the Γ-convergence. We have implemented the considered model to test the method, we have also made some comparisons with existing approaches to demonstrate its capabilities. We have underlined, in the presented tests, the good quality in the recovery of low-dimensional sets (edges, corners) and of the curvature (of edges) in the inpainted zones. We emphasize that in our approach the adaptive selection of the diffusion coefficients is: (i) Fully automatic, i.e., no external intervention on the algorithm is required. (ii) Objective and a posteriori, i.e., the scale of α (the rate of decreasing its values) is explicitly obtained from computable quantities (error indicators) which are sensitive to the singularities of the solution. Thus, no a priori guess is necessary (iii) Local, i.e., made at each location in the computation domain (finite-element cell, thus the pixel scale). (iv) Low-cost, in the sens of the number of degrees of freedom, thanks to the coarsening of the mesh in the homogeneous parts of the domain.

In addition, the method may be improved straightforwardly by considering α, the diffusion function, as a matrix which introduces some anisotropy in the models and the overall approach is easy to implement in the framework of variational methods of approximation.

Multiscale fourth-order model for image inpainting 17 Combined second-and fourth-order derivatives and β and we give an analysis of the limit behavior of the adaptive algorithm. The analysis in this case shows that the solution generated by the adaptive strategy applied on model (5.1) approximates, in the Γ-convergence sense, a solution of a new model which couples a Mumford-Shah functional, the L 2 -norm of the solution of u and H -1 fidelity part.

We present several numerical examples, in image inpanting and restoration problems. We give a qualitative comparison for different values of the weights a and b which assert the need of higher-order derivatives in order to inpaint large damaged domain. The numerical results demonstrate the capability of recovering the recovery of low dimensional sets (edges, corners) and curvature in the image.

1. Introduction. The image inpainting consists in recovering the lost or deteriorated pixels in a given digital image from the available part. It constitutes a central problem in the image processing field and many of its applications (e.g. movies restauration, peinture, . . . ). In the last decades, a wide range of mathematical and engineering approaches were developed to solve this problem. Many engineering methods are based on a "copy-past" protocols where local patches are filled using a "likelihood principle" [19,28]. Such methods may be very expensive as they are local ones, besides they perform very well particularly as far as textured images are considered. The Partial differential Equations (PDEs), used extensively in the last two decades, have proved to be accurate and with reasonable costs in this field, particularly for less textured images, i.e. images with blocky geometric structures (cartoon images), we refere without being exhaustive to [10,17,20,23,27]). Unlike the previous approaches, they are global methods which make them useful for giving a coherent vue to the entire image (e.g. connecting large broken lines, . . . ) but less accurate for small textures. The second-order PDEs (Laplacian based methods) were initially used in the first works [16,17,20], since they provide a complete mathematical framework which consist of a theoretical analysis of the operator used to fill the missed parts and a wide range of numerical tools providing stable and accurate algorithms to compute the solutions. However, mostly they fail to reproduce some features in the image (curvature, corners,...) and they usually disconnect edges over large distances (violating the connectivity principle). The main reason for this failure comes from the fact that such geometric features are poorly represented with such operators. The higher-order diffusion models, based on operators of fourth-order are more suited in this case, that is why, they were proposed to overcome the limitations of the secondorder PDEs. They have proved to perform better in the inpainting problems as well as in the image restoration (see [1,9,11,13,18,20,26,27]. In fact, they damp the high frequencies (noise) faster than second-order based diffusion models, fulfill the connectivity principle and preserve the curvatures.

In this article, we consider a fourth-order PDE model for the inpainting and the restoration. The model couples second-and fourth-order derivatives: let f denotes the initial image and u the one we want to reconstruct, we solve

(1.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 weighting constant parameters, ∆ β u = div(β(x)∇u). The function λ D is a Lagrange multiplier and it is chosen such that λ D = λ 0 >> 0 in the image restoration problem, while in the image inpainting problem, it is chosen such that λ D = λ 0 χ Ω\D where χ Ω\D is the characteristic function of the domain Ω \ D.

The system (1.1) presents an anology with the Euler's elastica model which reads (in the intrinsic geometry)

∫ a + bκ 2
where κ is the curvature. This functional received much attention in the image analysis community [1,20,24,25,27] since it provides a high-order equation for (the level lines of) the image u and allows to highlight their curvatures. the variational interpretation of this functional as well as the PDE system obtained from its minimization (Euler-Lagrange equation) is highly non-linear and its numerical resolution calls for sophisticated algorithms [5,20,21].

In the model (1.1), α and β are two spatially varying functions, which will be chosen carefully in order to preserve during the reconstruction process the relevant features of the image like the edges, corners, . . . . Thus, we show that despite the linearity of the model with respect to u, a judicious choices of those functions yields a multi-scale approach which performs like some widely used non-linear models. Actually, the values of α and β are obtained locally and in adaptive way from a criterium given by the residual error estimators and the a posteriori analysis. These error indicators provide the information on the magnitude of the gradients and high-order derivatives, which are large in the regions containing the singular set of u and for high curvatures, hence allows one to decrease the diffusion coefficients in the operators and to sharpen and highlight such zones. Alternating the resolution of a linear PDE (wrt u) and the selection of the parameters α and β we obtain a dynamical model adjusted to fit the geometry of the reconstructed solution u and taking into account, some relevant features. Such procedure of selection is completely local and a posteriori, in the sense that no guess on u is used but only its finite element approximation. Therefore, it conciliates the global performances of the PDE approach and the local approach of patches and constitutes a strong point in the approach. Moreover, the numerical solving of the PDE (linear wrt u) is very easy and efficient.

We state the model, analyze it in the framework of the Γ-convergence theory. We show the minimization of a sequence of discrete energies which Γ-converge to the Mumford-Shah functional results in the adaptive algorithm which is convergent. This allows us also to set on a solid mathematical foundation the adaptive approach for the control of the diffusion coefficients in the PDE's approach

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

2. Existence of stationary H 1 -solution. Let the operator ∆ -1 β be the inverse of the negative operator -∆ β with homogeneous boundary conditions, i.e., for g ∈

V 0 = { h ∈ L 2 (Ω); ∫ Ω h dx = 0 } , we set u = ∆ -1 β g is the unique solution in H1 (Ω) = V 0 ∩ H 1 (Ω) of the following Neumann problem: (2.1) { -∆ β u = g, in Ω, ∂u ∂n = 0, on ∂Ω.
We define the inner product:

⟨u, v⟩ -1,β = ⟨β 1 2 ∇∆ -1 β u, β 1 2 ∇∆ -1 β v⟩ 2 ,
with corresponding norm ∥u∥ -1,β = ∥β 1 2 ∇∆ -1 β u∥ 2 in the space H1 (Ω). Throughout the paper, we suppose that the domain Ω is partitioned into I disjoint sub-domains (Ω ℓ ) ℓ such that α is given by the piecewise constant scalar function:

α = α ℓ 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 β ℓ .
We have set β as diffusion coefficient in the second-order term in the system (1.1). This choice allows us to provide the existence of a H 1 -weak solution for (1.1) when applying the operator ∆ -1 β . For the sake of clarity in the rest of the paper, we omit the time dependence of functions and we use the notation u instead of u(t).

We shall prove the existence of a H 1 -weak solution for the following stationary equation:

(2.2)

{ a∆ β (∆ α u) -b∆ β u + λ D (u -f ) = 0, in Ω, ∂u ∂n = ∂∆ α u ∂n = 0, on ∂Ω.
We suppose here that f ∈ V = {H 1 (Ω) ; ∫ Ω\D u dx = 0} and we consider the space Ṽ = {v ∈ Ṽ ; ∥v∥ 2 ≤ R}. For v ∈ L 2 (Ω), we define the following energy:

(2.3) J (u, v) =    a ∫ Ω α(x) 2 |∇u| 2 dx + b 2 ∥u∥ 2 2 + 1 2λ 0 ∥F D (u, v)∥ 2 -1,β , if u ∈ Ṽ, +∞, otherwise,
where

F D (u, v) = λ 0 u -λ D f -(λ 0 -λ D )v.
Proposition 2.1. The functional J (•, v) admits a unique minimizer u ∈ Ṽ.

Proof. The proof can be handled by standard techniques of calculus of variations. In fact, the functional J (•, v) is strictly convex and weakly lower semi-continuous. We consider a minimizing sequence (u n ) n∈N of J (•, v), i.e.,

J (u n , v) -→ n→∞ inf u∈ Ṽ J (u, v) = L.
Then there is a constant M > 0 such that ∫ Ω α(x)|∇u n | 2 dx ≤ M for all n ≥ 0, and also we have ||u n || 2 ≤ R (otherwise J (u n , v) would not be finite). Then, using the boundedness of α(x), we get that the sequence (u n ) n∈N is uniformly bounded in H 1 (Ω). Therefore, there exists a subsequence, still denoted (u n ) n∈N , such that u n ⇀ n→∞ u weakly in H 1 (Ω) and u n -→ n→∞ u in L 2 (Ω), such that ||u|| 2 ≤ R. Thanks to the continuity of the operator ∆ -1 : H -1 (Ω) -→ L 2 (Ω), we get:

J (u, v) ≤ lim inf n→∞ J (u n , v).
The limit u is then a minimizer for J (•, v), and since Ṽ is weakly closed, we have u ∈ Ṽ.

The minimizer u of J (•, v) fulfills the following optimality conditions:

(2.4)

{ -a∆ α u + bu -∆ -1 β (F D (u, v)) = 0, in Ω, ∂u ∂n = ∂(∆ -1 β (F D (u,v))) ∂n = 0, on ∂Ω. Now, let T : L 2 (Ω) -→ L 2
(Ω) such that for a given v ∈ L 2 (Ω), we associate T (v) = u the unique solution of (2.4). Therefore, if exists a fixed point u = v of T , it will be a solution of the equation (2.2). In fact, it seems easy that a fixed point u = v of T verifies:

(2.5)

{ -∆ α u + bu -∆ -1 β (λ D (f -u)) = 0, in Ω, ∂u ∂n = ∂(∆ -1 β (λ D (f -u)) ∂n = 0, on ∂Ω,
or equivalently, the system (2.2).

The existence of a fixed-point u = v of the operator T can handled by means of Schauder's Fixed-point theorem and we have:

Proposition 2.2. The operator T admits a fixed point u ∈ Ṽ. Moreover, u is H 1 -weak solution of the stationary equation (2.2).

Proof. We prove here the existence of a fixed point of T . Let v ∈ B(0, R) (where B(0, R) denotes the ball in L 2 (Ω) with center 0 and radius R). From Proposition 2.1, the energy (2.3) admits a unique minimizer u = T (v) in the space H 1 (Ω) such that u ∈ B(0, R). Since the embedding H 1 (Ω) → L 2 (Ω) is compact, the operator T then maps L 2 (Ω) → K, where K is a compact subset of L 2 (Ω). Thus we have:

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

J (u k , v) ≤ J (0, v), i.e., J (u k , v) ≤ 1 2λ 0 ∥λ D f + (λ 0 -λ D )v k ∥ 2 -1 . Since L 2 (Ω) → H -1 (Ω), we get ∥v k ∥ -1 ≤ C∥v k ∥ L 2 (Ω) ≤ CR and also ∥λ D f ∥ -1 ≤ C ′
for some given constants C, C ′ > 0. Accordingly, we obtain the following estimate:

J (u , vk) ≤ C ′ + CR 2 , C > 0,
and then (u k ) k≥0 is uniformly bounded in H 1 (Ω). Thus, we can consider a convergent subsequence u

k j ⇀ j→∞ u ∈ H 1 (Ω) and u k j -→ j→∞ u in L 2 (Ω). Hence, the unique (weak) solution T (v k ) = u k of: (2.6) { -a∆ α u k + bu k -∆ -1 β (F D (u k , v k )) = 0, in Ω, ∂u k ∂n = ∂(∆ -1 β (F D (u k ,v k ))) ∂n = 0, on ∂Ω,
weakly converges to the unique weak solution u of (2.4). From the uniqueness of the minimizer of (2.3), 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.

2.1. The evolution equation. Now, observe that the previous existence of a stationary solution applies to show that the following unbounded operator (in L 2 ):

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

⟨A(u), u⟩ 2 ≥ I ∑ ℓ=1 β m α m ⟨∆u, ∆u⟩ L 2 (Ω ℓ ) + I ∑ ℓ=1 β m ⟨∇u, ∇u⟩ L 2 (Ω ℓ ) + ⟨λ D u, u⟩ L 2 (Ω) ,
which means that it is monotone. Thus it follows from the theory of maximal monotone operators [12], that the evolution problem (1.1) admits a unique solution

u ∈ L 2 (0, T ; Ṽ ∩ H 1 (0, T ; H -1 (Ω)) ∩ C([0, T ] ; L 2 (Ω)).
The problem (1.1) can be splitted into two second-order equations by introducing an auxiliary function w such that:

(2.7)

         ∂ 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) ∈ Ṽ × H1 (Ω), u(0, x) = f (x) such that:

(2.8) { ⟨∂ t u, ϕ⟩ + b⟨β∇u, ∇ϕ⟩ 2 + a⟨β∇w, ∇ϕ⟩ 2 + ⟨λ D u, ϕ⟩ 2 = ⟨λ D f, v⟩ 2 , ∀ϕ ∈ H 1 (Ω), ⟨α∇u, ∇ψ⟩ 2 -⟨w, ψ⟩ 2 = 0, ∀ψ ∈ H 1 (Ω).
It is easy to verify that the pair (u, -∆ α u) is a weak solution of (2.8) where u is the solution for the problem (1.1).

To prove uniqueness, let the pair (u 1 , w 1 ) ∈ Ṽ × H1 (Ω) be another solution of the system (2.8), we then have:

     ⟨∂ t (u -u 1 ), ϕ⟩ 2 + b⟨β∇(u -u 1 ), ∇ϕ⟩ 2 +⟨β∇(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 picking ψ = β ℓ α -1 ℓ ζ ℓ (ww 1 ), in the second equation, we have the identity: For each ℓ, (2.9)

∫ Ω ℓ α -1 ℓ β ℓ ζ ℓ (w -w 1 ) 2 dx = ∫ Ω ℓ β ℓ α ℓ ∇(u -u 1 ) α -1 ℓ ∇ζ ℓ (w -w 1 ) dx.
Multiplying by T and summing up, we get (2.10) a

I ∑ ℓ=1 α -1 ℓ ∫ Ω ℓ ζ ℓ (w -w 1 ) 2 dx = a ∫ Ω β∇(u -u 1 ) ∇(w -w 1 ) dx ≥ 0.
By choosing the test function ϕ = uu 1 in the first equation and using (2.10) 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. 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 . Then, both w 1 and w are solutions of the weak formulation: Find r ∈ H1 (Ω) such that:

⟨∂ t u, ϕ⟩ + b⟨β∇u, ∇ϕ⟩ 2 + a⟨β∇r, ∇ϕ⟩ 2 + ⟨λ D u, ϕ⟩ 2 = ⟨λ D f, v⟩ 2 , ∀ϕ ∈ H 1 (Ω),
which implies that w = w 1 . Thus, problems (1.1) and (2.7) are equivalent. Proof. For a given v ∈ L 2 (Ω), u old ∈ H 1 (Ω) and dt > 0, we define:

Λ D = λ D + 1 dt and f ′ = Λ -1 D (λ D f + u old dt ).
Therefore, the problem (2.11) can be equally rewritten as follows:

(2.12)

{ b⟨β∇u, ∇ϕ⟩ 2 + a⟨β∇w, ∇ϕ⟩ 2 + ⟨Λ D u, ϕ⟩ 2 = ⟨Λ D f ′ , ϕ⟩ 2 , ∀ϕ ∈ H 1 (Ω), ⟨α∇u, ∇ψ⟩ 2 -⟨w, ψ⟩ 2 = 0, ∀ψ ∈ H 1 (Ω).
Using the same techniques as for the equation (2.2), we can prove the existence of H 1 -solution of: and such that the pair (u, w = -∆ α u) satisfies the system (2.12). To prove uniqueness, let (u 1 , w 1 ) ∈ Ṽ × H1 (Ω) be another solution of the system (2.12), we then have: { b⟨β∇(uu 1 ), ∇ϕ⟩ 2 + a⟨β∇(ww 1 ), ∇ϕ⟩ 2 + ⟨Λ D (uu 1 ), ϕ⟩ 2 = 0, ∀ϕ ∈ H 1 (Ω), ⟨α∇(uu 1 ), ∇ψ⟩ 2 -⟨(ww 1 ), ψ⟩ 2 = 0, ∀ψ ∈ H 1 (Ω).

As it was down for the system (2.8) and using the inequality (2.10), we obtain: b⟨β∇(uu 1 ), ∇(uu 1 )⟩ 2 + ⟨Λ D (uu 1 ), (uu 1 )⟩ 2 = -⟨∇(uu 1 ), ∇(ww 1 )⟩ 2 ≤ 0.

From the positivity of Λ D , we get: 3. Discrete problem 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) } ∩ H 1 (Ω).
The discretized version of the splitted problem (2.8) leads to finding a pair (u h , w h ) ∈ X h × X h solution of (3. where f h is a finite element approximation of f associated with T h . Since X h ⊂ H 1 , the well-posedness of the problem (3.1) follows from (2.11).

3.1. Adaptive procedure. The spirit of proposed adaptive strategy is mesh adaptation techniques and local choice of diffusion coefficients. In fact, mesh adaptation in finite elements method has become an efficient tool for the computation of a numerical solution of partial differential equations and variational problems. It allows more accuracy on the approximation of the solution by placing more mesh elements in these regions and less elements in regions of small error. In addition, the diffusion functions encode generally different scales in the image and, a bad chosen parameters lead to an over-smoothed and blurred edges. For that reason, we will make a dynamical and local choice of theses parameters in order to control the amount of smoothing of the operator in order to enhance the quality of the reconstructed edges.

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 discrete energies:

(3.2) η K = α 1 2 K h 1 2 K ∥∇u h ∥ L 2 (K) and η ′ K = β 1 2 K h 1 2
K ∥∇w h ∥ L 2 (K) , as error indicators. In fact, η K captures the change in brightness and, because the magnitude of the gradient, it then provides an information about the "strength" of the edges (see [6]). The error indicator η ′ K plays the same role as α in capturing the change in the Laplacin ∆u which gives a supplemtary informations about discontinuities of second kind. Thus, the two quantities in (3.2) act as edge detectors (see also [2,3,4]) which makes it well suited to control and locally select the diffusion coefficient α using the following algorithm: Algorithm 1 1. Start with the initial grid T 0 h corresponding to the image. 2. Adaptive steps:

• Compute u α 0 ,h on T 0 h with a large constant α = α 0 . • Build an adapted mesh T 1 h (in the sense of the finite element method, i.e., with respect to the parameter h) with a metric error indicator. • In the locations where η K is large (with respect to its mean value), we perform a local choice of α(x) and β(x) on T 1 h to obtain new functions α 1 (x) and β 1 (x). 3. Go to steps 1. and 2. and compute u α 1 ,h on T 1 h . During the adaptation, we use the following formula for each triangle K;

α k+1 K = max      α k K 1 + κ * (( η K ∥η∥ ∞ ) -0.1 ) + , α trh      ,
where α trh 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 η ′ .

3.2. Comments on the algorithm and its limit behavior. The adopted mesh is generated as follows: we cut the element K, close the edges of the image (i. e., where the error (η K ) K∈T h is large) into a finite number of smaller elements to decrease the error and to fit the edges, while, faraway, the grid is coarsened. Then, we make an "optimal" choice of the functions α and β, following the map furnished by the error indicator (η K ) K∈T h , on each element K. We decrease the values of the diffusion parameters in order to approximate correctly the edges.

An analysis of the limit behavior of this adaptive algorithm was presented in [7] for optic flow estimation and in [8] for image inpainting and restoration problems. The algorithm is proven to be equivalent to the adaptive one introduced by Chambolle-Dal Maso [15] and Chambolle-Bourdin [14] where a straightforward approximation of Mumford-Shah functional was proposed based on mesh adaptation techniques and adaptive local choice of the diffusion functions.

A similar analysis in our case shows that the solution generated by the adaptive algorithm allows to approximate in the Γ-convergence sense a model which couples a Mumford-Shah functional, the L 2 -norm of u and H -1 fidelity part.

Numerical examples.

In this work, all the PDEs are solved with the open source software FreeFem++ [22]. In all examples, the damaged/missed regions are delimited by the red contour.

4.1. Image restoration. We begin by testing our approach for an image denoising problem. In the left-hand plot of Fig. 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 (1.1).

We initialized the algorithm with a large value of α = β = 50 and we performed 20 iterations of the adaptive algorithm. We plot in Fig. 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. 2, which emphasizes the low cost of the method.

We tested our model on gray-scale image (200 × 200 pixels) in Fig. 3 where we displayed the original (noisy) image and the restored one. In Fig. 4, we vary the level of noise in the input images. In the top row, we superpose the three noisy images corrupted with Gaussian noise having zero mean and standard deviation 25 In the bottom row, we show the restored images. Curvature inpainting. In Fig. 5, we tested the adaptive algorithm for different values of the ratio b a in equation (1.1), in order to show the effects of each term (the fourth-and second-order one). We plot in Fig. 5 the restored images at the end of algorithm. It appears that if more weight is set on the second-order derivatives in (1.1), then the inpainted edge tends to be a straight line as expected (the length term is enforced). In fact, if we consider equation (1.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 [7,15]. 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 [17]. In Fig. 5, 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 [27], which is highly nonlinear and numerically difficult to solve.

Blood vessel inpainting. The experiment in Fig. 6 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. 7, we display the mesh at different adaptation iterations. The mesh is coarsed in the homogeneous regions which makes the overall numerical procedure fast.

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 (1.1). More precisely, we want to see if our approach fulfills the connectivity principle. In Fig. 11, we give the restored images using Total variation, harmonic and biharmonic models, 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 stripe. We display the restored images using the model (1.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.

Noisy images inpainting. We illustrate in Fig. 11 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 perform the restoration and inpainting actions in the same algorithm and we can see that edges are connected 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 correpondent 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 produces good solutions that one might obtain by using more complex higher order PDEs. We perform numerical experiments by denoising images that have been corrupted by Gaussian noise, , An adaptive approach for segmentation and TV edge-enhancement in the optic flow estimation, Submitted, (2015).

[8] Z. Belhachmi, M. Kallel, M. Moakher, and A. Theljani, Weighted harmonic and complex ginzburg-landau equations for gray value image inpainting, submitted, (2015).

Conclusion and perspectives

In this thesis, we have considered the image inpainting problem. Our main contribution has three folds. In the first chapter, inpainting was modeled as a Cauchy problem when the damaged domain touches the border of the image. In this case, inpainting was achieved by an extrapolation of data, unlike the classical methods where the interpolation was usually considered. We rephrased the Cauchy problem as an optimization problem in the framework of game theory. Moreover, the existence and uniqueness of the solution of the direct problem as well as the convergence of the linearized problems have been proved. In addition, numerical results confirm and suggest the good performance of this approach as an inpainting method.

In addition, we have investigated an adaptive approach for image inpainting based on a local selection of the different parameters in the models, and on mesh adaptation techniques. We started with the formulation of a linear variational model, and detailed its numerical implementation. We gave an analysis of the proposed approach which shows that our model approximates, in the sense of the Γ-convergence, the Mumford-Shah functional. We also presented two nonlinear diffusion models which are well rated in image processing in order to make a comparison. Numerical experiments on different images were performed and showed the efficiency of the proposed method. We can also say that the linear continuous model presented in this work gives the same final result that one might expect by solving a nonlinear one. As expected, second-order PDE methods are unable to reconstruct and keep the curvature because of its highest order (third order). This shortcoming was corrected by the extension of the method to fourth-order linear PDEs in chapter 4. First, we have investigated a fourth-order PDE where Dirichlet boundary conditions are used for purely image inpainting, i.e., where the known parts of the image are not noisy. Second, we introduced a straightforward combination of second-and fourth-order derivatives in analogy with Euler-elastica model. This new model was dedicated to joint image inpainting and restoration and we have used Neumann boundary conditions which are usually considered for image restoration problems. We gave a reliable and efficient numerical resolution for the two proposed PDEs based on Mixed finite element method. Then, we presented an adaptive approach for regularization parameter choice and we analyzed their limits behavior using Γ-convergence tools. We recalled the Cahn-Hilliard inpainting model in order to make a comparison with our approach. In the following, we address some interesting improvements of the works presented in this thesis that one can consider in the future:

The use of the Nash game approach gives an encouragement results in the resolution of the nonlinear Cauchy problem arising from image inpainting. It might be very interesting to extend the method to other inverse problem and other operators, e.g, Cauchy-Stokes problem, the coupling between Darcy and heat equations, etc. For instance, one can consider

Conclusion

the Nash game theory for texture image inpainting and restoration. In fact, image is usually decomposed in two geometric and texture parts, which are not of the same natures. The approach introduced by game theory can be most natural and be an original treatment in this case.

The results presented in Chapter 3, Chapter 4 and Chapter 5 can be improved by the use of more theoretical error indicators and anisotropic regularization. Furthermore, we have investigated a H 1 -solution of the proposed fourth-order models in Chapter 4 and Chapter 5. It can be very interesting to look at these problems as fourth-order ones and to search a H 2 -solutions.

The obtained results using the adaptive algorithm suggested its performance for the nontexture images. However, the method is less relevant for texture images, because it is based on local Partial Differential equations which are suited only to reconstruct geometric parts of the image. We think that this shortcoming could be corrected by using nonlocal methods or by using decomposition techniques and applying the adaptive strategy only on the geometry parts while the texture ones are inpainted using well suited other methods.

In the adaptation algorithm, the varying diffusion parameter acts as the phase filed function and it allows to extract edges. This important feature of the method can be exploited in image segmentation problems.

Figure 1 . 1 :

 11 Figure 1.1: The continuous and the digital images.

Figure 1 . 2 :

 12 Figure 1.2: Superposition of three channels to obtain a color image.

Figure 1 . 4 :Figure 1 . 5 :

 1415 Figure 1.4: Minimal inpainting to a 15th century panel painting, http: // www. parmaconservation. com/ laboratory. html .In this thesis, we are mainly concerned with the mathematical formulations of inpainting problems. The term inpainting is derived from the ancient art of repairing images in museums. It refers to professional artists who restore "by hand" scratched images that have suffered from paint loss, weakened canvas, tears, water damage, fire damage, and insect damage. The latter necessitates the integration of any repairs using various and specific inpainting techniques in order to keep them looking good. The development of computers in the mid-nineties of the last century and their unavoidable and daily use encouraged scientists to imitate inpainting and perform it automatically, and hence Digital Image Inpainting emerged as a subfield of image restoration. It is a very interesting task and has covered amazing diversity of applications in various fields (medical imaging, augmented reality, seismic images . . . ). It differs from denoising problem by the size of the damaged region which is in Pixel-size in the latter, whereas in inpainting, we may have a region totally damaged. Loosely speaking, let Ω ⊂ R d (d = 2, 3) denotes the entire image domain, the basic idea is

Figure 1 . 6 :

 16 Figure 1.6: The map ϕ : D -→ Ω \ D.

Figure 1 . 7 :

 17 Figure 1.7: Illustration of Masnou and Morel's approach for two T -junctions.

1 2 ( 2 00 2 00 2 00 1 2

 22221 Γ c ) and φ ∈ (H 1 (Γ c )) are given functions. Here (H 1 (Γ c )) denotes the dual space of H 1 (Γ c ), which consists of functions in H (Γ c ), vanishing on Γ i . Neumann boundary conditions on Γ c can be computed from information available in Ω\D.

1 2 ( 2 00 2 00 2 00 1 2

 22221 Γ c ) and φ ∈ (H 1 (Γ c )) are given functions. Here (H 1 (Γ c )) denotes the dual space of H 1 (Γ c ), which consists of functions in H (Γ c ), vanishing on Γ i . We assume the following conditions on the function k(•):

Theorem 2 . 2 .

 22 Under the assumptions (a), (b), (c) and (d), we have: (i) The nonlinear problem (3) admits a unique solution u ∈ H 1 (D).

Remark 1 . 1 2

 11 Fixing the Neumann condition φ ∈ (H 1 2 00 (Γ c )) , we cannot use the same techniques to prove that M = {f ∈ H 1 2 (Γ c ); f and φ are compatible } is dense in H (Γ c ) as in the case of Dirichlet boundary condition. In that case, we reduced the Dirichlet condition u| Γ c = f to a homogeneous Dirichlet condition on Γ c by using the function u * . We cannot do the same thing for Neumann condition due to the nonlinear term k(|∇(u)| 2 ). Proposition 2. For compatible data f and φ, the Cauchy problem (1) has a unique solution in H 1 (D).

Figure 1 .

 1 Figure 1. Illustration of the two typical situations for the inpainting problem: (a) when D is totally inside Ω, (b) when ∂D ∩ ∂Ω = ∅.

Figure 2 .

 2 Figure 2. Left: L 2 -error of Nash strategies S (k) -S (k-1) as function of Nash iterations. Right: L 2 -error of missing Dirichlet data τ k -u| Γ i u| Γ i as function of Nash iterations.

Figure 3 .

 3 Figure 3. From left to right: Original, damaged and restored images using Nash game (M SE = 4.63 × 10 -4 , P SN R = 33, 34).

Figure 4 .

 4 Figure 4. Damaged and original images.

Figure 5 .

 5 Figure 5. From left to right: Reconstructed image for α = 1, 10 -2 and 10 -6 , respectively.

Figure 6 .

 6 Figure 6. The PSNR and MSE as functions of different values of α.

Figure 7 .

 7 Figure 7. From left to right: Original, damaged and reconstructed images (M SE = 1.12 × 10 -3 , P SN R = 29.47).

Figure 8 .

 8 Figure 8. In the top row and from left to right: Input noisy image for 1, 3 and 5% noise level, respectively. In the bottom row and from left to right: Reconstructed image for 1, 3 and 5% noise, respectively. The PSNR is 31.06, 30.2 and 27.16, respectively. The MSE is 7.9 × 10 -4 , 9.6 × 10 -4 , and 1.9 × 10 -3 , respectively.

and Fig. 12 . 6 .Figure 9 .

 1269 Figure 9. Nash game method: Reconstructed Dirichlet (τ N , left) and Neumann (η N , right) data over Γ i for α = 1, α = 10 -2 and 10 -8 , respectively.

Figure 10 .

 10 Figure 10. Kozlov's method: Reconstructed Dirichlet (τ N , left) and Neumann (η N , right) data over Γ i for α = 1, α = 10 -2 and 10 -8 , respectively.

1 2 00Figure 11 .

 1211 Figure 11. Nash game method: Sensitivity with respect to noise; reconstructed Dirichlet (τ N , left) and Neumann (η N , right) data over Γ i for 1%, 3% and 5% noise.

)Figure 12 .

 12 Figure 12. Kozlov's method: Sensitivity with respect to noise; reconstructed Dirichlet (τ N , left) and Neumann (η N , right) data over Γ i for 1%, 3% and 5% noise.

Volume 1 ,

 1 Number 1, Page 1 WEIGHTED HARMONIC AND COMPLEX GINZBURG-LANDAU EQUATIONS FOR GRAY VALUE IMAGE INPAINTING ZAKARIA BELHACHMI, MOEZ KALLEL, MAHER MOAKHER, AND ANIS THELJANI (Communicated by )

1 2Z

 1 . BELHACHMI, M. KALLEL, M. MOAKHER, AND A. THELJANI models are formulated as a constrained optimization problem: minimize (1) R(u) given u = f + n in Ω\D,

4 Z

 4 . BELHACHMI, M. KALLEL, M. MOAKHER, AND A. THELJANI

6 Z

 6 . BELHACHMI, M. KALLEL, M. MOAKHER, AND A. THELJANI (4) Go to step (1) and compute u α 1 ,h on T 1 h . During the adaptation, we use the following formulae: for each triangle K(12)

10 Z

 10 . BELHACHMI, M. KALLEL, M. MOAKHER, AND A. THELJANI

12 Z

 12 . BELHACHMI, M. KALLEL, M. MOAKHER, AND A. THELJANI

Figure 2 .

 2 Figure 2. The mesh at iterations 1, 10 and 20, respectively.

Figure 3 .

 3 Figure 3. From left to right: Damaged, Harmonic and Harmonic & adaptation images, respectively.

Figure 4 .

 4 Figure 4. Original and damaged images.

Figure 5 .

 5 Figure 5. From left to right: restored images using total variation, harmonic and harmonic & adaptation (indicator η K ).

Figure 10 .

 10 Figure 10. The mesh at iterations 1 and 2.

Figure 11 .

 11 Figure 11. The mesh at iterations 10 and 20.

Figure 12 .

 12 Figure 12. (Left:) The error log ER = log ||u k αu exact || 2 , and (Right:) log of the number of degrees of freedom, as functions of the iteration numbers.

Examples 4 .

 4 InFig 17 andFig 16, we compare different methods when the damaged region contains a corner. The anisotropic nonlinear model (see the left-hand plots ofFig 17) produces a well contrasted edges separating the three homogeneous areas of the inpainting domain. But, the corner it self is not well captured. The harmonic model with adaptation (see the middle plots ofFig 17), approximates better the solution near the corner (still not well captured) and the edges are improperly contrasted. In the left-hand plots ofFig 17, we display the solution of the complex Ginzburg-Landau equation with the adaptation. First, the complexification allows us to diffuse more than two colors (0 and 1). Second, this model with

Figure 13 .

 13 Figure 13. Top row: Original and damaged images. Middle row: Restored image (Harmonic & adaptation) at iterations 1 and 10, respectively. The difference between the original image and the restored one at iterations 1 and 10, respectively.

Figure 14 .

 14 Figure 14. Zoom on a damaged region (40 × 50 Pixels): Original and restored images, respectively.

Figure 15 .

 15 Figure 15. From left to right: Damaged and noisy -restored at iteration 5 -restored at iteration 20.

Figure 16 .

 16 Figure 16. The original and damaged images.

Figure 17 .

 17 Figure 17. From left to right: The restored images using anisotropic nonlinear, Harmonic & adaptation and complex Ginzburg-Landau & adaptation models, respectively.

Figure 18 .

 18 Figure 18. The mesh at iterations 0, 7 and 30, respectively.

Figure 19 .

 19 Figure 19. From left to right: Original, damaged and restored images (Harmonic & adaptation), respectively.

Fig. 1 .

 1 Fig. 1. Top row: Original and damaged images. Middle row: Restored image using the model (1.4) and adaptation at iterations 1, 5 and 20, respectively. Bottom row: Mesh evolution at iterations 1, 5 and 20, respectively.

Fig. 2 .

 2 Fig. 2. Form left to right: Original, damaged and restored images using model (1.4) and adaptation.

Fig. 3 .

 3 Fig. 3. Mesh evolution at iterations 2, 5 and 10, respectively.

Fig. 4 .

 4 Fig. 4. Top row: Original, damaged and restored images using model (1.4) and adaptation., respectively. Bottom row: The evolution of the mesh at iterations 1, 5 and 20, respectively.

Fig. 5 .

 5 Fig. 5. Zoom on the damaged region 2: Original and restored.

Fig. 7 .

 7 Fig. 7. Top row: Damaged and restored images using model (1.4) and adaption (iterations 1 and 5). Bottom row: Restored images using model (1.4) and adaption (iterations 10 and 20) and Cahn-Hilliard model.

Fig. 8 .

 8 Fig. 8. Top row: Damaged image and restored ones using total variation and harmonic models, respectively. Bottom row: Restored images using biharmonic model and our approach at iterations 5 and 20, respectively.

{

  a∆ β (∆ α u) -b∆ β u + Λ D (uf ′ ) = 0, in Ω, ∂u ∂n = ∂∆ α u ∂n = 0, on ∂Ω,

b⟨β∇(u -u 1 )

 1 , ∇(uu 1 )⟩ 2 + ⟨Λ D (uu 1 ), (uu 1 )⟩ 2 = 0 It follows from the strict positivity of Λ D and β that u = u 1 and ∇u = ∇u 1 . Therefore, both w 1 and w are solutions of the weak formulation: Find r ∈ H1 (Ω) such that:b⟨β∇u, ∇ϕ⟩ 2 + a⟨β∇r, ∇ϕ⟩ 2 + ⟨λ D u, ϕ⟩ 2 = ⟨Λ D f ′ , v⟩ 2 , ∀ϕ ∈ H 1 (Ω),which clearly admits a unique solution and consequently w = w 1 .

Fig. 1 .

 1 Fig. 1. Original, noisy and restored images, respectively.

Fig. 2 .

 2 Fig. 2. The diffusion function α, the mesh and the error indicator η at convergence.

Fig. 3 .

 3 Fig. 3. The noisy and the restored images.

Fig. 4 .

 4 Fig. 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 image.

Fig. 5 .

 5 Fig. 5. Model (1.1) and adaptation: the ratio b a = 5, 1, 0.2 and 0.1, respectively.

Fig. 6 .

 6 Fig. 6. Blood vessel inpainting: the damaged image and restored images at iteration 5 and 20, respectively.

Fig. 7 .

 7 Fig. 7. Mesh at iteration 0, 5 and 20, respectively.

Fig. 8 .

 8 Fig. 8. Connectivity across large gaps: Total variation, harmonic and biharmonic models, respectively.

Fig. 9 .

 9 Fig. 9. Our approach: Connectivity across large gaps for the ratio b a is: 10, 5 and 1, respectively.

Fig. 10 .

 10 Fig. 10. Our approach: Connectivity across large gaps for the ratio b a = 0.2 and 0.1, respectively.

Fig. 11 .

 11 Fig. 11. Connectivity across large gaps for noisy image.

[ 6 ]

 6 Z. Belhachmi and F. Hecht, Control of the effects of regularization on variational optic flow computations, Journal of Mathematical Imaging and Vision, 40 (2011), pp. 1-19. [7]

  -Shah. Ensuite, on étend cette technique pour l'équation de Ginzburg-Landau et on présente des résultats numériques et des comparaisons avec un autre modèle non linéaire et anisotrope. Les EDP considérées dans le deuxième chapitre sont d'ordre deux, et il est connu que de tels modèles sont incapables de reconstruire des courbures et de rejoindre les contours de part et d'autre d'une partie endommagée ayant une taille importante. Pour remédier à ces inconvénients, on étend dans le chapitre 4 la technique du choix du paramètre de régularisation pour des EDP d'ordre quatre en désocclusion et restauration d'images. L'avantage des EDP d'ordre-quatre par rapport a celles d'ordre-deux est la connexion des contours de l'image de part et d'autre des parties manquantes, tout en préservant les courbures car elles fournissent plus d'informations sur leur géométrie. Contrairement au modèles classiques d'ordre supérieur, provenants des énergies de type Willmore, qui sont fortement non linéaires et difficiles à résoudre numériquement, on propose un modèle simple, et donc facile à résoudre numériquement. Le modèle est définie par le système d'équations suivant :

  Abstract. Image inpainting or disocclusion, which refers to the process of restoring a damaged image with missing information, has many applications in different fields. Different techniques can be applied to solve this problem.

	Inverse Problems and Imaging	doi:10.3934/ipi.2015.9.853
	Volume 9, No. 3, 2015, 853-874
	THE CAUCHY PROBLEM FOR A NONLINEAR ELLIPTIC
	EQUATION: NASH-GAME APPROACH AND APPLICATION TO
	IMAGE INPAINTING
	Moez Kallel, Maher Moakher and Anis Theljani
	Laboratoire LAMSIN, ENIT, University of Tunis El Manar
	B.P. 37, 1002 Tunis-Belvédère, Tunisia
	(Communicated by William Rundell)

  is a linear elliptic problem and has a unique solution which is the trivial solution. Hence, u 1 ≡ u 2 .

	858	Moez Kallel, Maher Moakher and Anis Theljani
	Inverse Problems and Imaging	Volume 9, No. 3 (2015), 853-874

  ⟩ 2 + b⟨β∇u h , ∇ϕ⟩ 2 + a⟨β∇w h , ∇ϕ⟩ 2 + ⟨λ D u h , ϕ h ⟩ 2 = ⟨Λ D f h , ϕ h ⟩ 2 , ∀ϕ h ∈ X h , ⟨α∇u h , ∇ψ h ⟩ 2 -⟨w, ψ h ⟩ 2 = 0, ∀ψ h ∈ X h ,

	1)    ⟨ u h -u old h ∆t	, ϕ h

Mathematics Subject Classification. Primary: 65M32, 68U10; Secondary: 35G25, 91A10. Key words and phrases. Inverse problem, Cauchy problem, game theory, Nash equilibrium, image inpainting.

c 2015 American Institute of Mathematical Sciences

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

BELHACHMI et al.

Acknowledgements

First, I would like to express my sincere gratitude to my supervisors Professor Zakaria BELHACHMI and Professor Maher MOAKHER and my co-supervisor Professor Moez KALLEL for the continuous support of my Ph.D study and related research, for their patience, motivation, and immense knowledge. They introduced me to interesting and Acknowledgements ii

Contents

Abstract i which produces a blurring of the edges at the first iterations. However, the simultaneous adaptive choice of α and λ allows us to recover the edges. We present the evolution of the restored image for iterations 5 and 20. In the 5th iteration (middle), the image is smoothed even in the regions where the data is available. λ is increased during the process following the formula (20) in order to fit the data-term.

"All of exact science is dominated by the idea of approximation."

Bertrand Russell (1872-1970).

Chapter 4

Multiscale fourth-order model for image inpainting and low-dimensional sets recovery

The chapter is devoted to the study of an adaptive choice of regularization parameters for fourth-order variational models in image inpainting and restoration problems, with emphasis on the recovery of low-order sets (edges, corners) and the curvature. It is based on the work presented in [12].

Summary

Second-order PDE-based methods are unable to connect edges over large distances (Connectivity Principle), and to reproduce features of higher-order (curvature, corners, etc) of the initial image due to the lack of information about the geometry of edges. These shortcomings gave rise to a new class of higher-order diffusion models which in general perform better, and they are effective and produce images with a high visual quality, specially nonlinear ones. However, most higher-order models in image restoration and inpainting make generally use of complex systems of nonlinear PDEs, e.g, Euler-elastica, Mumford-Shah-Euler, Blake and Zesserman, etc. In this chapter, we present a simple and effective linear fourth-order inpainting model given by the following equation:

where ∆ α u = ∇(α(x)∇u). The previous equation is completed with Dirichlet boundary conditions for the case of free-noise image inpainting. We give a theoretical analysis of the evolution equation and its correspondent stationary one. We prove the existence of H 1solutions and we give a simple and reliable discretization based on mixed finite elements. Then, we apply the adaptive approach detailed in Chapter 3 and we prove that, even the proposed PDE belongs to linear diffusion models which are not well-suited for capturing sharp edges, this approach allows us to recover fine features and to preserve curved edges in the image.

Multiscale fourth-order model for image inpainting 5

Image inpainting problem

We study in this section the image inpainting problem by considering the system (1.4). For the sake of clarity, in the rest of the paper we omit the time dependence of functions and we use the notation u instead of u(t).

H 1 -weak solution of the stationary equation

Assume that Ω is a bounded open set of R 2 with piecewise smooth boundary ∂Ω. We recall that the operator ∆ -1 is the inverse of the negative Laplacian with homogeneous Dirichlet boundary conditions, i.e., u = ∆ -1 g is the unique solution of:

Let H -1 (Ω) be the dual space of H 1 0 (Ω) with corresponding norm ∥u∥ -1 = ∥∇∆ -1 u∥ 2 and inner product ⟨u, v⟩ -1 = ⟨∇∆ -1 u, ∇∆ -1 v⟩ 2 , where ⟨•, •⟩ 2 and ∥ • ∥ 2 are the classical inner product and its corresponding norm in the space L 2 (Ω).

We consider the solution of the following stationary problem:

which can be seen as stationary solution of (1.4). Problem (2.2) can be rewritten as follows:

or equivalently as a coupled elliptic system:

(2.4)

The weak formulation of problem (2.3) is:

Throughout the paper, we suppose that he domain Ω is partitioned into I disjoint sub-domains (Ω ℓ ) ℓ such that α is given by the piecewise constant scalar function:

We denote α m = min Note that even if f ∈ L 2 (Ω), from the elliptic regularity results, one has w ∈ H 2 ∩ H 1 0 (Ω) however, the solution u might be only H 1 0 (Ω) since α is only bounded

Multiscale fourth-order model for image inpainting 7

Furthermore, J (•, v) is weakly lower semi-continuous in H 1 (Ω). We consider a minimizing sequence (u n ) n∈N of J (•, v), i.e.,

Then, there is a constant M > 0 such that ∫ Ω α(x)|∇u n | 2 dx ≤ M for all n ≥ 0, and also we have ∥u n ∥ 2 ≤ R (otherwise J (u n , v) would not be finite). Then, using the boundedness of α(x), we get that the sequence (u n ) n∈N is uniformly bounded in H 1 (Ω). Therefore, there exists a subsequence, still denoted (u n ) n∈N , such that u n ⇀ n→∞ u weakly in H 1 (Ω) and u n -→ n→∞ u in L 2 (Ω), with ∥u∥ 2 ≤ R. Thanks to the continuity of the operator ∆ -1 : H -1 (Ω) -→ L 2 (Ω), we get:

The limit u is then a minimizer for J (•, v). The uniqueness is guaranteed by the strict convexity of J (•, v).

The Euler-Lagrange equation corresponding to (2.6) reads:

Its weak formulation is:

(2.9)

Let T : L 2 (Ω) -→ L 2 (Ω) be the operator such that T (v) = u where u is the unique solution of (2.9). Therefore, if there exists a fixed point u = v of the operator, it will be a solution of problem (2.5).

Proposition 2.4. The operator T admits a fixed point u ∈ V . Moreover, u is H 1 -weak solution of the equation (2.2).

Proof. Let v ∈ B(0, R) (where B(0, R) denotes the ball in L 2 (Ω) with center 0 and radius R). From Proposition 2.3, the minimization problem (2.7) admits a unique minimizer u = T (v) in the space H 1 (Ω) such that u ∈ B(0, R). Since the embedding H 1 (Ω) → L 2 (Ω) is compact, the operator T then maps L 2 (Ω) → K, where K is a compact subset of L 2 (Ω). Thus we have:

where K is a compact and convex subset of L 2 (Ω). To apply Schauder's fixed-point theorem, it remains to prove that T is continuous in B(0, R). Let (v k ) k≥0 be a sequence which converges to v ∈ L 2 (Ω) 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:

Proposition 2.5. For a fixed u old ∈ H 1 0 (Ω), the problem (2.14) admits a solution

Proof. For a given v ∈ L 2 (Ω), u old ∈ H 1 0 (Ω) and ∆t > 0, we define:

) .

We then consider the following problem:

Similarly to Proposition 2.3, we can prove the existence of H 1 -solution of (2.15) such that the pair (u, -∆ α u) satisfies the system (2.14). To prove the uniqueness, let (u 1 , w 1 ) ∈ H 1 0 (Ω) × H 1 0 (Ω) be another solution of the system (2.14), we then have:

Following the same lines than for the system (2.11) and using the inequality (2.13), we obtain:

From the nonnegativity of λ ′ D , we get: ⟨λ ′ D (uu 1 ), (uu 1 )⟩ 2 = 0. Therefore u = u 1 and consequently w = w 1 .

Fully discrete problem and the adaptive strategy

We assume that the domain Ω is polygonal and we consider a regular family of triangulations T h made of element 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:

The discretized version of the splitted problem (2.11) leads to finding a pair

) where f h is a finite element approximation of f associated with T h . Since X h ⊂ H 1 0 (Ω), the well-posedness of the problem (3.1) follows from (2.14).

"Simplicity is the ultimate sophistication."

Leonardo da Vinci (1452-1519).

Chapter 5

Combined second-and fourth-order derivatives model for image inpainting and restoration

In this chapter, we propose a combined second-and fourth-order derivatives model for image inpainting and restoration problems. The material of this chapter is the subject of the work presented in [12], which is summarized bellow:

Summary

In this chapter, we propose a fourth-order PDE which couples second-and fourth-order derivatives in image inpainting and denoising. The model is given by the following equation:

where a > 0 and b ≥ 0 are two constants.

Model (5.1) might be considered as a simplified version of the Euler's elastica model, where the curvature and the length terms are replaced by the fourth-and second-order derivatives, respectively. The parameters a and b are used to control the trade off between the length and curvature in analogy with Euler's elastica model. Contrary to the choice of Dirichlet boundary conditions for the previous model detailed in Chapter 4, we used here Neumann boundary conditions which are most often used in the restoration of images. In fact, in some practical situations, the image to be inpainted suffers from noise in the intact regions which makes the PDE-based image inpainting algorithms very sensitive to noise. As the basic idea is to diffuse the known data accessible in neighborhood of the damaged regions, the noise present in the know portions will certainly affect the inpainted zones. Then, we must meet the twin objective of noise reduction and image completion by incorporating inpainting and denoising in parallel.

We give an theoretical analysis of the model and we present its discretization which is relatively simple contrary to the elastica model where the numerics are difficult and the algorithms are very slow (see [21,33]). Then, we make the adaptive choice of the functions α