Noémie Debroux

Carole Le Guyader
email: carole.le-guyader@insa-rouen.fr

Luminita A Vese
email: lvese@math.ucla.edu

Mila Searcher

Noémie Nikolova

Debroux

A Nonlocal Laplacian-Based Model for Bituminous Surfacing Crack Recovery and its MPI Implementation

Keywords: Γ -convergence, viscosity solutions, nonlocal second order operators, Message Passing Interface (MPI) parallelization Mathematics Subject Classification (2010) 35D40, 45E, 49J45, 65D18, 68W10, 68U10, 90C06

This paper is devoted to the challenging problem of fine structure detection with applications to bituminous surfacing crack recovery. Drogoul (2014) shows that such structures can be suitably modelled by a sequence of smooth functions whose Hessian matrices blow up in the perpendicular direction to the crack, while their gradient is null. This observation serves as the basis of the introduced model that also handles the natural dense and highly oscillatory texture exhibited by the images: we propose weighting

, u denoting the reconstructed image, by a variable that annihilates great expansion of this quantity, making then a connection with the elliptic approximation of the Blake-Zisserman functional. Extending then the ideas developed in the case of first order nonlocal regularization to higher order derivatives, we derive and analyze a nonlocal version of the model, and provide several theoretical results among which a Γ -convergence result as well as a detailed algorithmic approach and an MPI implementation based on a natural domain decomposition approach.

Introduction

Segmentation is a cornerstone step in image processing that attempts to reproduce the ability of human beings to track down significant patterns and automatically gather them into relevant and identified structures. More precisely, image segmentation consists in identifying meaningful constituents of a given image (e.g., homogeneous regions, edges, textures, etc.) for quantitative analysis or visualization purposes. As emphasized by Zhu et al. ([46]), this task is challenging and ill-posed since the definition of an 'object' encompasses various acceptations: it can be something material -a thing -or a periodic pattern, this heterogeneity entailing the design of suitable methodologies for each specific application. For a relevant overview of the existing literature on this topic, we refer the reader to [START_REF] Zhu | Beyond pixels: A comprehensive survey from bottom-up to semantic image segmentation and cosegmentation[END_REF] in which an exhaustive classification of segmentation methods into three main categories (fully supervised methods, semi/weakly supervised methods, unsupervised methods ) is provided, combined with a description / analysis of each methodology. In this paper, we focus on unsupervised methods in a variational formulation and phrased in the continuous domain. While a contour is classically viewed as the boundary of a non zero volume object and is thus defined as the set of intensity discontinuities with jump (detected using the spatial gradient information carried by the image), fine structures exhibit discontinuities without jump. To exemplify this latter statement, we can think of a curve γ ⊂ Ω, Ω being an open bounded subset of R 2 , and a related image I defined by I = 1 on γ and I = 0 on Ω\γ. The discontinuity related to the object of zero Lebesgue measure {I = 1} is a discontinuity without jump (see [START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF][START_REF] Drogoul | Méthode du gradient topologique pour la détection de contours et de structures fines en imagerie[END_REF] for further details). These thin/filament like structures cannot be captured efficiently with differential operators of order 1. In [START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF][START_REF] Drogoul | Méthode du gradient topologique pour la détection de contours et de structures fines en imagerie[END_REF], Drogoul provides a heuristic illustration of this observation as follows. Let f be the function defined by f (x) = 0 if x = 0 and f (x) = 1 if x = 0. This function can be approximated by the sequence of functions (f ε ) such

that f ε (x) = 0 if x ≥ ε 2 ε 3 |x| 3 -3 ε 2 |x| 2 + 1 if x < ε
. On the one hand, f ε (0) = 0, showing that the differential operator of order 1 does not recover the singularity at 0. On the other hand, f ε (0) = -6 ε 2 , meaning that f ε becomes singular at 0. This observation is then mathematically formalized in the two dimensional case through the following lemma taken from [START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF][START_REF] Drogoul | Méthode du gradient topologique pour la détection de contours et de structures fines en imagerie[END_REF] that we recall for the sake of completeness. Let D(R 2 ) be the space of C ∞ functions with compact support in R 2 and let D (R 2 ) be the space of distributions on R 2 .

Lemma 1 (Lemma 2.1 of [START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF]) Let ϕ : R 2 → R be a Lipschitz continuous function, and let (g h ) h>0 be a sequence of functions defined by g h (x) = 1 θ1(h) e and spec ∇ 2 g h (x) = -2 √ πh 3/2 , 0 , with spec(M ) denoting the eigenvalues of the matrix M . The associated eigenvectors on Γ are ∇ϕ(x), ∇ϕ(x) ⊥ , where ∇ϕ(x) = n is the unit normal of Γ .

Since filaments/cracks can be modelled by a Dirac distribution supported by a smooth curve Γ , it results from this lemma that any such structure can be approximated by a sequence of smooth functions whose Hessian matrices blow up in the perpendicular direction to Γ while their gradient are null. This observation serves as the core of the introduced material that falls within second order variational models, includes an additional variable encoding the singularities of the quan-

tity ∂ 2 u ∂x 2 1 2 + ∂ 2 u ∂x 2 2 2
(object we aim to recover) -the cracks being mainly oriented in the vertical or horizon-tal direction-, and that also handles the natural dense and highly oscillatory texture exhibited by the images.

More precisely, we propose weighting

∂ 2 u ∂x 2 1 2 + ∂ 2 u ∂x 2 2 2
, u denoting the reconstructed image function, by a variable v 2 with values in [0, 1] annihilating great expansion of this quantity so that the set of cracks is now modelled by {v 2 = 0}. This leads to a variant of the elliptic approximation of the Blake-Zisserman functional ( [START_REF] Blake | Visual Reconstruction[END_REF]) given by Ambrosio, Faina and March ( [START_REF] Ambrosio | Variational Approximation of a Second Order Free Discontinuity Problem in Computer Vision[END_REF]). Before introducing and analyzing in depth the proposed model, we review some prior works related to thin structure detection and enhancement. Recently, the question of differentiating objects within an image by contrast, size or structure -the size being the feature of interest in our case -has been addressed in [START_REF] Hait | Spectral Total-Variation Local Scale Signatures for Image Manipulation and Fusion[END_REF]. It relies on the spectral total-variation transform that yields a continuous, multi-scale, fully edge-preserving, local descriptor, referred to as spectral total variation local scale signatures. These signatures demonstrate nice properties such as sensitivity to size, local contrast and composition of structures, as well as invariance to rotation, translation, flip and linear illumination change. In [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], Frangi et al. tackle the issue of vessel enhancement with the eventual goal of vessel segmentation in a multiscale approach. It is based on the construction of a vessel likeliness measure that relies on the computation of the eigenvalues of the Hessian. In [START_REF] Aubert | Topological Gradient for a Fourth Order Operator Used in Image Analysis[END_REF][START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF], Aubert and Drogoul introduce a topological gradient-based method for the detection of fine structures in 2D, the major difference with our model lying in the introduction of an auxiliary variable encoding the crack-type singularities. Still in a variational setting, Bergounioux and Vicente ( [START_REF] Bergounioux | Parameter Selection in a Mumford-Shah Geometrical Model for the Detection of Thin Structures[END_REF]) propose a model derived from the Mumford-Shah functional and constrained by geometrical priors to perform the segmentation of tube-like structures with small diameter. A limitation of this method rests on its inability to handle junctions of tubes. Other variational approaches have been investigated, dedicated to specific applications, such as automatic road network extraction in [START_REF] Rochery | Higher Order Active Contours[END_REF] in which nonlocal regularizers enforcing straightness are applied, or detection/completion of thin filaments in noisy blurred 2D or 3D images in [START_REF] Baudour | Detection and Completion of Filaments: A Vector Field and PDE Approach[END_REF], relying on the Ginzburg-Landau energy. In [START_REF] Aubert | Detecting Codimension-Two Objects in an Image with Ginzburg-Landau Models[END_REF], Aubert et al. propose detecting image singularities of codimension greater or equal to 2, inspired again by Ginzburg-Landau models. Obviously, the range of the methods addressing this issue of thin structure detection does not limit to the variational framework: morphological approaches are developed in [START_REF] Zana | Segmentation of Vesell-Like Patterns Using Mathematical Morphology and Curvature Evaluation[END_REF] for automatic detection of vessel-like patterns, though sensitive to the noise type and timeconsuming. In [START_REF] Stoica | A Gibbs Point Process for Road Extraction from Remotely Sensed Images[END_REF], stochastic methods in which a thin network is simulated by a point process penalizing disconnected segments and favoring aligned pieces are elaborated. We emphasize that this paper is an extension of [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF] not reduced to a simple supplement, containing no overlap with [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF] and including significant theoretical and computational aspects. It completes [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF] since it includes substantial results, such as a Γ -convergence result (main result) relating the local basis functional with the derived nonlocal one, and an exhaustive and substantive PDE analysis, in the viscosity solution theory framework, of the resulting evolution equation satisfied by the variable Q introduced hereafter. It enriches [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF] with the derivation of an MPI parallelization of the designed algorithm inspired by the Nonlocal Means ( [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]), this algorithmic part being entirely new and numerically analyzed (in particular, the choice of the weights is discussed). The next section (Section 2) is dedicated to the mathematical material. The original local functional intertwining decomposition and fine structure detection is first introduced. Inspired then by the works of Bourgain et al. ( [START_REF] Bourgain | Another look at Sobolev spaces[END_REF]), Dávila ([17]) or Ponce ([35]) devoted to the design of nonlocal counterparts of Sobolev and BV semi-norms using radial mollifiers and integral operators, we extend their results to the case of second order operators (also inspired by the work of Lellmann et al. ( [START_REF] Lellman | Analysis and application of a non-local Hessian[END_REF])) and rephrase our original model in the obtained nonlocal setting. The main theoretical result -Γ -convergence result -connecting the nonlocal minimization problem to the local one is provided. Section 3 deals with the nonlocal algorithm that is devised, motivated by the Nonlocal Means approach ( [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]), to include additional information related to the image content via some weights. Optimality conditions are obtained, and a thorough analysis, in the viscosity solution framework, of the PDE resulting from the evolution equation in variable Q, including existence/uniqueness and regularity in space and time is carried out. Some limitations regarding the theoretical expression of the nonlocal component are then highlighted, which leads us to introduce a numerical compromise removing in particular the assumption of radiality of the kernel. The resulting algorithm is depicted and to improve the computation efficiency, an MPI parallelization of the code is proposed. This parallelization is driven by the natural geometry of the problem, making the partition of the image domain into subdomains supporting simultaneous local computations relevant. The paper concludes with Section 4 that intends to qualitatively assess the proposed algorithm on several experiments in comparison to prior works.

Mathematical modelling of the problem

Mathematical preliminaries

Let us first introduce some notations that will be useful throughout the paper. The image domain is denoted by Ω which is assumed to be a connected, bounded, open subset of R 2 of class C 1 . Let f : Ω → R be the two dimensional image. For theoretical purposes, we suppose that it is in L ∞ (Ω), which is not restrictive in practice. Let (e 1 , e 2 ) be the canonical basis of R 2 . We use dx (x = (x 1 , x 2 )) for integration with respect to the Lebesgue measure on R 2 and dt, ds, dh for various integrations with respect to the Lebesgue measure on R. The differentiation indices will be a pair α = (α 1 , α 2 ), where α i is the order of the partial derivative in the variable x i , and the total order of the derivative is denoted by |α| = α 1 + α 2 . We will use the shortened notation

D α u = ∂ |α| u ∂x1 α 1 ∂x2 α 2 .
Given an integer j ≥ 0, we define the family of spaces C j b (R

C j b (R 2 ) = {u ∈ C j (R 2 ) | ∀α ∈ N 2 , |α| ≤ j, ∃K α , D α u L ∞ (R 2 ) ≤ K α }.
For a positive real number 0 < λ < 1, the subspace

C j,λ b (R 2 ) consists of the functions in C j b (R 2 ) such that if |α| ≤ j, then ∃C α,λ , ∀x, y ∈ R 2 , |D α u(x) -D α u(y)| ≤ C α,λ |x -y| λ .
We also recall the definition of Sobolev spaces for k ∈ N, and 1 ≤ p ≤ ∞, W k,p (Ω) = {u ∈ L p (Ω) | D α u ∈ L p (Ω), ∀|α| ≤ k}, and the subspaces W k,p 0 (Ω) = {u ∈ W k,p (Ω) | γ 0 (u) = 0}, with γ 0 the trace operator. We then introduce the Hilbert space H(div) defined as 2 . Finally, we need the operator called normalized infinity Laplacian and defined as

H(div) = {σ ∈ (L 2 (Ω)) 2 | divσ ∈ L 2 (Ω)} endowed with the scalar product σ 1 , σ 2 H(div) = σ 1 , σ 2 (L 2 (Ω)) 2 + divσ 1 , divσ 2 L 2 (Ω) , ∀(σ 1 , σ 2 ) ∈ H(div)
∆ ∞ (u) = Du |Du| , D 2 u Du |Du|
. We now turn to the motivation and the depiction in depth of our model.

Original local basis model

As emphasized in the introduction, the bituminous surfacing images we consider in this work exhibit a natural dense and highly oscillatory texture (corresponding to white and black spots) that appears to be unnecessary and superfluous for the crack detection task. We thus want to discriminate and separate the texture component v that catches redundant and oscillatory patterns from the simplified version of the image denoted by u in which useless information has been removed, while keeping the crack. We therefore introduce a joint decomposition/segmentation model in which the thin-structure recognition task operates on the component u. Texture modelling consists in finding the best functional space to represent the oscillatory patterns and has been extensively studied, either on its own or as a close counterpart of image denoising (see [START_REF] Gao | Infimal convolution of oscillation total generalized variation for the recovery of images with structured texture[END_REF], [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF], [START_REF] Kim | Image recovery using functions of bounded variation and Sobolev spaces of negative differentiability[END_REF], [START_REF] Lieu | Image restoration and decomposition via bounded total variation and negative Hilbert-Sobolev spaces[END_REF], [START_REF] Meyer | Oscillating Patterns in Image Processing and Nonlinear Evolution Equations[END_REF], [START_REF] Ono | Cartoon-Texture Image Decomposition Using Blockwise Low-Rank Texture Characterization[END_REF], [START_REF] Osher | Image decomposition and restoration using total variation minimization and the H -1 norm[END_REF], [START_REF] Rudin | Nonlinear Total Variation Based Noise Removal Algorithms[END_REF], [START_REF] Schaeffer | A texture model based on a concentration of measure[END_REF], [START_REF] Schaeffer | A Low Patch-Rank Interpretation of Texture[END_REF] or [START_REF] Vese | Modeling Textures With Total Variation Minimization And Oscillating Patterns In Image Processing[END_REF]).

In this work, to model the texture, we propose to use the space G(R 2 ) introduced by Meyer ( [START_REF] Meyer | Oscillating Patterns in Image Processing and Nonlinear Evolution Equations[END_REF]), space of distributions v that can be written as v = divg where

g = (g 1 , g 2 ) ∈ (L ∞ (R 2
)) 2 and endowed with the norm defined by

v G(R 2 ) = inf g 2 1 + g 2 2 L ∞ (R 2 ) |v = divg .
This space is rather large to capture the texture exhibited by the bituminous surfacing images, while its norm is quite easily approximated by a numerical efficient expression. A further justification to use this functional space is its link with the notion of scale presented by Strong et al. in [START_REF] Strong | Scale Recognition, Regularization Parameter Selection, and Meyer's G Norm in Total Variation Regularization[END_REF]:

if v ∈ G(R 2 ), then v G = sup E⊂Ω | E v| P (E,Ω) and scale(x) = |E| |∂E| for x ∈ E ⊂ Ω, with P (E, Ω) denoting the perimeter of E in Ω.
Therefore, the stronger the penalization of v G is, the smaller the scale of the features kept in v is. Based on the following observation: by modelling the crack as a rectangle of k × 1 pixels with k 1, its scale behaves like 1/(2n) for an image of size n × n while small features of a pixel size have a scale behaving like 1/(4n), one can discriminate these small details from the crack by choosing adequately the penalization of the G norm. To approximate this norm, we rely on the work of Osher et al. [START_REF] Osher | Image decomposition and restoration using total variation minimization and the H -1 norm[END_REF] and introduce an auxiliary variable naturally stemming from Helmholtz-Hodge decomposition as follows: g = ∇Q + P, with P a divergence-free vector disregarded afterwards. An L 2 penalization term ensures the closeness between g and ∇Q and the minimization of the L ∞ norm is applied to ∇Q, yielding a problem related to the absolutely minimizing Lipschitz extension and to the infinity Laplacian. Equipped with this material, we propose combining the decomposition and Laplacian-based extraction tasks into a single variational framework as follows:

inf F ε (u, g, Q, v 1 , v 2 ) = f -u -div g 2 L 2 (Ω) + µ |∇Q| L ∞ (Ω) + γ 2 |g -∇Q| 2 L 2 (Ω) + ρ Ω (v 2 2 + κ ε ) |D (2,0) u| 2 + |D (0,2) u| 2 dx (1) + ξ ε Ω (v 2 1 + ζ ) |∇u| 2 dx + (α -β) G ε (v 1 ) + β G ε (v 2 ), with G ε (l) = Ω ε |∇l| 2 + (l-1) 2 4ε dx, κ ε , ξ ε , ζ suit-
able infinitesimals, and ε a given parameter. The first three terms cope with the decomposition of f into its texture component v such that v = divg, g being close to ∇Q, and a simplified version of it, u, including the crack, while the useless small oscillatory patterns have been removed. The remaining terms are devoted to the crack detection process on u. v 2 is close to 1 almost everywhere yielding a regularization of u, except where |D (2,0) u| 2 + |D (0,2) u| 2 is large, v 2 being close to 0 in this case. This term thus encodes the fine structures we aim to recover, that is to say the set of cracks is now modelled as {v 2 = 0}. Similarly, v 1 is close to 1 almost everywhere, yielding a regularization of u except in a neighborhood of the jump set of u where |∇u| 2 is large, v 1 being close to 0 in this case to preserve it. It encodes the discontinuity set of u, that is to say, the edges with jumps by {v 1 = 0}. The regularization on v 1 and v 2 ensures their smoothness and that they are close to 1 almost everywhere, while the terms combining v 2 2 and |D (2,0) u| 2 + |D (0,2) u| 2 (resp. v 2 1 and |∇u| 2 ) assure that they are close to 0 on thin structures (resp. on the jump sets). The fidelity term guarantees similarity between u and the observation f , while the component Q regulates the size of the features captured by divg. The minimisers are searched in the respective following functional spaces:

-g in H(div). -Q in Q ∈ W 1,∞ (Ω)| Ω Q dx = 0 -we are inter- ested in ∇Q rather than in Q itself, so this condition is not restrictive -. -v 1 in W 1,2 (Ω, [0, 1]) = {f ∈ L 2 (Ω, [0, 1]), | ∀α such that |α| ≤ 1, D α f ∈ L 2 (Ω)}. -u in W 1,2 0 (Ω) ∩ W 2,2
(Ω) enabling u to be extended by 0 outside the domain Ω without being too restrictive.

-v 2 in W 1,2 (Ω, [0, 1]).
The reasons for such requirements will be made clearer in the following. Nevertheless, these assumptions are reasonable and not restrictive if we assume for instance that the observed image f is with compact support. Existence of minimizers is guaranteed as stated in [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF]Theorem 11].

Towards a nonlocal related model

Inspired by the strength and robustness of nonlocal methods exemplified in many image processing tasks in [START_REF] Jung | Nonlocal Mumford-Shah regularizers for color image restoration[END_REF] (color image denoising, color image deblurring in the presence of Gaussian or impulse noise, color image inpainting, color image super-resolution, and color filter demosaicing), and also motivated by their mathematical aspects, we introduce a nonlocal version of our functional. Based on prior related works by Bourgain, Brezis and Mironescu [START_REF] Bourgain | Another look at Sobolev spaces[END_REF] introducing a new nonlocal characterization of the Sobolev spaces W 1,p , 1 < p < +∞, Aubert and Kornprobst [START_REF] Aubert | Can the nonlocal characterization of Sobolev spaces by Bourgain et al. be useful for solving variational problems?[END_REF] questioning the utility of this characterization to actually solve variational problems, Dàvila [START_REF] Dávila | On an open question about functions of bounded variation[END_REF] and Ponce [START_REF] Ponce | A new approach to Sobolev spaces and connections to Γ -convergence[END_REF] dedicated to the design of nonlocal counterparts of first order Sobolev and BV seminorms, we introduce a sequence of radial mollifiers {ρ n } n∈N satisfying:

∀n ∈ N, ∀x ∈ R, ρ n (x) = ρ n (|x|); ∀n ∈ N, ρ n ≥ 0; ∀n ∈ N, R ρ n (x) dx = 1; ∀δ > 0, lim n→+∞ +∞ δ
ρ n (r) dr = 0, and extend the previous results to second order operators. Our nonlocal extension is also deeply inspired by the formulation of a nonlocal Hessian and the novel characterization of higher Sobolev and BV-spaces introduced in [START_REF] Lellman | Analysis and application of a non-local Hessian[END_REF]. This new operator is derivative free and only requires the considered function u to be in an L p space as in our model. The main difference relies in our choice to treat independently the directional derivatives and to neglect the cross derivatives based on the empirical observation that the crack is most likely horizontal, making it not too restrictive to consider only the canonical directions. Besides, thanks to the theory of tempered distributions, one can show that if u,

∂ 2 u ∂x 2 1 , ∂ 2 u ∂x 2 2 ∈ L 2 (R 2 ), then u ∈ W 2,2 (R 2
). We thus introduce a sequence of functionals F ε,n depending on n and such that the component Ω (v 2 2 + κ ε ) |D (2,0) u| 2 + |D (0,2) u| 2 dx is approximated by an integral operator involving a differential quotient and the radial mollifiers. The resulting model inherits fine analytical properties, has the advantage of being numerically more tractable than [START_REF] Lellman | Analysis and application of a non-local Hessian[END_REF] and of being straightforwardly connected to our imaging problem. We would like to point out that the qualifying term "nonlocal" might sound inadequate in the sense that the "nonlocal counterpart" includes a parameter n, via the mollifier ρ n , that is destined to tend to infinity, concentrating the measure around the point of interest and removing in some way the nonlocal nature of the component. Nevertheless, it takes on mathematical interests since, to the best of our knowledge, this kind of approximation with an independent treatment of the directional derivatives for second order operators has not been investigated yet except in [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF] where it has been initiated. It also represents a good compromise, shown by the new Γ -convergence result, between our local model and the actual numerical model we implement which is introduced and motivated later, and falls within the "true nonlocal algorithms". In practice, for "true nonlocal methods", the computations are re-stricted to a small area around the point of interest (like the NL-means algorithm [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]) and it makes these methods quite similar to our model from our point of view. We now introduce the related nonlocal problem whose construction relies on a derivative free nonlocal formulation of the L 2 -norms R 2 |D (2,0) u| 2 dx and R 2 |D (0,2) u| 2 dx respectively, and more precisely on the following theorem ([18, Theorem 10]).

Theorem 1 Let u ∈ W 2,2 (R 2 ). Then R 2 R |u(x + he i ) -2u(x) + u(x -he i )| 2 |h| 4 ρ n (h) dh dx -→ n→+∞ D (2,0) u 2 L 2 (R 2 ) if i = 1, D (0,2) u 2 L 2 (R 2 ) if i = 2
.

Equipped with this theoretical result and characterization, we are henceforth able to derive a nonlocal counterpart of our local problem [START_REF] Adam | Denoising by inpainting[END_REF], phrased in terms of the functional F ε,n defined below and for which an existence result is provided (see [START_REF] Debroux | A second order free discontinuity model for bituminous surfacing crack recovery and analysis of a nonlocal version of it[END_REF]Theorem 12]).

Theorem 2 Let Ω be a regular bounded open subset of R 2 with boundary of class C 2 . Let us assume that the functions t → ρ n (t), t → t q ρ n (t) are non-increasing for t ≥ 0 and q ∈]0, 1[. (Such a function ρ n exists: for instance, with q ∈]0, 1[, ρ(t) = e -|t| |t| q and ρ n (t

) = C n ρ(nt) with C = 1 R ρ(t) dt ). With κ ε , ξ ε , ζ ε > 0, for any n ∈ N * , inf F ε,n (u, g, Q, v 1 , v 2 ) = f -u -div g 2 L 2 (Ω) + µ |∇Q| L ∞ (Ω) + γ 2 |g -∇Q| 2 L 2 (Ω) + ρ R 2 (v 2 2,e (x) + κ ε ) (2) 2 i=1 R |u e (x + he i ) -2u e (x) + u e (x -he i )| 2 |h| 4 ρn(h) dh dx + ξ ε Ω (v 2 1 + ζ ) |∇u| 2 dx + (α -β) G ε (v 1 ) + β G ε (v 2 ),
where v 2,e and u e are respectively the extensions of v 2 according to [12, Theorem IX.7, p. 158] -by construction, 0 ≤ v 2,e ≤ 1 a.e. -and of u on R 2 by 0 (-with the regularity assumed on Ω, v 2,e and u e are in W 1,2 (R 2 ) -), admits minimizers

(u n = u ε,n , g n = g ε,n , Q n = Q ε,n , v 1,n = v 1,ε,n , v 2,n = v 2,ε,n ) on W 1,2 0 (Ω)∩W s,2 (Ω)× H(div) × Q ∈ W 1,∞ (Ω) | Ω Q dx = 0 ×W 1,2 (Ω, [0, 1]) ×W 1,2 (Ω, [0, 1]), with s ∈ 3 2 , 2 .

Main theoretical result

Given the above material, we can now state the main theoretical result, a Γ -convergence one, relating the nonlocal minimization problem to the local one as n tends to +∞.

Theorem 3 (Γ -convergence) Let (u n , g n , Q n , v 1,n , v 2,n ) ∈ H 1 0 (Ω)∩W s,2 (Ω)×H(div)× {Q ∈ W 1,∞ (Ω)| Ω Q dx = 0} × W 1,2 (Ω, [0, 1]) ×W 1,2 (Ω, [0, 1]) with s ∈ [ 3 2
, 2[ be a sequence of minimizers of (2) for each n ∈ N * . Let us assume additionally for technical purposes, that

v 2,n ∈ {v 2 ∈ W 1,2 (Ω, [0, 1])} ∩ W 1,∞ (Ω) with sup n∈N * ∇v 2,n L ∞ (Ω) ≤ C 1 < ∞. Then there exist a sub- sequence still denoted by (u n , g n , Q n , v 1,n , v 2,n ) and a minimizer (ū, ḡ, Q, v1 , v2 ) ∈ W 2,2 (Ω)∩H 1 0 (Ω)×H(div) ×{Q ∈ W 1,∞ (Ω)| Ω Q dx = 0} × W 1,2 (Ω, [0, 1]) ×W 1,2 (Ω, [0, 1]) of (1) such that u n n→+∞ ū in W 3 2 ,2 (Ω), g n n→+∞ ḡ in H(div), Q n * n→+∞ Q in W 1,∞ (Ω), v 1,n n→+∞ v1 in W 1,2 (Ω), v 2,n n→+∞ v2 in W 1,2 (Ω) and Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) -→ n→+∞ Fε (ū, ḡ, Q, v1 , v2 ). Proof
The proof is subdivided into three main steps that we sketch hereafter. The details can be found in Appendix A.

Step 1. We have proved in what precedes that for any n ∈ N * , there exists a solution to [START_REF] Ambrosio | Variational Approximation of a Second Order Free Discontinuity Problem in Computer Vision[END_REF]. We consider a sequence of such minimizers (u n ,

g n , Q n , v 1,n , v 2,n ) ∈ H s (Ω) ∩ H 1 0 (Ω) × H(div) × {Q ∈ W 1,∞ (Ω)| Ω Q dx = 0}×W 1,2 (Ω, [0, 1])×W
) weakly- * converges to Q in W 1,∞ (Ω) with Ω Q dx = 0, (v 1,n ) weakly converges to v1 in W 1,2 (Ω) with v1 ∈ W 1,2 (Ω, [0, 1]), (v 2,n ) weakly converges to v2 in W 1,2 (Ω) with v2 ∈ W 1,2 (Ω, [0, 1]), (g n ) weakly converges to ḡ in H(div), (u n ) weakly converges in W 1,2 (Ω) to ū with ū ∈ W 1,2
0 (Ω), and (u n,e ) weakly converges in H 3 2 (R 2 ) to ūe . Some additional regularity is then obtained showing that ū ∈ W 2,2 (Ω) ∩ W 1,2 0 (Ω).

Step 2. In a second stage, we demonstrate that

lim sup n→+∞ F n,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤ F ε (ū, ḡ, Q, v1 , v2 ).
Step 3. The last step consists in showing that

F ε (ū, ḡ, Q, v1 , v2 ) ≤ lim inf n→+∞ F n,ε (u n , g n , Q n , v 1,n , v 2,n ), (ū, ḡ, Q, v1 , v2
) being a minimizer of (1).

Numerical implementation

Optimality conditions

We first derive the Euler-Lagrange equations according to each unknown (-for the sake of simplicity, we did not include the interpolation constraints-), with x = (x 1 , x 2 ) and n = (n x1 , n x2 ), the unit outward normal to the boundary. Making use of the absolutely minimizing Lipschitz extensions ( [START_REF] Aronsson | A Tour of the Theory of Absolutely Minimizing Functions[END_REF]) for the equation in Q, we get:

                     v 1 = α-β 2ε + 2(α -β)ε∆v 1 2ξ ε |∇u| 2 + α-β 2ε , v 2 = β 2ε +2βε∆v 2 2ρ 2 i=1 R |ue(x+he i )-2ue (x)+ue (x-he i )| 2 |h| 4 ρn(h) dh+ β 2ε
,

g 1 = ∂ x1 Q -2 γ ∂ x1 (f -u -div g) , g 2 = ∂ x2 Q -2 γ ∂ x2 (f -u -div g) , -µ∆ ∞ Q -γ∆Q + γdiv g = 0,
combined with suitable boundary conditions. Focusing now on the variable u, denoting by η a test function, by a real number, and setting

J( ) = f -u -η -div g 2 L 2 (Ω) + ρ R 2 (v 2 2,e (x) + κ ε ) 2 i=1 R |ue (x+he i )+ η(x+he i )-2ue(x)-2 η(x)+ue (x-he i )+ η(x-he i )| 2 |h| 4 ρ n (h) dh dx + ξ ε Ω (v 2 1 + ζ ) |∇u + ∇η| 2 dx,
we obtain, with suitable boundary conditions:

J ( ) = -2 Ω (f -u-η-divg)η dx + ρ R 2 (v 2 2,e (x) + κ ε ) 2 i=1 R 2(ue(x+he i )+ η(x+he i )-2ue (x)-2 η(x)+ue (x-he i )+ η(x-he i )) |h| 4 (η(x+hei)-2η(x)+η(x-hei))ρn(h) dh dx + 2ξ ε Ω (v 2 1 + ζ ) ∇u + ∇η, ∇η dx, J (0) = - Ω 2(f -u-divg)η dx + ρ R 2 (v 2 2,e (x) + κ ε ) 2 i=1 R 2(ue(x+he i )-2ue (x)+ue (x-he i ))(η(x+he i )-2η(x)+η(x-he i )) |h| 4 ρ n (h) dh dx + 2ξ ε Ω (v 2 1 + ζ ) ∇u, ∇η dx = 0, J (0) = - Ω (f -u-divg)η dx -ξ ε Ω div((v 2 1 +ζ )∇u)η dx +ρ R 2 2 i=1 R -2(v 2 2,e (x)+κε)(ue(x+he i )-2ue(x)+ue(x-he i )) |h| 4 + (v 2 2,e (x-he i )+κε)(ue(x)-2ue(x-he i )+ue (x-2he i )) |h| 4 + (v 2 2,e (x+he i )+κε)(ue (x+2he i )-2ue(x+he i )+ue (x)) |h| 4
ηρn(h) dh dx = 0, which yields the desired equation in u. The equation in Q encodes the most technical elements (nonlinearity, singularity at ∇Q = 0). The next subsection is thus devoted to the derivation and analysis of the evolution equation satisfied by Q and involving the infinity Laplacian, in the viscosity solution theory framework. This theory applies to certain partial differential equations and allows merely continuous functions to be solutions of fully nonlinear equations. We refer the reader to [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF] for a general introduction including the main material and definitions.

3.2 Theoretical results: existence, uniqueness and regularity of the solution of the evolution equation in Q

In that purpose, let us introduce an artificial time and embed the equation in Q in a time-dependent setting.

Let T > 0 be given. The parabolic equation in the unknown Q (which is the one we solve in practice) is thus given by For the sake of conciseness and using the normalized version of the infinity Laplacian, the evolution equation is now written in the form

∂Q ∂t = µ∆ ∞ Q + γ∆Q -γdiv g on R 2 × (0, T ), Q(x, 0) = Q 0 (x) on R 2 , (EE) with Q 0 ∈ W 1,∞ (R
∂Q ∂t + G(x, t, ∇Q, ∇ 2 Q) = 0, with G : R 2 × [0, T ) × R 2 × S 2 (S 2
being the set of symmetric 2 × 2 matrices equipped with its Loewner partial order) defined by

G(x, t, p, X) = -γtrace(X) -µ p |p| , X p |p| + γdiv g = -γtrace(X) -µtrace p ⊗ p |p| 2 X + γdiv g, = E(X) + F (p, X) + γdiv g,
⊗ representing the Kronecker product, and with the following properties :

1. The operators G, E : X → -γtrace(X) and F :

(p, X) → -µ trace p⊗p |p| 2 X are independent of Q and are elliptic, i.e., ∀X, Y ∈ S 2 , ∀p ∈ R 2 \ {0 R 2 }, if X Y then F (p, X) ≥ F (p, Y ) since F (p, X) - F (p, Y ) = -µtrace p⊗p |p| 2 X + µtrace p⊗p |p| 2 Y = -µtrace p⊗p |p| 2 (X -Y ) = -µ p |p| , (X -Y ) p |p| ≥ 0 as X Y . 2. F is locally bounded on R 2 × S 2 , continuous on R 2 \ {0 R 2 } × S 2
, and F * (0, 0) = F * (0, 0) = 0, where F * (resp. F * ) is the upper semicontinuous (usc) envelope (resp. lower semicontinuous (lsc) envelope) of F . Indeed, using Rayleigh quotient and its properties, it is not difficult to see that for nonzero vector p,

λ min (X) ≤ trace p⊗p |p| 2 X = p, Xp p, p = R(X, p) ≤ λ max (X), λ min (resp. λ max ) denoting the smallest (resp. biggest) eigenvalue of X.
We begin our analysis by establishing a comparison principle stating that if a sub-solution and a supersolution are ordered at the initial time, then they are ordered at any time.

Theorem 4 (Comparison principle, adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]) Assume (H) and let q -: R 2 × [0, T ) → R be a bounded upper semicontinuous sub-solution and q + : R 2 × [0, T ) → R be a bounded lower semicontinuous super-solution of (EE). Assume that q -(x, 0)

≤ Q 0 (x) ≤ q + (x, 0) in R 2 , then q -≤ q + in R 2 × [0, T ).
Proof We refer the reader to Appendix B.

We now construct barriers in order to use Perron's method.

Theorem 5 (Construction of barriers, adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]) Assume (H) and let

Q 0 ∈ W 1,∞ (R 2 ). Then q + = sup x∈R 2 Q 0 (x) + γ div g L ∞ (R 2 ×[0,T )) t and q -= inf x∈R 2 Q 0 (x)-γ div g L ∞ (R 2 ×[0,T ))
t are respectively super-and sub-solution of (EE).

Proof The proof of this theorem is straightforward. Indeed, since Q 0 ∈ W 1,∞ (R 2 ), then q -and q + are twice differentiable in space and once differentiable in time and are bounded. Besides q

-(x, 0) = inf x∈R 2 Q 0 (x) ≤ Q 0 (x), q + (x, 0) = sup x∈R 2 Q 0 (x) ≥ Q 0 (x), ∀x ∈ R 2 . Furthermore, ∀(x, t) ∈ R 2 × [0, T ), q - t (x, t) = -γ div g L ∞ (R 2 ×[0,T )) , Dq -(x, t) = 0 and D 2 q -(x, t) = 0; q + t (x, t) = γ div g L ∞ (R 2 ×[0,T )) , Dq + (x, t) = 0 and D 2 q + (x, t) = 0. It follows then that ∀(x, t) ∈ R 2 × [0, T ), q - t (x, t) + G * (t, x, q -(x, t), Dq -(x, t), D 2 q -(x, t)) = -γ div g L ∞ (R 2 ×[0,T )) + γdivg(x, t) ≤ 0,
and

q + t (x, t) + G * (t, x, q + (x, t), Dq + (x, t), D 2 q + (x, t)) = γ div g L ∞ (R 2 ×[0,T )) + γdivg(x, t) ≥ 0.
Hence q -is a sub-solution and q + a super-solution of (EE).

We are now able to use Perron's method and to give the following existence result.

Theorem 6 (Existence and uniqueness of a solution) Assume (H) and

Q 0 ∈ W 1,∞ (R 2 ).
Then there exists a unique bounded continuous solution of (EE) in

R 2 × [0, T ).
Proof We follow Perron's method ([15, Theorem 4.1]). Indeed, we have constructed a bounded sub-solution and a bounded super-solution which fall within the comparison principle.

Then function Q = sup{w, w sub-solution of (EE) and

∀(x, t) ∈ R 2 × [0, T ), q -(x, t) ≤ w(x, t) ≤ q + (x, t)} is a solution of (EE) potentially discontinuous.
Clearly, Q is bounded since q + and q -are bounded and since Q is a solution of (EE), then Q * is a super-solution and Q * is a sub-solution. Applying the comparison principle yields

Q * ≤ Q * . But by definition, Q * ≤ Q * , so we finally get Q = Q * = Q * meaning that Q is continuous on R 2 × [0, T ).
We now analyze the regularity of the solution. Let us first consider the regularity in space.

Theorem 7 (Lipschitz regularity in space, adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]) Assume (H) and that ∇Q 0 L ∞ (R 2 ) ≤ B 0 with B 0 > 0. Then the solution of (EE) is Lipschitz continuous in space and satisfies ∇Q(., t)

L ∞ (R 2 ) ≤ B(t), with B(t) = γκ g t + B 0 .
Proof We refer the reader to Appendix C.

Besides, by analyzing the smoothness in time of the solution, we can show that this solution is uniformly continuous in time.

Theorem 8 (adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]) Assume (H), and that div g is uniformly continuous in time with ω div g its modulus of continuity. Then the solution of (EE) is uniformly continuous in time.

Proof We refer the reader to Appendix D.

This concludes the analysis of the evolution equation in Q in the viscosity solution theory setting.

Handling of the nonlocal component and resulting algorithm:

Some limitations can be identified regarding the theoretical expression of the nonlocal component, which leads us to introduce a numerical compromise. In particular, the properties of the kernels ρ n are not all welladapted for imaging problems. We recall that they sa-

tisfy ∀n ∈ N, ∀x ∈ R, ρ n (x) = ρ n (|x|); ∀n ∈ N, ρ n ≥ 0; ∀n ∈ N, R ρ n (x) dx = 1; ∀δ > 0, lim n→+∞ +∞ δ
ρ n (r)dr = 0). While the last three properties are consistent with the discrete setting as in practice the weights are nonnegative, normalized and concentrated near the current window center, the radiality property is not relevant for imaging problems. In fact, we would like to define a nonlocal version of our second order operator at a given point x that grants more weight to points that belong to the same region as x not only based on the difference between intensities but for instance, on the difference between patches around those points and thus not favoring spatial proximity. This first observation led us to reconsider the definition of our nonlocal weights. We would like also to stress that, to the best of our knowledge, this assumption of radiality cannot be removed in the theoretical analysis (see [START_REF] Bourgain | Another look at Sobolev spaces[END_REF], [START_REF] Dávila | On an open question about functions of bounded variation[END_REF] and [START_REF] Ponce | A new approach to Sobolev spaces and connections to Γ -convergence[END_REF]). We thus redefine the weights ρn(h) |h| 4 involved in the theoretical expression, inspired by the NL-means algorithm ( [START_REF] Buades | A Review of Image Denoising Algorithms, with a New One[END_REF]). Integrating additional information related to the image content f is pertinent here since it strengthens/favors similar behaviors. More precisely, we put more weight to neighboring pixels that have similar edges/ creases identified by the gradient / Hessian behaviour as seen in Lemma 1 and that are geographically close. In that purpose, we consider the following nonlocal weights ρn(h)

|h| 4 ≈ w f ,i,x (h) = exp - d f ,x,i (h) α 2 where d f ,x,i (h) = R 2 G a (t) ∇f (x+t)-∇f (x+t+he i ) 2 dt or d f ,x,i (h) = R 2 G a (t) ∇ 2 f (x+t)-∇ 2 f (x+t+he i ) 2 dt, G a being a
Gaussian kernel with standard deviation a controlling the patch size, and α is the filtering parameter. The data structure summed-area table also known as integral image is used to efficiently compute the involved sums. Later in the text is a pseudo-code associated with the computation of the nonlocal weights (Algorithm 1).

The resulting nonlocal algorithm relies on an alternating strategy in which we solve the Euler-Lagrange equations related to each unknown using specifically, a timedependent scheme for Q = Q(x 1 , x 2 , t) (nonlinear overrelaxation method, see [14, Section 4]), and a stationary semi-implicit fixed-point scheme in

v 1 = v 1 (x 1 , x 2 ), v 2 = v 2 (x 1 , x 2
) and g = g(x 1 , x 2 ). Neumann boundary conditions for u and Q, and Dirichlet boundary conditions for v 1 -1, v 2 -1 and g are applied. We denote by ∆x 1 = ∆x 2 = h = 1 the space steps in each direction, by ∆t the time step, and by for all pixels x = (x 1 , x 2 ) do 3. Compute the distances (we do not make use in practice of the Gaussian kernel G a , yielding equal weights for the contributions)

f i,j , u n i,j , Q n i,j , v n 1,i,j , v n 2,i,j , g n i,j = (g n 1,i,j , g n 2,i,j ) T the discrete versions of f , u, Q, v 1 , v 2 and g at iteration n ≥ 0, for 1 ≤ i ≤ M , 1 ≤ j ≤ N .
d f ,x,1 (y) = p-1 2 i=- p-1 2 p-1 2 j=- p-1 2 ∇f (x 1 +i,x 2 +j)-∇f (x 1 +y+i,x 2 +j) 2 2 + ∇f (x 1 +i,x 2 +j)-∇f (x 1 +i-y,x 2 +j) 2 2 or d f ,x,1 (y) = p-1 2 i=- p-1 2 p-1 2 j=- p-1 2 ∇ 2 f (x 1 +i,x 2 +j)-∇ 2 f (x 1 +y+i,x 2 +j) 2 2 + ∇ 2 f (x 1 +i,x 2 +j)-∇ 2 f (x 1 +i-y,x 2 +j) 2 2
between all patches centered at x + ye 1 of size p with 0 ≤ y ≤ iw-1 2 and the patch centered at the current x. Compute the distance

d f ,x,2 (y) = p-1 2 i=- p-1 2 p-1 2 j=- p-1 2 ∇f (x 1 +i,x 2 +j)-∇f (x 1 +i,x 2 +y+j) 2 2 + ∇f (x 1 +i,x 2 +j)-∇f (x 1 +i,x 2 +j-y) 2 2 .
or

d f ,x,2 (y) = p-1 2 i=- p-1 2 p-1 2 j=- p-1 2 ∇ 2 f (x 1 +i,x 2 +j)-∇ 2 f (x 1 +i,x 2 +y+j) 2 2 + ∇ 2 f (x 1 +i,x 2 +j)-∇ 2 f (x 1 +i,x 2 +j-y) 2 2 .
between all patches centered at x + ye 2 of size p with 0 ≤ y ≤ iw-1 2 and the patch centered at the current x. end for for all pixels x do 4.1. Sort the previous distances in ascending order for each direction and keep only the lowest N bN eigh-1 values with the corresponding shift (y) for the first direction in indice 1 and for the second direction in indice 2 . 4.2. Add the closest neighbor in geographical sense and in each direction to make the weights more similar to the theoretical ones. end for for all pixels x do 5.1. Compute w f ,x,1 by the following formula: w f ,x,1 (y) = 0 if y does not belong to the previous list of neighbors, w f ,x,1 (y) = exp -

d f ,x,1 (y) h 2
otherwise. Compute w f ,x,2 by the following formula: w f ,x,2 (y) = 0 if y does not belong to the previous list of neighbors, w f ,x,2 (y) = exp -

d f ,x,2 (y) h 2 otherwise. end for Return w f ,x,1 , w f ,x,2 , indice 1 , indice 2 .
Algorithm 2 Nonlocal alternating algorithm for crack recovery.

[Initialization step]:

u 0 = f , g 0 = 0, Q 0 = 0, v 0 1 = 1 and v 0 2 = 1. 2. [Main step]: for n=1 to n=500 do if n%100==0 then 2.1 Compute the nonlocal weights associated to u n : [w f ,x,1 , w f ,x,2 , indice 1 , indice 2 ] = compute weights(u n ). end if 2.
2 Compute for all pixels (i, j):

|∇u n | 2 i,j = (u n i+1,j -u n i,j ) 2 + (u n i,j+1 -u n i,j ) 2 , v n+1 1,i,j = (α-β) 2ε +2(α-β)ε (v n 1,i+1,j +v n 1,i-1,j +v n 1,i,j+1 +v n 1,i,j-1 -4v n 1,i,j ) 2ξ ε |∇u n | 2 i,j + (α-β) 2ε , nlnormu 2 i,j = h∈indice 1 (u n (i + h, j) -2u n (i, j) + u n (i -h, j)) 2 w f ,i,j,1 (h)) + h∈indice 2 (u n (i, j + h) -2u n (i, j) + u n (i, j -h)) 2 w f ,i,j,2 (h), v n+1 2,i,j = β 2ε +2βε (v n 2,i+1,j +v n 2,i-1,j +v n 2,i,j+1 +v n 2,i,j-1 -4v n 2,i,j ) 2ρ nlnormu 2 i,j + β 2ε , u n+1 i,j -u n i,j ∆t = (f i,j -u n i,j - g n 1,i+1,j -g n 1,i-1,j 2 - g n 2,i,j+1 -g n 2,i,j-1 2 ) + ξ ε (v n 1,i,j ) 2 (u n i,j+1 -u n i,j ) -(v n 1,i,j-1 ) 2 (u n i,j -u n i,j-1 ) + ξ ε (v n 1,i,j ) 2 (u n i+1,j -u n i,j ) -(v n 1,i-1,j ) 2 (u n i,j -u n i-1,j ) -h∈indice 1 ρw f ,i,j,1 (h) (v n 2,i+h,j ) 2 (u n i+2h,j -2u n i+h,j +u n i,j ) -2 (v n 2,i,j ) 2 (u n i+h,j -2u n i,j +u n i-h,j )+(v n 2,i-h,j ) 2 (u n i,j -2u n i-h,j +u n i-2h,j ) -h∈indice 2 ρw f ,i,j,2 (h) (v n 2,i,j-h ) 2 (u n i,j -2u n i,j-h +u n i,j-2h ) +(v n 2,i,j+h ) 2 (u n i,j+2h -2u n i,j+h +u n i,j )-2 (v n 2,i,j ) 2 (u n i,j+h -2u n i,j +u n i,j-h ) ,
using symmetry if it does not belong to the image domain, and equations derived in the same way for g n 1 , g n 2 and Q n . end for Remark 1 Provided that we solve exactly each subproblem, a result of convergence of the alternating strategy to existing partial minimisers can be proved.

In order to improve the computation efficiency, we propose an MPI parallelization of our code, which motivates our choice of a rather simple alternating minimization method for which a decomposition domain approach is well-suited.

MPI parallelization

The parallelization of the C code is motivated by the natural geometry of the problem -an image is defined on a rectangle domain Ω -making the partition of the image domain into subdomains supporting simultaneous local computations relevant. Note also that the computational complexity increases with the image size (in practice we have worked with some images of size 2248 × 4000), requiring more memory to store the data and the results, this fact being particularly marked in the nonlocal model that involves the resolution of a nonlocal partial differential equation.

The meshing is made of ntx interior points in the row direction (we removed the first and last rows of points) and nty interior points in the column direction (we removed the leftmost and rightmost columns of points). The implementation revolves around the following steps:

(i) we generate a Cartesian topology (see Figure 1 for an example), each subdomain comprising w-1 (-w is the window size -) rows of ghost cells above, w -1 rows of ghost cells below, and similarly for the columns, in order to store the data exchanged with neighboring subdomains. Either the developer selects the number of nodes in each direction, or it is left to the MPI library. Some latitude is also given to the user in terms of periodicity (-a periodicity can be applied on the grid in each direction if required thanks to the array periods -) and reorganization (-if the user wants the processes to keep the same rank as in the original communicator -). (ii) For each subdomain, we recover the bounds with respect to the original image reference frame of the indices i and j that are then stored in the 1d array tab bounds: tab bounds[0]=sx, tab bounds [START_REF] Adam | Denoising by inpainting[END_REF]=ex, tab bounds [START_REF] Ambrosio | Variational Approximation of a Second Order Free Discontinuity Problem in Computer Vision[END_REF]=sy and tab bounds [START_REF] Aronsson | A Tour of the Theory of Absolutely Minimizing Functions[END_REF]=ey. The numbering of the original image reference frame starts at 0 and the origin is the top left corner. The indices (sx -w + 1, . . . , sx -1, ex + 1, . . . ,ex + w -1, sy -w + 1, . . . , sy -1, ey + 1, . . . , ey + w -1) are used to store the data sent by the 8-connected neighboring subdomains. Created function : void domaine(MPI Comm comm2d, int rang, int * coords, int ntx, int nty, int * dims, int * tab bounds, int * periods, int reorganization, int nb procs) (iii) For each subdomain, the neighboring subdomain ranks are returned. This is achieved thanks to the routine voisinage and these ranks are stored respectively in the 1d array voisin (for the 4-connected blocks) and voisin diagonale (for the diagonally connected blocks). The necessity of storing the ranks of diagonally connected blocks arises from the numerical schemes used to discretize the partial differential equations satisfied by u, Q, g 1 and g 2 that involve for instance components like Q i-1,j+1 (vii) The values of

u n+1 , Q n+1 , v n+1 1 , v n+1 2
, g n+1 1 and g n+1 2 are computed using the above mentioned finite difference schemes (v 1 and v 2 have been initialized to 1, g 1 , g 2 and Q to 0, while u has been set to the values of the original image) and the question of boundary conditions is addressed. For the sake of simplicity, we have assumed homogeneous Neumann boundary conditions for u and Q, simulating reflection of the array through its boundaries, and Dirichlet boundary conditions for v 1 -1, v 2 -1 and g. We thus identify subdomains with at least one edge included in the image domain boundary. As an illustration, for the variable u, we replicate the row of index iw-1 (resp. ex-sx+iw-1) of the local matrix in the row of index iw -2 (resp. ex -sx + iw), similarly for the columns. The diagonal components are processed using point symmetry. We no longer comment on the optimization of the computation code but now focus on the quantitative and qualitative analysis of the results.

Numerical experiments

Before showing some numerical simulations on real datasets resulting from the application of the above algorithm (with different weights) and discussing its ac-curacy, we would like to give some insight on how to choose the parameters involved in our model. The more sensitive ones are α, β, ρ. Smaller values for these parameters induce less regularization and consequently more information (edges/creases) in the v 1 and v 2 components. Parameter µ weights the L ∞ -norm of |∇Q| and thus influences the size of the texture captured in v = divg. A higher parameter µ leads to smaller scale features caught by v = divg. ε plays on the thickness of the contours in v 1 and v 2 : larger values of ε imply thicker contours.

We now apply the proposed algorithm to crack detection on concrete walls, both on Figure 3 [START_REF] Drogoul | Numerical Analysis of the Topological Gradient Method for Fourth Order Models and Applications to the Detection of Fine Structures in Imaging[END_REF]). The cracks are correctly enhanced while the oscillatory patterns are well captured by the v = div g component. We also observe that the crack is smoother with Aubert and Drogoul's method and our method seems to detect more accurately the center of the crack. The contrast is the highest for the nonlocal method with weights based on Hessian information and the connectivity of the crack appears clearer too. Again, the role of the decomposition part of the algorithm is highlighted (Figure 3) by depicting the obtained v 2 component when decomposition is turned off in our local implementation (spurious details are visible on the top of the image). Also, the piecewise linear nature of the component u in Figure 4 is properly returned. We conclude the paper with two applications dedicated to crack detection on bituminous surfacing Figure 5 (size 231 × 650) and 6 (size 201 × 640), courtesy of CEREMA (Centre d' Études et d'Expertise sur les Risques, l'Environnement et l'Aménagement), France, that motivated our work. The two considered slices of bitumen, in addition to long and thin cracks, exhibit high oscillatory patterns and white spots of varying sizes, which makes the straight application of our algorithm difficult. Indeed, in terms of scale, the crack and some of these spots could be comparable and could not be properly discriminated, resulting in superfluous information in the v 2 component. Think for instance of a white spot assimilated to a ball of radius 2 pixels (-if the image domain is the n × n discretized unit square, then the scale behaves like 1 n -), and of a long thin crack of width 2 pixels and length k pixels (k 1) leading to a similar scale. To circumvent this issue, a pre-processing step is applied. It consists in apprehending the problem first as an inpainting one ( [START_REF] Adam | Denoising by inpainting[END_REF]), and by considering these white spots as missing parts of the image that need to be filled. This is achieved with the MATLAB R function imfill (https://fr.mathworks. com/help/images/ref/imfill.html -to fill holes in a grayscale image) applied to the inverse image, yielding an image that serves as input of our algorithm. In both cases, the cracks are well recovered in the v 2 component which does not include superfluous information. The edge detector v 1 also recovers parts of the crack but contains spurious information regarding the problem we address, such as asphalt defect boundaries. It thus justifies the use of a second order method. Besides, Figure 5-( oscillatory texture of the images and to discriminate it from the crack. To further improve the quality of the detection, a nonlocal version of the model has been theoretically studied with a new Γ -convergence result, and numerically examined by providing a detailed algorithmic approach and an efficient MPI implementation.

The method has been tested on real images given by the CEREMA and shows good performance with still room for improvement. A first lead would consist in diffusing along the obtained cracks in order to highlight the contrast and fill the missing data.

Appendix A -Proof of Theorem 3: Γ -convergence result

Step 1.

We have proved in what precedes that n ∈ N * , there exists a solution to [START_REF] Ambrosio | Variational Approximation of a Second Order Free Discontinuity Problem in Computer Vision[END_REF]. Let us consider a sequence of such minimizers (u n ,

g n , Q n , v 1,n , v 2,n ) ∈ H s (Ω) ∩ H 1 0 (Ω)× H(div) × {Q ∈ W 1,∞ (Ω)| Ω Q dx = 0} × W 1,2 (Ω, [0, 1])× W 1,2 (Ω, [0, 1]), with s ∈ [ 3 2 , 2[. Let u ∈ W 2,2 (Ω)∩W 1,2 0 (Ω), g ∈ H(div), Q ∈ {Q ∈ W 1,∞ (Ω)| Ω Q dx = 0}, v 1 ∈ W 1,2 (Ω, [0, 1]), v 2 ∈ W 1,2 (Ω, [0, 1]), then ∀n ∈ N * , Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤ Fn,ε (u, g, Q, v 1 , v 2 ) ≤ ρ(1+κ ε ) u 2 W 2,2 (Ω) +ξ ε Ω (v 2 1 (x)+ ζ ε )|∇u(x)| 2 dx + f -u -divg 2 L 2 (Ω) + µ ∇Q L ∞ (Ω) + γ 2 g -∇Q 2 L 2 (Ω) + (α -β) Ω (v 1 (x)-1) 2 4ε + ε|∇v 1 (x)| 2 dx + β Ω (v 2 (x)-1) 2 4ε +ε|∇v 2 (x)| 2 dx = C < +∞ using [18, Lemma 2].
We deduce that: -∇Q n L ∞ (Ω) is uniformly bounded with respect to n and since Ω Q n dx = 0 for any n ∈ N * , then using Poincaré-Wirtinger inequality, we get that (Q n ) is uniformly bounded in W 1,∞ (Ω). Thus we can extract a subsequence still denoted by

(Q n ) weakly- * converging to Q in W 1,∞ (Ω). Besides, the compact embedding W 1,∞ (Ω) c C 0 ( Ω) leads to Ω Q dx = 0. -∇v 1,n L 2 (Ω)
is uniformly bounded with respect to n and since v 1,n ∈ W 1,2 (Ω, [0, 1]) for any n ∈ N * and Ω being bounded, then by using Poincaré-Wirtinger inequality, we get that (v 1,n ) is uniformly bounded in W 1,2 (Ω). We can therefore extract a subsequence still denoted by (v 1,n ) weakly converging to v1 in W 1,2 (Ω) and thanks to the compact embedding W 1,2 (Ω) c L 2 (Ω), the convergence is also pointwise almost everywhere up to a subsequence, leading to v1 ∈ W 1,2 (Ω, [0, 1]). -Using the same arguments, we get that v 2,n n→+∞ v2 in

W 1,2 (Ω) with v2 ∈ W 1,2 (Ω, [0, 1]).
-∇u n L 2 (Ω) is uniformly bounded with respect to n and since u n ∈ W 1,2 0 (Ω) for any n ∈ N * , then using Poincaré inequality, we can extract a subsequence still denoted by (u n ) weakly converging in W 1,2 (Ω) to ū. By continuity of the trace operator, we have ū

∈ W 1,2 0 (Ω). -C ≥ 1 4 divg n 2 L 2 (Ω) -1 2 u n 2 L 2 (Ω) -f 2 L 2 (Ω) + γ 4 g n 2 L 2 (Ω) -γ 2 ∇Q n 2 L 2 (Ω)
. So, (g n ) is uniformly bounded in H(div) and we can extract a subsequence still denoted by (g n ) weakly converging to ḡ in H(div). Now, let us show that (u n,e ) weakly converges in H 3 2 (R 2 ) to ūe .

We set E n (h) = R 2 |u n,e (x + he 1 ) -2u n,e (x) + u n,e (xhe 1 )| 2 dx. It satisfies E n (2h) ≤ 16E n (h), ∀h ∈ R. By using Fubini-Tonelli theorem, we have

R 2 R |u n,e (x+he 1 )-2u n,e (x)+u n,e (x-he 1 )| 2 |h| 4 ρ n (h) dh dx = R E n (h) |h| 4 ρ n (h) dh = 2 ∞ 0 E n (h) |h| 4 ρ n (h) dh ≤ C. We then apply [5, Lemma 3.1] with M = δ = 1, g(t) = E n (t)
t q+1 , k(t) = t q-3 ρ n (t) and we get C(1)

1 0 t q-3 ρ n (t) dt 1 0 E n (t) t q+1 dt ≤ 1 0 E n (t)
t 4 ρ n (t) dt. (We will see further that the condition of monotonicity on k is fulfilled). We now need g to verify the assumption of this lemma, that is to say, g( t 2 ) ≥ g(t). We know that g(

t 2 ) = E n ( t 2 )
t q+1 2 q+1 ≥ 2 q-3 g(t). Thus if q ≥ 3, this condition is fulfilled. By using the properties of ρ n , we deduce first that 1 0 E n (t) t q+1 dt ≤ C with C independent of n for q = 3 and n large enough since then

1 0 t q-3 ρ n (t) dt = 1 0 ρ n (t) dt = 1 -∞ 1 ρ n (t) dt, ∀n ∈ N * with lim n→+∞ ∞ 1 ρ n (t) dt = 0 and so for n large enough, 1 0 ρ n (t) dt ∈ 1 2 , 1 . Then ∞ 1 E n (t) t q+1 dt ≤ C u n,e 2 L 2 (R 2 )
∞ 1 1 t q+1 dt, C being a constant and the last integral being convergent since q = 3, resulting in the uniform boundedness of ∞ 0 E n (t)

t q+1 dt. Besides, R E n (t) |t| q+1 dt = R 1 |h| q+1 R 2 |u n,e (x+he 1 )-2u n,e (x)+u n,e (x-he 1 )| 2 ρ n (h) dx dh = R 1 |h| q+1 τ he 1 u n,e -2u n,e + τ -he 1 u n,e 2 L 2 (R 2 ) dh = R 1 |h| q+1 R 2 |e 2iπhξ 1 -2 + e -2iπhξ 1 | 2 |F (u n,e )| 2 (
ξ) dξ dh by Plancherel's theorem (τ denoting the usual translation operator and ξ = (ξ 1 , ξ 2 )). Then one can prove that R E n (t)

|t| q+1 dt = C R 2 |F (u n,e )(ξ)| 2 |ξ 1 | q R sin 4 (u)
|u| q+1 du dξ ≤ C (the constant C hereafter may change line to line). The generalized integral in u converges if and only if q ∈ [3, 4[ so for q = 3. By using the same arguments in the other direction (e 2 ), we get that | • | q 2 F (u n,e )(•) ∈ L 2 (R 2 ) and so u n,e ∈ H 3 2 (R 2 ) (being a Hilbert space) and is uniformly bounded for the associated norm for n large enough. There exists a subsequence still denoted by (u n,e ) weakly converging to ũ in H 3 2 (R 2 ). Besides, we know that u n,e = u n on Ω and D (1,0) u n,e = (D (1,0) u n ) e = D (1,0) u n on Ω, and D (0,1) u n,e = (D (0,1) u n ) e = D (0,1) u n on Ω. Thus, u n 2

W 3 2 ,2 (Ω) = u n 2 W 1,2 (Ω) + Ω Ω |∇u n (y) -∇u n (x)| 2 |x -y| 3 dx dy ≤ C + R 2 R 2 |∇u n,e (y)-∇u n,e (x)| 2 |x-y| 3
dx dy, with C independent of n. From [START_REF] Demengel | Functional Spaces for the Theory of Elliptic Partial Differential Equations[END_REF]Lemma 4.33,p. 200], we know that

R 2 R 2 |∇u n,e (x) -∇u n,e (y)| 2 |x -y| 3 dx dy < ∞ ⇔ R R 2 |∇u n,e (x) -∇u n,e (x + he 1 )| 2 |h| 2 dx dh < ∞ and R R 2 |∇u n,e (x) -∇u n,e (x + he 2 )| 2 |h| 2 dx dh < ∞. Let us now prove that R R 2 |∇u n,e (x) -∇u n,e (x + he 1 )| 2 |h| 2 dx dh < ∞ and R R 2 |∇u n,e (x) -∇u n,e (x + he 2 )| 2 |h| 2 dx dh < ∞ independen- tly of n. We have R 1 |h| 2 τ he1 (∇u n,e ) -∇u n,e 2 L 2 (R 2 ) dh = C R sin 2 (u) |u| 2 R 2 |ξ 1 | (|ξ 1 | 2 + |ξ 2 | 2 ) |F (u n,e )(ξ)| 2 dξ du ≤ C u n,e 2
H 3 2 (R 2 )
by using Plancherel's theorem and with C independent of n. By doing the same computations in the other direction, we prove that u n W ,2 (Ω) ⊂ W 1,2 (Ω) and the mapping T : W ,2 (Ω) → W 1,2 (Ω), u → T u = u is linear and continuous for the strong topology, so for the weak topology ([19, Theorem III.9, p. 39]). As W (Ω) and everywhere on Ω. Then ũ = ū on Ω and ũ = 0 on ∂Ω, by uniqueness of the weak limit. Now,

H 3 2 (R 2 ) L 2 (R 2 ) S (R 2 ) D (R 2 ) with continuous embeddings. ∀ϕ ∈ D(R 2 ), R 2 (u n,e -ũ)ϕ dx -→ n→+∞ 0 = Ω (u n -ū)ϕ dx -→ n→+∞ 0 + R 2 \Ω (u n,e -ũ)ϕ dx. Consequently, ∀ϕ ∈ D(R 2 ), R 2 \Ω ũϕdx = R 2 \ Ω ũϕdx = 0 since ũ ∈ H 3 2 (R 2 ) C(R 2 ). In particular, ∀ϕ ∈ D(R 2 \ Ω), R 2 \
Ω ũϕ dx = 0, meaning that ũ = 0 on R 2 \ Ω in the sense of distributions. Due to the continuity of ũ, we deduce that ũ = 0 everywhere on R 2 \ Ω and so ũ = ūe ∈ H

3 2 (R 2 ). Now, let us prove that ū ∈ H 2 (Ω)∩H 1 0 (Ω) and ūe ∈ H 2 (R 2 ). Since u n -→ n→+∞ ū in L 2 (Ω) then u n,e -→ n→+∞ ūe in L 2 (R 2 ). Let ϕ ∈ C ∞ 0 (R 2 )
. We denote by ūe and by u n,e the extensions by 0 of ū and u n for any n

∈ N * on R 2 . Since Ω is of class C 2 then ūe , u n,e ∈ W 1,2 (R 2 ), for any n ∈ N * . We have R 2 R ϕ(x + he 1 ) -2ϕ(x) + ϕ(x -he 1 ) h 2 ρ n (h) dh(u n,e (x)
+ (ū e (x) -u n,e (x))) dx , ≤ R 2 R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh(u n,e (x)) dx + R 2 R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh( ūe (x)-u n,e (x)) dx , ≤ R 2 R (ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 )) 2 h 4 ρ n (h) dh dx 1 2 . R 2 R
|u n,e (x) -ūe (x)| 2 ρ n (h) dh dx

1 2 + R 2 R
un,e (x+he 1 )-2un,e (x)+un,e (x-he 1 )

h 2 ρ n (h) dh(ϕ(x)) dx ,
owing to Hölder's inequality with respect to the measure ρ n (h) dh dx. Since ϕ ∈ C ∞ 0 (R 2 ) (justification of the integrability and the possibility to change the order of integration is postponed -in fact, it is not difficult to check that Fubini-Lebesgue theorem applies),

≤ D (2,0) ϕ L 2 (R 2 ) u n,e -ūe L 2 (R 2 ) + R 2 R (u n,e (x+he 1 )-2u n,e (x)+u n,e (x-he 1 )) 2 h 4 ρ n (h) dh dx 1 2 • R 2 R ϕ 2 (x)ρ n (h) dh dx 1 2 , ≤ D (2,0) ϕ L 2 (R 2 ) u n,e -ūe L 2 (R 2 ) + R 2 R
(un,e (x+he 1 )-2un,e (x)+un,e (x-he 1 )) 2

h 4 ρ n (h) dh dx 1 2 ϕ L 2 (R 2 ) , ≤ D (2,0) ϕ L 2 (R 2 ) u n,e -ūe L 2 (R 2 ) + √ C ϕ L 2 (R 2 ) . Since ϕ ∈ C ∞ 0 (R 2 ), then for almost every x ∈ R 2 , R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh ūe (x) -→ n→+∞ D (2,0) ϕ(x)ū e (x) thanks to [18, Theorem 9]. Besides, ∀n ∈ N * and almost every x ∈ R 2 , R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dhū e (x) ≤ |ū e (x)| • R |ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 )| |h| 2
ρ n (h) dh by applying Jensen's inequality with respect to the measure ρ n (h) dh (see [START_REF] Spector | Characterization of Sobolev and BV Spaces[END_REF]Theorem 139,p. 56

]). But ūe ∈ L 2 (R 2 ) with ūe L 2 (R 2 ) = ū L 2 (Ω) ≥ c ū L 1 (Ω) = c ūe L 1 (R 2 ) , since ū is extended by 0 outside Ω. As ϕ ∈ C ∞ 0 (R 2 ), then R |ϕ(x + he 1 ) -2ϕ(x) + ϕ(x -he 1 )| |h| 2 ρ n (h) dh ∈ L 1 (R 2 ) and R |ϕ(x + he 1 ) -2ϕ(x) + ϕ(x -he 1 )| |h| 2 ρ n (h) dh ≤ R 1 0 1 0 |D (2,0) ϕ(x + (t + s -1)he 1 )| ρ n (h) ds dt dh ≤ D (2,0) ϕ L ∞ (R 2 )
. We also have

R 2 R ϕ(x + he 1 ) -2ϕ(x) + ϕ(x -he 1 ) h 2 ρ n (h) dh 2 dx ≤ R 2 R (ϕ(x + he 1 ) -2ϕ(x) + ϕ(x -he 1 )) 2 h 4 ρ n (h) dh dx ≤ D (2,0) ϕ 2 L 2 (R 2 ) < +∞. Thus R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh ∈ L 2 (R 2 ), and ūe being L 2 (R 2 ), it leads to R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh ūe (x) ∈L 1 (R 2 ) with R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh ūe (x) ≤ |ū e (x)| D (2,0) ϕ L ∞ (R 2 ) ∈ L 1 (R 2 ). From the dominated con- vergence theorem, we can conclude that lim n→+∞ R 2 ūe (x) R ϕ(x+he 1 )-2ϕ(x)+ϕ(x-he 1 ) h 2 ρ n (h) dh dx = R 2 ūe (x)D (2,0) ϕ(x) dx.
By letting n tend to infinity in the previous inequality, we get R 2 ūe (x)D (2,0)

ϕ(x) dx ≤ √ C ϕ L 2 (R 2 ) . Eventually, D (2,0) ūe ∈ L 2 (R 2 ) using [12, Proposition VIII.3]
. By applying the same reasoning, we also get D (0,2) ūe ∈ L 2 (R 2 ). Since ūe ∈ L 2 (R 2 ) then it is a tempered distribution and so are its successive derivatives. We can now take the Fourier transform D (1,1) u e (ξ) = -ξ 1 ξ 2 4π 2 u e (ξ), D (0,2) u e (ξ) = -ξ 2 ξ 2 4π 2 u e (ξ), D (2,0) u e (ξ) = -ξ 1 ξ 1 4π 2 u e (ξ). Since D (0,2) u e ∈ L 2 (R 2 ) and D (2,0)

u e ∈ L 2 (R 2 ), then -ξ 2 ξ 2 u e (ξ) ∈ L 2 (R 2 ) and -ξ 1 ξ 1 u e (ξ) ∈ L 2 (R 2 ). But R 2 4π 2 ξ 2 1 ξ 2 2 u e 2 (ξ) dξ ≤ 2π 2 R 2 ξ 4 1 u e 2 (ξ) dξ + 2π 2 R 2 ξ 4 2 u e 2 (ξ) dξ < +∞,
which means that D (1,1) u e (ξ) ∈ L 2 (R 2 ) and by Plancherel's theorem we can conclude that D (1,1)

u e ∈ L 2 (R 2 ). This proves that ūe ∈ W 2,2 (R 2 ) since ūe ∈ W 1,2 (R 2 ) by construction.
As Ω ∈ C 2 and ūe is the extension of ū by 0 outside

Ω, then ū ∈ W 2,2 (Ω) ∩ W 1,2 0 (Ω).
Step 2.

By definition of the sequence (u

n , g n , Q n , v 1,n , v 2,n ), we have ∀n ∈ N * , Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤ Fn,ε (ū, ḡ, Q, v1 ,
v2 ) and by taking the lim sup when n → +∞, it yields lim sup

n→+∞ Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤ Fε (ū, ḡ, Q, v1 , v2 ). In- deed, thanks to [18, Theorem 9], we know that (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -→ n→+∞ (v 2
2,e (x) + κ ε )|D (2,0) u| 2 (x) everywhere in R 2 and for all u ∈ C 4 c (R 2 ). Without loss of generality, we assume that supp(u) ⊂ B(0, R) with R > 0. We now aim to prove that ∀ε > 0, ∃L = L(ε) > 1 (we believe that the confusion with the ε from the elliptic approximation is not possible) such that

sup n B(0,LR) c (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4
ρ n (h) dh dx ≤ ε. We have that

B(0,LR) c (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh dx≤ (1+κ ε ) B(0,LR) c R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh dx ≤ (1 + κ ε )4 (L -1) 4 R 4 u 2 L 2 (R and the conclusion follows. As R 2 (v 2 2,e (x) + κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx = B(0,LR) (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx + B(0,LR) c (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx, ≤ B(0,LR) (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx + ε. We know that (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -→ n→+∞ (v 2
2,e (x) + κ ε )|D (2,0) u| 2 (x) everywhere in B(0, LR) and

(v 2 2,e (x) + κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u| 2 (x) ≤ 2(1+κ ε ) D (2,0) u 2 L ∞ (R 2 ) ∈ L 1 (B(0, LR)).
Using the dominated convergence theorem, we get that

B(0,LR) (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx -→ n→+∞ 0 and by letting ε tend to 0, we conclude that R 2 (v 2 2,e (x)+κ ε ) R |u(x+he 1 )-2u(x)+u(x-he 1 )| 2 |h| 4 ρ n (h) dh -(v 2 2,e (x) + κ ε )|D (2,0) u(x)| 2 dx -→ n→+∞ 0. We now extend this result to u ∈ H 2 (R 2 ). Let ε > 0. By density, there exists v ε ∈ C ∞ 0 (R 2 ) such that D (2,0) u -D (2,0) v ε L 2 (R 2 ) ≤ ε. Let us set u n (x,h)=(v 2 2,e (x)+κ ε ) 1 2 u(x+he 1 )-2u(x)+u(x-he 1 ) h 2 ρ 1 2 n (h) ∈L 2 (R 2 ×R) and v n,ε (x,h)=(v 2 2,e (x)+κ ε ) 1 2 vε (x+he 1 )-2vε (x)+vε (x+he 1 ) h 2 ρ 1 2 n (h)∈L 2 (R 2 ×R).
We have

u n -v n,ε L 2 (R 2 ×R) ≤ (1 + κ ε ) 1 2 D (2,0) u -D (2,0) v ε L 2 (R 2 ) ≤ (1 + κ ε ) 1 2 ε,
yielding, thanks to the second triangle inequality,

u n L 2 (R 2 ×R) -v n,ε L 2 (R 2 ×R) ≤ (1 + κ ε ) 1 2 ε. Now, u n L 2 (R 2 ×R) - R 2 (v 2 2,e (x) + κ ε ) |D (2,0) u| 2 dx 1/2 ≤ u n L 2 (R 2 ×R) -v n,ε L 2 (R 2 ×R) (A.1) + v n,ε L 2 (R 2 ×R) - R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε | 2 dx 1/2 + ( R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε | 2 dx) 1/2 -( R 2 (v 2 2,e (x)+κ ε ) |D (2,0) u| 2 dx) 1/2 . Also, R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε +D (2,0) u| 2 dx= R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε | 2 dx + R 2 (v 2 2,e (x)+κ ε ) |D (2,0) u| 2 dx+2 R 2 (v 2 2,e (x)+κ ε ) D (2,0) u D (2,0) v ε dx, ≤ R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε | 2 dx+ R 2 (v 2 2,e (x)+κ ε ) |D (2,0) u| 2 dx+ +2 ( R 2 (v 2 2,e (x)+κ ε ) |D (2,0) u| 2 dx) 1/2 ( R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε | 2 dx) 1/2 ,
from Cauchy-Schwarz inequality. It follows that

R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε + D (2,0) u| 2 dx 1/2 ≤ R 2 (v 2 2,e (x) + κ ε ) |D (2,0) u| 2 dx 1/2 + R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε | 2 dx 1/2 .
It results from this that:

R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε | 2 dx 1/2 =( R 2 (v 2 2,e (x)+κ ε ) |D (2,0) v ε -D (2,0) u+D (2,0) u| 2 dx) 1/2 , ≤ R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε -D (2,0) u| 2 dx 1/2 + R 2 (v 2 2,e (x) + κ ε ) |D (2,0) u| 2 dx 1/2
.

Exchanging the role of v ε and u, we obtain:

R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε | 2 dx 1/2 - R 2 (v 2 2,e (x) + κ ε ) |D (2,0) u| 2 dx 1/2 ≤ R 2 (v 2 2,e (x) + κ ε ) |D (2,0) v ε -D (2,0) u| 2 dx 1/2 , ≤ (1 + κ ε ) 1 2 ε.
Combining these different elements in (A.1) yields the desired result. Since Ω ∈ C 2 , then ūe is the extension by 0 of ū from W 2,2 (Ω) to W 2,2 (R 2 ). So, D (2,0) ūe = (D (2,0) ū) e and D (2,0) ūe (x) = 0 almost everywhere on R 2 \ Ω. We finally have

R 2 R (v 2,e (x) 2 +κ ε )|ū e,n (x, h)| 2 dh dx -→ n→+∞ R 2 (v 2,e (x) 2 + κ ε )|D (2,0) ūe (x)| 2 dx ≤ Ω (v 2,e (x) 2 + κ ε )|D (2,0) ū(x)| 2 dx. Step 3 It remains to prove that Fε (ū, ḡ, Q, v1 , v2 ) ≤ lim inf n→+∞ Fn,ε (u n , g n , Q n , v 1,n , v 2,n ). From what precedes, it suffices to prove that Ω (v 2 (x) 2 + κ ε )|D (2,0) ū(x)| 2 dx ≤ lim inf n→+∞ R 2 (v 2 2,n,e (x) + κ ε ) R |u n,e (x+he 1 )-2u n,e (x)+u n,e (x-he 1 )| 2 h 4 ρ n (h) dh dx and Ω (v 2 2 (x) + κ ε )|D (0,2) ū(x)| 2 dx ≤ lim inf n→+∞ R 2 (v 2 2,n,e (x) + κ ε ) R |u n,e (x+he 2 )-2u n,e (x)+u n,e (x-he 2 )| 2 h 4 ρ n (h) dh dx. Let η ∈ C ∞ 0 (R 2
) be a non-negative radial function satisfying R 2 η dx = 1, supp(η) ⊂ B(0, 1). We then define a regularized function associated with f δ (x) = 1 δ 2 R 2 f (y)η( x-y δ ) dy, ∀x ∈ R 2 . We focus on the direction e 1 and compute with δ > 0 and δ > 0:

R 2 (v 2 2,n,e,δ (x) + κ ε ) R |u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )| 2 |h| 4 ρ n (h) dh dx = R 2 (v 2 2,n,e,δ (x) + κ ε ) R | B(0,δ)
un,e (x+he 1 -z)-2un,e (x-z)+un,e (x-he 1 -z)

δ 2 η( z δ ) dz| 2 |h| 4
ρ n (h) dh dx.

We apply Jensen's inequality with respect to the measure

η( z δ ) δ 2 dz, yielding ≤ R 2 (v 2 2,n,e,δ (x) + κ ε ) R B(0,δ)
|un,e (x+he 1 -z)-2un,e (x-z)+un,e (x-he 1 -z)| 2

δ 2 |h| 4 η( z δ )ρ n (h) dz dh dx.
We make the following change of variable: Since we assumed that ρ n (h) dh dv ρ n (h) dh dx (2,0) ūe,δ (x)| 2 dx, this last equality being obvious. We have ρ n (h) dh dx.

v = x -z, ≤ R 2 R B(0,δ) (v 2
v 2,n ∈ W 1,∞ (Ω, [0, 1]) with sup n∈N * v 2,n,e W 1,∞ (R 2 ) ≤ sup n∈N * v 2,n W 1,∞ (Ω) ≤ C 1 < ∞ and so sup n∈N * ∇v 2,n,e L ∞ (R 2 ) ≤ C 1 + 1 where v 2,n,e is the exten- sion of v 2,n from W 1,2 (Ω, [0, 1])∩W 1,∞ (Ω) to W 1,2 (R 2 , [0, 1])∩ W 1,∞ (R 2 )
+ 2(C 1 + 1)δC. As W 1,∞ (R 2 ) C(R 2 ), then v 2,
= R 2 (v 2 2,e (x)+κ ε )|D (2,0) ūe,δ (x)| 2 dx = lim n→+∞ R 2 (v 2 2,n,e (x)+ κ ε )|D
R 2 R (v 2
Using Jensen's inequality, we get (2,0) u n,e,δ (x+(t+s-1)he 1 )-D (2,0) ūe,δ (x) ds dt (2,0) u n,e,δ (x+(t+s-1)he 1 )-D (2,0) ūe,δ (x)| 2 ds dt (2,0) u n,e,δ (x + (t + s -1)he 1 ) -D (2,0) ūe,δ (x + (t + s -1)he 1 )| 2 ds dt ρ n (h) dh dx (2,0) ūe,δ (x+(t+s-1)he 1 )-D (2,0) ūe,δ (x)| 2 ds dt ρ n (h) dh dx.

R 2 R 0 1 0 1 0 D
2 ρ n (h) dh dx ≤ R 2 R 0 1 0 1 0 |D
ρ n (h) dh dx, ≤ 2 R 2 R 0 1 0 1 0 |D
+ 2 R 2 R 0 1 0 1 0 |D
We then use Fubini's theorem, Jensen's inequality, and a Taylor's development, and get (2,0) u n,e,δ (x + (t + s -1)he 1 ) -D (2,0) ūe,δ (x + (t + s -1)he 1 )| 2 ds dt ρ n (h) dh dx

≤ 2 R 0 1 0 1 0 R 2 |D
+ 2 R 2 R 0 1 0 1 0 1 0 |D (3,0) ūe,δ (x + k(t + s -1)he 1 )| 2 |h| 2 |t + s -1| 2 dk ds dt ρ n (h) dh dx, ≤ 2 D (2,0) u n,e,δ -D (2,0) ūe,δ 2 L 2 (R 2 ) + 1 6 D (3,0) ūe,δ 2 L 2 (R 2 ) R 0 h 2 ρ n (h) dh, ≤ 2 D (2,0) u n,e,δ -D (2,0) ūe,δ 2 L 2 (R 2 ) + R 2 6 D (3,0) ūe,δ 2 L 2 (R 2 ) .
As u n,e -→ n→+∞ ūe in L 2 (R 2 ) and D (2,0) u n,e,δ = D (2,0) ρ n (h) dh dx

+2 R 2 +∞ R (1 + κ ε )|D (2,0) ūe,δ (x)| 2 ρ n (h) dh dx, ≤ 2 R 2 +∞ R (1 + κ ε ) 1 |h| 4
|u n,e,δ (x + he 1 ) -ūe,δ (x + he 1 ) -2u n,e,δ (x) +2ū e,δ (x) + u -he 1 ) -ūe,δ (x -he 1 )

+ ūe,δ (x+he 1 )-2 ūe,δ (x)+ ūe,δ (x-he 1 )| 2 ρ n (h) dh dx

+ 2(1 + κ ε ) D (2,0) ūe,δ (x) 2 L 2 (R 2 ) +∞ R ρ n (h) dh, ≤ 4 R 2 +∞ R (1 + κ ε ) 1 |h| 4
|u n,e,δ (x + he 1 ) -ūe,δ (x + he 1 )| 2 +4|u n,e,δ (x) -ūe,δ (x)| 2 +4|ū e,δ (x + he 1 ) -2ū e,δ (x) + ūe,δ (x -he 1 )| 2 +4|u n,e,δ (x -he 1 ) -ūe,δ (x -he 1 )| 2 ρ n (h) dh dx ρ n (h) dh dx

+2(1 + κ ε ) D (2,0) ūe,δ (x) 2 L 2 (R 2 ) +∞ R ρ n (h) dh, ≤ 36(1 + κ ε ) R 4 u n,e,δ -ūe,δ 2 L ∞ (R 2 ) + 16(1 + κ ε ) D (2,0) ūe,δ 2 L 2 (R 2 ) + 2(1 + κ ε ) D (2,0) ūe,δ (x) 2 L 2 (R 2 ) +∞ R ρ n (h) dh. Since +∞ R ρ n (h) dh -→ n→+∞ 0, then R 2 +∞ R (v 2
+ 2C(1 + C 1 )δ. Since ūe ∈ W 2,2 (R 2 ) then D (2,0) ūe,δ -→ δ→0 D (2,0) ūe in L 2 (R 2 )
and we deduce by letting δ tend to 0 that ρ n (h) dh dx.

R 2 (v 2 2,e (x) + κ ε )|D (2,0) ūe (x)| 2 dx ≤ lim inf n→+∞ R 2 R (v 2
We finally get lim

n→+∞ Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) = Fε (ū, ḡ, Q, v1 , v2 ) and ∀(u, g, Q, v 1 , v 2 ) ∈ W 2,2 (Ω) ∩W 1,2 0 (Ω) × H(div) × {Q ∈ W 1,∞ (Ω)| Ω Q dx = 0} ×W 1,2 (Ω, [0, 1]) × W 1,2 (Ω, [0, 1]), Fn,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤ Fn,ε (u, g, Q, v 1 , v 2 )
, and by letting n tend to infinity Fε (ū, ḡ, Q, v1 , v2 ) ≤ Fε (u, g, Q, v 1 , v 2 ) thanks to Theorem 1. This concludes the proof.

Appendix B -Proof of Theorem 4: Comparison principle

We follow the arguments of [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF] for parabolic equations and [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]. We first remark that for any λ > 0, q-= q --λ T -t is also a sub-solution of (EE) and satisfies:

q- t + G * (x, t, D q-, D 2 q-) ≤ - λ (T -t) 2 ≤ - λ T 2 , since q- t = q - t - λ (T -t) 2
, D q-= Dq -and D 2 q-= D 2 q -. As q -≤ q + comes from q-≤ q + in the limit when λ tends to 0, it is sufficient to prove the comparison under the additional assumption:

(i) q - t + G * (x, t, Dq -, D 2 q -) ≤ -λ T 2 (ii) lim t→T q -(x, t) = -∞ . Let us set M = sup (x,t)∈R 2 ×[0,T )
q -(x, t) -q + (x, t). We want to prove that M ≤ 0. To do so, we argue by contradiction and assume that M > 0 so there exists (x * , t * ) ∈ R 2 × [0, T ) such that q -(x * , t * ) -q + (x * , t * ) > 0.

We introduce M 0 defined by M 0 = sup

(x,y,t)∈R 2 ×R 2 ×[0,T ) f (x, y, t) = q -(x, t)-q + (y, t)-1 4 |x- y| 4 -α 2 (|x| 2 +|y| 2 ) , > 0.
We notice that M 0 ≥ q -(x * , t * )q + (x * , t * ) -α|x * | 2 > 0 and so M 0 > 0 for α small enough. This supremum is reached owing to the term α 2 (|x| 2 + |y| 2 ), M q -the bound above of q -and -m q + the bound above of -q + , the upper semicontinuity of f , the fact that f is proper and lim

t→T ,|x|→+∞,|y|→+∞ f (x, y, t) = -∞ since f (x, y, t) ≤ M q --m q + -α 2 (|x| 2 + |y| 2
). We denote by (x 0 , y 0 , t 0 ) ∈ R 2 × R 2 × [0, T ) a maximum. Consequently, we get the following lemma:

Lemma 2 (adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF])

Let M = lim h→0 sup |x-y|≤h,t∈[0,T ) (q -(x, t) -q -(y, t)). Then 1. lim α→0 αx 0 = lim α→0 αy 0 = 0. 2. lim →0 |x 0 -y 0 | 4 = 0. 3. lim →0 lim α→0 M 0 = M . 4. lim →0 lim α→0 1 |x 0 -y 0 | 4 = 0. 5. lim →0 lim α→0 α(|x 0 | 2 + |y 0 | 2 ) = 0.
boundedness of q -and q + , the function (x, y) → q -(x, t) -q + (y, t) is bounded for any t ∈ [0, T ). Besides, we assume that M 0 > 0 (α is assumed small enough) . We also

have |x 0 -y 0 | 4 4 + α 2 (|x 0 | 2 + |y 0 | 2 ) ≤ C, C
depending on the bounds of q -and q + . We deduce that α

Let us now set

M h = sup |x-y|≤h,t∈[0,T ) (q -(x, t) -q -(y, t)). Let (x h n , y h n , t h n ) ∈ R 2 × R 2 × [0, T ) be a sequence such that ∀n ∈ N * , q -(x h n , t h n ) -q + (y h n , t h n ) ≥ M h -1 n with |x h n - y h n | ≤ h. This sequence is independent of α. We then get M h -1 n -h 4 4 -α 2 (|x h n | 2 + |y h n | 2 ) ≤ q -(x h n , t h n ) -v(y h n , t h n ) - |x h n -y h n | 4 4 -α 2 (|x h n | 2 + |y h n | 2 ) ≤ M 0 ≤ q -(x 0 , t 0 ) -q + (y 0 , t 0 ). Let α tend to 0 in what precedes, then M h -1 n -h 4 4 ≤ lim inf α→0 (q -(x 0 , t 0 )-v(y 0 , t 0 )) ≤ lim sup α→0 (q -(x 0 , t 0 )-v(y 0 , t 0 )).
Now let h tend to 0, we get M -1 n ≤ lim inf α→0 (q -(x 0 , t 0 )q + (y 0 , t 0 )) ≤ lim sup α→0 (q -(x 0 , t 0 ) -q + (y 0 , t 0 )). Finally, we let tend to 0 and obtain:

M - 1 n ≤ lim inf →0 lim inf α→0 (q -(x 0 , t 0 ) -q + (y 0 , t 0 )) ≤ lim inf →0 lim sup α→0 (q -(x 0 , t 0 ) -q + (y 0 , t 0 )) ≤ lim sup →0 lim sup α→0 (q -(x 0 , t 0 ) -q + (y 0 , t 0 )) ≤ lim sup →0 lim sup α→0 sup |x-y|≤C 1 4 ,t∈[0,T ) (q -(x, t) -q + (y, t)) ≤ lim sup h→0 sup |x-y|≤h,t∈[0,T ) (q -(x, t) -q + (y, t)) = M .
So, lim →0 lim α→0 (q -(x 0 , t 0 )-q + (y 0 , t 0 )) = M . In the same way, we get lim

→0 lim α→0 M 0 = M . Thus lim →0 lim α→0 |x 0 -y 0 | 4 4 + α 2 (|x 0 | 2 + |y 0 | 2 ) = lim →0 lim α→0 (q -(x 0 , t 0 ) -q + (y 0 , t 0 ) -M 0 ) = 0. So    lim →0 lim α→0 |x 0 -y 0 | 4 4 = 0 lim →0 lim α→0 α 2 (|x 0 | 2 + |y 0 | 2 ) = 0 .
We then distinguish two cases.

1. ∀ > 0, ∃α ∈ (0, ) such that t 0 = 0. Then there exists n -→ n→+∞ 0 and α n -→ n→+∞ 0 with t n = 0 and 0

< M 0 ≤ q -(x 0,n , 0)-q + (y 0,n , 0) ≤ Q 0 (x 0,n )-Q 0 (y 0,n ) ≤ B 0 |x 0,n -y 0,n | where B 0 is the Lipschitz constant of Q 0 .
We obtain a contradiction since lim n→+∞ |x 0,n -y 0,n | = 0.

2. ∃ > 0 such that ∀α ∈ (0, ), t 0 > 0. We can choose small enough (otherwise we use the first argument), i.e., such that |x 0 -y 0 | 4 ≤ λ 2T 2 γκ g . We consider q-(x, t) = q -(x, t) -α 2 |x| 2 and q+ (x, t) = q + (x, t) + α 2 |x| 2 and so M 0 = sup

(x,y,t)∈R 2 ×R 2 ×[0,T ) q-(x, t) -q+ (y, t) -|x-y| 4 4 .
Let us take the test function ψ(x, y, t) = |x-y| 4 4 and set p 0 = x 0 -y 0 . By using the parabolic version of Ishii's lemma ([15, Theorem 8.3.]) with the same notations, it comes that τ = 0,

p 1 = |p 0 | 2 p 0 = p 2 , A = 2 |p 0 | 2 Z -Z -Z Z with Z = I 2 + p 0 ⊗p 0 |p 0 | 2
and for each β > 0 such that βA < I, there exist (X, Y ) ∈ (S 2 ) 2 and two reals τ 1 and τ 2 such that

τ 1 -τ 2 = 0, (τ 1 , |p 0 | 2 p 0 , X) ∈ P+ ũ(x 0 , t 0 ), (τ 2 , |p 0 | 2 p 0 , Y ) ∈ P-ṽ(y 0 , t 0 ), -1 β I 0 0 I ≤ X 0 0 -Y ≤ (I -βA) -1 A.
The last inequality implies that X Y . We then use the following lemma. Lemma 3 (adapted from [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]) We have the following estimations on the matrix A:

1 A A I, A ≤ 6|p 0 | 2 , If δ = 1 2 A , then (I -δA) -1 A 2 A I ≤ 12
ξ = ξ 1 ξ 2 ∈ R 4 , ξ i ∈ R 2 , i = {1, 2}.
We have:

Aξ • ξ = 2 |p 0 | 2 Z(ξ 1 -ξ 2 ) • (ξ 1 -ξ 2 ), ≤ 4 |p 0 | 2 Z |ξ| 2 .
So it suffices to show that Z ≤ 3 2 .

Zξ 1 • ξ 1 = |ξ 1 | 2 2 + p 0 ⊗ p 0 |p 0 | 2 ξ 1 • ξ 1 ≤ |ξ 1 | 2 2 + |ξ 1 | 2 , ≤ 3 2 |ξ 1 | 2 .
For the last point, it suffices to notice that if B ≥ 0 and C ≥ 0 with B, C ∈ S 2 such that BC = CB, then CB ≥ 0. Indeed, CBξ•ξ = CB By using this lemma with β = 1 2 A , the inequality becomes

-12|p 0 | 2 I 0 0 I X 0 0 -Y 12|p 0 | 2 I 0 0 I . Thus X and Y are bounded independently of α (since p 0 is boun- ded independently of α having |x 0 -y 0 | 4 ≤ λ 2T 2 γκ g ). In particular, we have -12|p 0 | 2 I X 12|p 0 | 2 I. So there ex- ists a sequence (α n ) such that α n -→ n→+∞ 0, t 0 -→ n→+∞ t ∞ , p 0 -→ n→+∞ p ∞ , X -→ n→+∞ X ∞ and Y -→ n→+∞ Y ∞ .
Furthermore, since q -a sub-solution, q + is a solution and using the additional assumption, we get:

τ 1 + G * (x 0 , t 0 , |p 0 | 2 p 0 + α n x 0 , X + α n I) ≤ -λ T 2 , τ 2 + G * (y 0 , t 0 , |p 0 | 2 p 0 -α n y 0 , Y -α n I) ≥ 0.
Then using the matrix inequality X Y and the ellipticity of the equation, we obtain τ 2 + G * (y 0 , t 0 , |p 0 | 2 p 0α n y 0 , X -α n I) ≥ 0. Subtracting these two inequalities yields successively

λ T 2 + G * x 0 , t 0 , |p 0 | 2 p 0 + α n x 0 , X + α n I ≤ G * y 0 , t 0 , |p 0 | 2 p 0 -α n y 0 , X -α n I , λ T 2 ≤ -γdivg(x 0 , t 0 ) + γdivg(y 0 , t 0 ) -E(X + α n I) -F * |p 0 | 2 p 0 + α n x 0 , X + α n I + E(X -α n I) + F * |p 0 | 2 p 0 -α n y 0 , X -α n I .
Let n tend to infinity. As divg is Lipschitz continuous uniformly in time with κ g the Lipschitz constant, we get:

γκ g |p ∞ | -F * |p ∞ | 2 p ∞ , X ∞ + F * |p ∞ | 2 p ∞ , X ∞ ≥ λ T 2 .
We distinguish here two cases:

-First case: p ∞ = 0. From the matrix inequality

-12|p 0 | 2 I 0 0 I X 0 0 -Y 12|p 0 | 2 I 0 0 I , it yields X ∞ = Y ∞ = 0.
Having F * (0, 0) = F * (0, 0) = 0 leads to 0 ≥ λ T 2 , which is absurd. . Letting tend to 0 raises again a contradiction.

-Second case: p ∞ = 0. Since F is continuous on R 2 \ 0 R 2 ×S 2 , then F * |p ∞ | 2 p ∞ , X ∞ = F * |p ∞ | 2 p ∞ , X ∞ = F |p ∞ | 2 p ∞ , X ∞ . Now, |p 0 | ≤ λ 2T 2

Appendix C -Proof of Theorem 7: Lipschitz regularity in space

We follow the arguments of [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]Lemma 4.15]. We have proved that Q is bounded and continuous on R 2 × [0, T ). Let us set Φ (x, y, t) = B(t)(|x -y| 2 + 2 ) 1 2 . We aim to show that Q(x, t)-Q(y, t) ≤ Φ (x, y, t) for any (x, y, t) ∈ R 2 ×R 2 ×[0, T ).

We introduce M = sup (x,y)∈R 2 ×R 2 ,t∈[0,T ) (Q(x, t) -Q(y, t) -Φ (x, y, t)).

We assume that M > 0. We denote by M = sup We now prove that t > 0 by contradiction. Let us assume that t = 0. Then Q(x, 0) -Q(ȳ, 0) -Φ (x, ȳ, 0) > 0, that is to say, Q 0 (x) -Q 0 (ȳ) > B(0)(|x -ȳ| 2 + 2 ) 1 2 > B 0 |x -ȳ|, which is absurd since ∇Q 0 L ∞ (R 2 ) ≤ B 0 . We set p = D x Φ (x, ȳ, t) = -D y Φ (x, ȳ, t) = B( t)(xȳ)(|x-ȳ| 2 + 2 ) -1 2 = 0, Z = D 2 x Φ (x, ȳ, t) = D 2 y Φ (x, ȳ, t) = B( t) (|x-ȳ| 2 + 2 ) - As the matrix inequality implies X Y and using the ellipticity of the equation, we deduce that τ 2 + divg(ȳ, t) + E(X -αI) + F * (p -αȳ, X -αI) ≥ 0.

Subtracting both inequalities and using again the Lipschitz continuity of divg, we get successively λ (T -t) 2 + γκ g (|x -ȳ| 2 + 2 ) 1 2 -γκ g |x -ȳ| + E(X + αI) -E(X -αI) + F * (p + αx, X + αI) -F * (p -αȳ, X -αI) ≤ 0, λ T 2 + E(X + αI) -E(X -αI) +F * (p + αx, X + αI) -F * (p -αȳ, X -αI) ≤ 0.

Letting α tend to 0 (p and X are bounded independently of α so we can extract convergent subsequences) and using similar arguments to those in [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF] to handle F * (p + αx, X + αI) -F * (p -αȳ, X -αI) yields λ T 2 ≤ 0, which is absurd. Then Q(x, t) -Q(y, t) ≤ Φ (x, y, t). Letting tend to 0, we get Q(x, t) -Q(y, t) ≤ B(t)|x -y|. By exchanging x and y in what precedes, we get |Q(x, t) -Q(y, t)| ≤ B(t)|x -y|, which concludes the proof.

Appendix D -Proof of Theorem 8 : Uniform continuity in time We follow again the arguments of [START_REF] Forcadel | Dislocation Dynamics with a Mean Curvature Term: Short Time Existence and Uniqueness[END_REF]Lemma 4.15]. Let us set δ > 0. For any (x, t) ∈ R 2 × (0, T ) such that t + δ ≤ T , we set q1(x, t) = Q(x, t + δ). Then q1 is a sub-solution of w tω div g (δ) + divg(x, t) + E(D 2 w) + F (Dw, D 2 w) = 0 on R 2 × (0, T -δ) with w(x, 0) = Q(x, δ) on R 2 in the sense of viscosity solutions theory. Indeed, we have q 1,t (x, t) = Q t (x, + δ), q 1,t + divg(x, t + + E(D 2 q 1 ) + F (Dq 1 , D 2 q 1 ) = 0 and divg(x, t + δ) ≥ -ω divg (δ) + divg(x, t) leading to q 1,tω divg (δ)+divg(x, t)+E(D 2 q 1 )+F (Dq 1 , D 2 q 1 ) ≤ 0. Besides, q = Q + sup x∈R 2 ((Q(x, δ) -Q 0 (x)) + ) + ω divg (δ)t (-(α) + = max{0, α} -) is a super-solution with q1(x, 0) ≤ q(x, 0) for any x ∈ R 2 . Using the comparison principle, we get, ∀(x, t) ∈ R 2 ×(0, T -δ), Q(x, t+δ)-Q(x, t) ≤ sup x∈R 2 ((Q(x, δ)-Q 0 (x)) + )+ ω divg (δ)t. We will now focus on an auxiliary problem. Let us first assume that Q 0 ∈ C 2 b (R 2 ). We set q +/-= Q 0 (x) + / -

C 1 t with C 1 = inf x∈R 2 {-E(D 2 Q 0 ) -F * (DQ 0 , D 2 Q 0 ), E(D 2 Q 0 ) + F * (DQ 0 , D 2 Q 0 )}.
We easily check that q + is a super-solution and q -a sub-solution of (SP ): u t + E(D 2 u) + F (Du, D 2 u) = 0, u(x, 0) = Q 0 (x). Then there exists a unique solution q of this problem and by the comparison principle, the following holds: ∀t ∈ [0, T ), ∀x ∈ R 2 , |q(x, t) -Q 0 (x)| ≤ C 1 t. We then set q1(x, t) = q(x, t + h). So q1 is also a solution to the problem (SP ) with q1(x, 0) = q(x, h), ∀x ∈ R 2 and by the comparison principle, we get q(x, t + h) -q(x, t) ≤ sup x∈R 2 {q(x, h) -Q 0 (x)} ≤ C 1 h. Similarly, we have q(x, t)q(x, t+h) ≤ sup x∈R 2 {q(x, h)-Q 0 (x)} ≤ C 1 h and so |q(x, t+h)q(x, t)| ≤ sup x∈R 2 {q(x, h) -Q 0 (x)} ≤ C 1 h. We now assume that Q 0 ∈ W 1,∞ (R 2 ). We set Q 0 = Q 0 * ρ , where ρ is a regularizing sequence satisfying ρ = 1 2 ρ( • ), ρ ∈ C ∞ c (R 2 , R), ρ ≥ 0, supp(ρ) ⊂ B(0, 1), and R 2 ρ(x) dx = 1. Then We denote by q the solution with initial condition Q 0 . Then, by the comparison principle, q (., t)-q (., t)

Q 0 ∈ C 2 c (R 2 ) and DQ 0 L ∞ (R 2 ) ≤ B 0 , D 2 Q 0 L ∞ (R 2 ) ≤ B 0 C 2 .
L ∞ (R 2 ) ≤ Q 0 - Q 0 L ∞ (R 2 )
, and so (u ) uniformly converges to q the solution of (SP ) with initial condition Q 0 since (Q 0 ) converges uniformly to Q 0 . We then have by the comparison principle that q (., t) -q(., t)

L ∞ (R 2 ) ≤ Q 0 -Q 0 L ∞ (R 2
) and deduce that q(., t + h) -q(., t)

L ∞ (R 2 ) ≤2 Q 0 -Q 0 L ∞ (R 2 )
+ q (., t + h) -q(., t) L ∞ (R 2 ) ,

≤2B 0 + C 1 B 0 , B 0 C 2 h.
By taking the minimum on , we obtain the modulus of continuity of q, called ω F which depends only on B 0 . And so q satisfies Q 0 (x) -ω F (t) ≤ q(x, t) ≤ Q 0 (x) + ω F (t). We now prove that Dq L ∞ (R 2 ×(0,T )) ≤ DQ 0 L ∞ (R 2 ) . Indeed, let us consider the function q h (x, t) = q(x+h, t)+ DQ 0 L ∞ (R 2 ) |h|. Then q h is still a solution of (SP ) with initial condition q h (x, 0) = Q 0 (x + h) + DQ 0 L ∞ (R 2 ) |h| ≥ Q 0 (x). So, by the comparison principle, we have q h ≥ q. We deduce that for every h ∈ R 2 , q(x, t) -q(x + h, t) ≤ DQ 0 L ∞ (R 2 ) |h|, and so Dq L ∞ (R 2 ×(0,T )) ≤ DQ 0 L ∞ (R 2 ) . We then set Q + (x, t) = q(x, t) + divg L ∞ (R 2 ×[0,T )) t and come back to the initial problem. Then DQ + L ∞ (R 2 ×(0,T )) ≤ Dq L ∞ (R 2 ×(0,T )) ≤ B 0 and Q + is a super-solution of (EE) since Q + (x, 0) = q(x, 0) = Q 0 (x) and Q + t (x, t) + G * (x, t, Q + , DQ + , D 2 Q + ) = q t (x, t) + divg L ∞ (R 2 ×[0,T )) + divg(x, t) + E(D 2 q(x, t)) +F * (Dq(x, t), D 2 q(x, t)) = divg L ∞ (R 2 ×[0,T )) +divg(x, t) ≥ 0 and satisfies Q + (x, t) ≤ Q 0 (x)+ω F (t)+ divg L ∞ (R 2 ×[0,T )) t. Similarly, we construct a sub-solution Q -such that Q -(x, t) ≥ Q 0 (x) -ω F (t) -divg L ∞ (R 2 ×[0,T )) t by setting Q -(x, t) = q(x, t) -divg L ∞ (R 2 ×[0,T )) t. By applying the comparison principle, we get that

Q 0 (x) -ω F (t) -divg L ∞ (R 2 ×[0,T )) t ≤ Q -(x, t) ≤ Q(x, t) ≤ Q + (x, t) ≤ Q 0 (x) + ω F (t) + divg L ∞ (R 2 ×[0,T )) t.

Thus sup

x∈R 2 ((Q(x, δ)-Q 0 (x)) + ) ≤ ω F (δ)+ divg L ∞ (R 2 ×[0,T )) δ
and so Q(x, t + δ) -Q(x, t) ≤ ω F (δ) + divg L ∞ (R 2 ×[0,T )) δ + ω divg (δ)T . Similarly, v is a super-solution of w t + ω divg (δ) + divg(x, t)+E(D 2 w)+F (Dw, D 2 w) = 0 with w(x, 0) = Q(x, δ) on R 2 and ũ = Q(x, t)-sup

Algorithm 1

 1 Computation of the nonlocal weights inspired by the NL-means algorithm. Input: Initial image f Output: Weights in the first direction w f ,x,1 , the shifts of the selected neighbors in the first direction indice 1 , weights in the second direction w f ,x,2 , the shifts of the selected neighbors in the second direction indice 2 . 1. Define w := window size, p := patch size, h := 0.25, N bN eigh := number of actual required neighbors including the closest one. 2. Compute the extended image by symmetry.

Fig. 1 :

 1 Fig. 1: Cartesian topology in the local case with 4 processes on a (10 × 10) image before reversing the row axis: Process 0: (sx = 2, ex = 4, sy = 5, ey = 7); Process 1: (sx = 5, ex = 7, sy = 5, ey = 7); Process 2: (sx = 2, ex = 4, sy = 2, ey = 4); Process 3: (sx = 5, ex = 7, sy = 5, ey = 7).

Fig. 2 :

 2 Fig. 2: Phantom cells and communications, example on process 0 for the local case.

  Original image (b) u local (c) v = div g local (d) v 2 local (e) u nonlocal (weights based on gradient) (f) v 2 nonlocal (weights based on gradient) (g) u nonlocal (weights based on Hessian) (h) v 2 nonlocal (weights based on Hessian) (i) u nonlocal (only 4 closest neighbors) (j) v 2 nonlocal (only 4 closest neighbors)

( k )Fig. 3 :

 k3 Fig. 3: Crack detection: µ = 0.001, ξ ε = 3.5, α = 0.14, β = 0.07, ρ = 3.5, ε = 1, γ = 0.5, 50 iterations; nonlocal (weights based on Hessian): µ = 0.0001, ξ ε = 3.5, α = 0.14, β = 0.07, ρ = 3.5, ε = 1.5, γ = 1.0, w = 15, p = 7, 550 iterations; nonlocal (only 4 closest neigbors): µ = 0.0001, ξ ε = 3.5, α = 0.06, β = 0.03, ρ = 3.5, ε = 1.5, γ = 1.0, 550 iterations; nonlocal (weights based on gradient): µ = 0.0001, ξ ε = 3.5, α = 0.06, β = 0.03, ρ = 3.5, ε = 1.0, γ = 1.0, w = 9, p = 5, 1000 iterations.

  g)-(h)-(i) and 6-(g)-(h)-(i) are the results obtained without considering the decomposition part. Thanks to Figure 5-(j)-(k)-(l) and 6-(j)-(k)-(l)showing the absolute difference between both results, we observe that u is less noisy with our method, v 1 and v 2 also exhibit better contrast with less superfluous information. The results obtained with the three versions of our nonlocal algorithm are comparable with the ones obtained with our local implementation.

  Original image (b) u local (c) v = div g local (d) v 2 local (e) u nonlocal (weights based on gradient) This work was devoted to crack detection on bituminous surfacing images. To overcome the limitation of classical first order methods, namely the inability to capture singularities without jumps, and motivited by Drogoul et al.'s observation, a second order appproach has been proposed based on a modified elliptic approximation of the Blake-Zisserman functional with a variable encoding the fine structures. A decomposition model has been added to handle the dense and highly

( k )Fig. 4 :Fig. 5 :Fig. 6 :

 k456 Fig. 4: Crack detection: µ = 0.001, ξ ε = 3.5, α = 0.14, β = 0.07, ρ = 3.5, ε = 1, γ = 0.5, 50 iterations; nonlocal (weights based on Hessian): µ = 0.0001, ξ ε = 1.0, α = 0.06, β = 0.03, ρ = 1.0, ε = 1.5, γ = 1.0, w = 15, p = 7, 750 iterations; nonlocal (only 4 closest neighbors): µ = 0.0001, ξ ε = 1.0, α = 0.06, β = 0.03, ρ = 1.0, ε = 1.5, γ = 1.0, 330 iterations; nonlocal (weights based on gradient): µ = 0.0001, ξ ε = 3.5, α = 0.06, β = 0.03, ρ = 3.5, ε = 1.0, γ = 1.0, w = 9, p = 5, 1000 iterations.

  , with λ < 1 2 , then (u n ) strongly converges to ū in C 0,λ b

2 1δ 2 D 2 1δ 2 D 2 ( 2 R (v 2 2 ,

 2222222 and v2,e is the extension of v2 from W1,2 (Ω, [0, 1]) to W 1,2 (R 2 , [0, 1]), then we have thanks to[22, Theorem 1 p.123] that ∀x ∈ R 2 , D (1,0) v 2,n,e,δ (x) = η δ * D (1,0) v 2,n,e (x) = R (1,0) v 2,n,e (y)η( x-y δ ) dy ≤ D (1,0) v 2,n,e L ∞ (R 2) , D (0,1) v 2,n,e,δ (x) = η δ * D (0,1) v 2,n,e (x) = R (0,1) v 2,n,e (y)η( x-y δ ) dy ≤ D (0,1) v 2,n,e L ∞ (R 2 ) and so ∇v 2,n,e,δ L ∞ (R 2 ) ≤ ∇v 2,n,e L ∞ (R 2 ) ≤ C 1 + 1. Eventually, we get:R v 2,n,e,δ (x) 2 + κ ε ) R |u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )| 2 |h| 4 ρ n (h) dh dx ≤ R n,e,δ (v) + κ ε ) |u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4

2 R 2 +∞R

 22 2,n,e (x) + κ ε )|u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )-h 2 D (2,0) ūe,δ (x)u n,e,δ (x+(t+s-1)he 1 )-D(2,0) ūe,δ (x) ds dt 2 ρ n (h) dh dx + (v 2 2,n,e (x) + κ ε ) |u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )-h 2 D (2,0) ūe,δ (x)| 2 |h| 4

2 RD 2 |h| 4 ρ n (h) dh dx ≤ 2 R 2 +∞R( 1 +

 224221 η δ * u n,e , then D (2,0) n,e,δ -D (2,0) ūe,δ L 2 (R 2 ) -→ n→+∞ 0. Also by letting R tend to 0, lim n→+∞ R (2,0) u n,e,δ (x + (t + s -1)he 1 ) -D (2,0) ūe,δ (x) ds dt 2 ρ n (h) dh dx = 0. Fur-thermore, ,e (x) + κ ε ) |u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )-h 2 D (2,0) ūe,δ (x)| κ ε ) |u n,e,δ (x+he 1 )-2u n,e,δ (x)+u n,e,δ (x-he 1 )| 2 |h| 4

2 (

 2 |x 0 | 2 + |y 0 | 2 ) ≤ C and |x 0 -y 0 | 4 ≤ 4C leading to lim →0 |x 0 -y 0 | = 0 and lim α→0 αx 0 = lim α→0 αy 0 = 0.

|p 0 | 2

 2 ξ•ξ , the • denoting the inner product in R 4 . Proof By definition of A , we have Aξ•ξ ξ•ξ ≤ A for any ξ ∈ R 4 \ {0 R 4 } and so Aξ•ξ A Iξ • ξ, which gives the first result of the lemma. Let

1 2 ξ•B 1 2

 11 ξ ≥ 0 since C is nonnegative and B is symmetric non-negative. So, if B C and D ≥ 0 with D(B -C) = (B -C)D then DB DC. It thus suffices to prove that A 2 A (I -δA) = 2 A I -A, i.e., A A I which is true by definition of the norm.

(x,y)∈R 2

 2 ×R 2 ,t∈[0,T ) {Q(x, t)-Q(y, t)-α 2 (|x| 2 +|y| 2 )-λ T -t -Φ (x, y, t)}.For λ and α small enough, we have M > 0. As Q is bounded and continuous with lim|x|→+∞,|y|→+∞,t→T Q(x, t)-Q(y, t) -α 2 (|x| 2 + |y| 2 ) -λ T -t -Φ (x, y, t) = -∞,the maximum is attained at (x, ȳ, t) with x = ȳ and α 2 (|x| 2 +|ȳ| 2 ) ≤ C owing to M > 0 and Q bounded. Therefore, lim

  )

2 + 2 ,

 22 γκ g (|x -ȳ| 2 + 2 ) 1 (τ 1 , p + αx, X + αI) ∈ P+ Q(x, t), (τ 2 , p -αȳ, Y -αI) ∈ P-Q(ȳ, t), -βA) -1 A.Thus we haveτ 1 + divg(x, t) + E(X + αI) + F * (p + αx, X + αI) ≤ 0, τ 2 + divg(ȳ, t) + E(Y -αI) + F * (p -αȳ, Y -αI) ≥ 0.

2 DQ 0 2 |DQ 0 2 ρ 2 DQ 0 2 |DQ 0 B 0 2 |Q 0 ρ

 202022020020 Indeed using properties stated in[START_REF] Evans | Measure Theory and Fine Properties of Functions[END_REF] Theorem 1, p. 123], we have|DQ 0 (x)| = |DQ 0 * ρ (x)| = | R (x -y)ρ (y) dy|, ≤ R (x -y)ρ (y)| dy ≤ B 0 R (y) dy, ≤ B 0 ,and|D 2 Q 0 (x)| = |DQ 0 * Dρ (x)| = | R (x -y)Dρ (y) dy|, ≤ R (x -y)Dρ (y)| dy, Dρ L 1 (R 2 ) . Moreover, Q 0 -Q 0 L ∞ (R 2 ) ≤ B 0 since |Q 0 (x) -Q 0 (x)| ≤ R (x) -Q 0 (x -y)|ρ (y) dy, (y) dy = B 0 .

x∈R 2 ( 2 (

 22 (Q(x, δ)-Q 0 (x)) -)-ω divg (δ)t ((α) -= min{0, α}) is a sub-solution. So, by the comparison principle, we have Q(x, t) -Q(x, t + δ) ≤ sup x∈R (Q(x, δ) -Q 0 (x)) -) + ω divg (δ)t, ≤ω F (δ) + divg L ∞ (R 2 ×[0,T )) δ + ω divg (δ)T.And so,|Q(x, t) -Q(x, t + δ)| ≤ω F (δ) + divg L ∞ (R 2 ×[0,T )) δ + ω divg (δ)T,which achieves the proof.

  . Once the Cartesian topology is created, one needs to distribute the data file (image data) to each subdomain in parallel (more precisely, the portion of the data file that must be visible for the related process). The general file manipulation function MPI File open is called. We then create a datatype MPI Datatype mysubarray describing a two-dimensional subarray (the portion of the image related to the current subdomain) of a bigger twodimensional array (the image here) (routine MPI Type create subarray). The MPI File set view routine allows to change the process view of the data in the file: the beginning of the data accessible in the file through that view is set to 0, the type of data is set to MPI DOUBLE, and the distribution of data to processes is set to mysubarray. Then the MPI File read routine enables us to read the file starting at the specified location. (v) Derived datatypes are created, describing the rows, ey) is sent to the eastern neighbor that receives data in u local mat(:, sy -w + 1 : sy -1); (e) u local mat(sx : sx + 1, ey -1 : ey) is sent to the northeast neighbor that receives data in u local mat(ex + 1 : ex + 2, sy -2 : sy -1); (f) u local mat(ex -1 : ex, ey -1 : ey) is sent to the southeast neighbor that receives data in u local mat(sx -2 : sx -1, sy -2 : sy -1); (g) u local mat(ex -1 : ex, sy : sy + 1) is sent to the southwest neighbor that receives data in u local mat(sx -2 : sx -1, ey + 1 : ey + 2);

		For subdomains with at least one edge included in
		the image domain boundary, diagonal communica-
		tions do not occur and are replaced by communica-
		tions with the ghost cells of the involved contiguous
	MPI communications: MPI Datatype type (w-1)colonnes, type ligne, type colonne mono, type (w-1)lignes among others. In that purpose, the routines MPI Type contiguous -creating a contiguous datatype, here a single row of data -and MPI Type vector -general constructor that al-lows replication of a datatype into locations that consist of equally spaced blocks, here a single co-lumn, a group of (w-1) adjacent columns and a group of (w-1) adjacent rows -are used. For the commu-nications with diagonally connected subdomains, a datatype describing a two-dimensional subarray of size 2×2 of a bigger two-dimensional array is created for each spatial configuration: top left corner, top right corner, bottom left corner and bottom right corner. (iv) columns and 2×2 diagonal arrays involved in the ey sy	Created function: void voisinage(MPI Comm comm2d, int * voisin, int * voisin diagonale, int * coords, int * dims) subdomain. Created function : void communication(double ** u, double ** v1 local mat, double ** v2 local mat, double ** g1 local mat, dou-ble ** g2 local mat, double ** Q local mat, int * tab bounds, MPI Comm comm2d,int * voisin, int * voisin diagonale) Received from 2 R e c e i v e d f r o m 1 3 3 Send to 2 S e n d to 1 sx ex

(vi)

The communications are then handled with the routine MPI Sendrecv. This send-receive operation combines in one call the sending of a message to one destination and the receiving of another message from another process. As an illustration (see also Figure

2

), u local mat denoting the local array describing u : (a) u local mat(sx : sx + w -2, :) is sent to the northern neighbor that receives data in u local mat(ex + 1 : ex + w -1, :); (b) u local mat(ex -w + 2 : ex, :) is sent to the southern neighbor that receives data in u local mat(sx -w + 1 : sx -1, :); (c) u local mat(:, sy : sy + w -2) is sent to the western neighbor that receives data in u local mat(:, ey + 1 : ey + w -1); (d) u local mat(:, ey -w + 2 : (h) u local mat(sx : sx + 1, sy : sy + 1) is sent to the northwest neighbor that receives data in u local mat(ex + 1 : ex + 2, ey + 1 : ey + 2).

  The newly computed values are then transmitted to neighboring subdomains and a new time step is achieved. (viii) We finally write the result (component v 2 encoding the fine structure we aim to recover) in a file using again MPI I/O. The general file manipulation function MPI File open is called. We then create a new datatype describing the two-dimensional subarray extracted from the original local array v2 local mat when removing the ghost cells. The routine MPI File set view changes the process view of the data in the file, while the routine MPI File write all writes the file, starting at the locations specified by individual file pointers. Created function : void ecrire mpi(double * v2 local vect,int ntx,int nty,int * tab bounds, MPI Comm comm2d) (ix) The computation of the nonlocal weights is done by every process on its associated subdomain. The computations have been made with the supercomputer Myria operated by the CRIANN (Centre Régional Informatique et d'Applications Numériques de Normandie, http://www.criann.fr/). Myria is an ATOS BULL solution with 11144 computing cores, with a power of 403 TFlops Xeon, 170 TFlops GPU and 27 TFlops Xeon Phi KNL. Myria also has a fast storage space of about 2.5 Po. Submission of the work is done through the SLURM software.

	The local algorithm has been applied to a bituminous
	surfacing image of size 2248 × 4000, requiring 500 time
	step iterations both in the case of the sequential algo-
	rithm, execution time: 289.6 seconds, and parallelized
	algorithm with 224 tasks, execution time: 1.7 seconds.
	The efficiency is of 74% from 1 to 224 tasks and of
	100% from 28 to 224 tasks. The nonlocal algorithm has
	been applied to the same bituminous surfacing image
	of size 2248 × 4000, requiring 360 time step iterations
	both in the case of the sequential algorithm, execution
	time: 596.582511 seconds, and parallelized algorithm
	with 224 tasks, execution time: 4.409619 seconds.
	The numerical applications presented below have been
	obtained in this MPI setting, with the same orders of re-
	duction factor when passing from a single to 224 tasks.

  (size 501 × 501) and Figure 4 (size 285 × 429), courtesy of A. Drogoul. We depict the three main components of the decomposition/segmentation process, i.e., u local, v = div g local, v 2 local, for the local model (related to problem (1) and for the sake of comparison with the nonlocal model), and u nonlocal and v 2 nonlocal, for the three versions of our nonlocal algorithm (weights based on gradient stands for the algorithm in which the weights are com-

puted thanks to the distance between image gradients, weights based on Hessian is for the algorithm in which the weights are given by comparing the Hessians of the image, and only 4 closest neighbors refers to our nonlocal implementation in which we consider only the four closest neighbors), as well as the results obtained with Aubert and Drogoul's topological gradient method

([6,

  L ∞ (R 2 ) δC, from the coercivity inequality.

					≤	R 2 R	(v 2 2,n,e,δ (v) + κ ε )
					|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4	ρ n (h) dh dv
					+	R 2 R B(0,1)	( ∇v 2 2,n,e,δ L ∞ (R 2 ) δ|u|)
					|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4	η(u)ρ n (h) du dh dv,
					≤	R 2 R	(v 2 2,n,e,δ (v) + κ ε )
					|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4	ρ n (h) dh dv
					+ ∇v 2 2,n,e,δ L ∞ (R 2 )
					δ R 2 R	|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4	ρ n (h) dh dv,
					≤	R 2 R	(v 2 2,n,e,δ (v) + κ ε )
					|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 |h| 4	ρ n (h) dh dv
					+ ∇v 2 2,n,e,δ
			2,n,e,δ (v + z) + κ ε )
	|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 δ 2 |h| 4	η( z δ )ρ n (h) dz dh dv.
	Since v 2,n,e,δ ∈ C ∞ (R 2 ), then v 2 2,n,e,δ (v+z) = v 2 2,n,e,δ (v)+
	1 0 ∇v 2 2,n,e,δ (v + sz), z ds and by introducing it in the pre-vious inequality, we get:
	≤	R 2 R B(0,δ)	(v 2 2,n,e,δ (v) + κ ε )
	|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 δ 2 |h| 4	η( z δ ) ρ n (h) dz dh dv
	+	R 2 R B(0,δ)	0 1	( ∇v 2 2,n,e,δ (v + sz), z ) ds
	|u n,e (v+he 1 )-2u n,e (v)+u n,e (v-he 1 )| 2 δ 2 |h| 4	η( z δ )ρ n (h) dz dh dv.
	We now integrate with respect to z in the first integral and use
	Cauchy-Schwarz inequality, and the change of variable u = z δ in the second one, after bounding the component above by
	∇v 2 2,n,e,δ L ∞ (R 2 ) :

  2,n,e (x) + κ ε ) |u n,e (x+he 1 )-2u n,e (x)+u n,e (x-he 1 )| 2 |h| 4

  1 2 I -(|x-ȳ| 2 + 2 ) -3 2 (x-ȳ)⊗(x-ȳ) and A = Z -Z -Z Z .Then by the parabolic version of Ishii's lemma ([15, Theorem 8.3.]) applied to ũ(x, t) = Q(x, t) -, for every β such that βA < I, there existτ 1 ∈ R, τ 2 ∈ R and X ∈ S 2 , Y ∈ S 2 such that τ 1 -τ 2 = λ (T -t) 2 + B ( t)(|x -ȳ| 2 + 2

α 2 |x| 2 , ṽ(x, t) = Q(x, t) + α 2 |x| 2 , and Φ(x, y, t) = Φ (x, y, t) + λ T -t

ϕ 2 (x) θ 2(h) , where θ i : R + → R + and lim h→0 θ i (h) = 0.Let Γ be a smooth closed curve or a smooth infinite curve of R 2 delimiting two subdomains R 2- Γ and R 2+ Γ forming a partition of R 2 . Let then take ϕ as the signed distance function to Γ defined by ϕ(x) = dist(x, R 2- Γ )dist(x, R 2+