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

This paper is devoted to the challenging problem of ﬁne structure detection with applications to bituminous surfacing crack recovery. Drogoul (SIAM J Imag Sci 7(4):2700–2731, 2014) shows that such structures can be suitably modeled 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 (cid:2)(cid:2)(cid:2)(cid:2) ∂ 2 u ∂ x 21 (cid:2)(cid:2)(cid:2)(cid:2) 2 + (cid:2)(cid:2)(cid:2)(cid:2) ∂ 2 u ∂ x 22 (cid:2)(cid:2)(cid:2)(cid:2) 2 , 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 ﬁrst-order nonlocal regularization to higher-order derivatives, we derive and analyze a nonlocal version of the model, and provide several theoretical results among which there are 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. [43], 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 [43] in which an exhaustive classification of segmentation methods into three main categories (fully supervised methods, semi-/weakly supervised methods and 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 nonzero 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 [18,19] for further details). These thin/filament-like structures cannot be captured efficiently with differential operators of order 1. In [18,19], 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 . 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 twodimensional case through the following lemma taken from [18,19] 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 . 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+ Γ ). R 2+ Γ (resp. R 2− Γ ) identifies to the subdomain {ϕ > 0} (resp. {ϕ < 0}). Taking the following scalings θ 1 (h) = √ π h and θ 2 (h) = h, one has g h −→ Since filaments/cracks can be modeled 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 quantity ∂ 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 horizontal direction-and that also handles the natural dense and highly oscillatory texture exhibited by the images. More precisely, we propose weighting ∂ 2 u , 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 modeled by {v 2 = 0}. This leads to a variant of the elliptic approximation of the Blake-Zisserman functional ( [9]) given by Ambrosio,Faina and March ([2]).
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 [25]. It relies on the spectral total variation transform that yields a continuous, multiscale, 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 [22], 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 [6,18], 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 ( [8]) propose a model derived from the Mumford-Shah functional and constrained by geometrical priors to perform the segmentation of tubelike 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 [34] in which nonlocal regularizers enforcing straightness are applied, or detection/completion of thin filaments in noisy blurred 2D or 3D images in [7], relying on the Ginzburg-Landau energy. In [4], 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 [42] for automatic detection of vessel-like patterns, though sensitive to the noise type and time-consuming. In [39], stochastic methods in which a thin network is simulated by a point pro-cess penalizing disconnected segments and favoring aligned pieces are elaborated.
We emphasize that this paper is an extension of [16] not reduced to a simple supplement, containing no overlap with [16] and including significant theoretical and computational aspects. It completes [16] since it includes substantial results, such as a Γ -convergence result (main result) relating the local basis functional to 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 [16] with the derivation of an MPI parallelization of the designed algorithm inspired by the nonlocal means ( [12]), this algorithmic part being entirely new and numerically analyzed. (In particular, the choice of the weights is discussed.) The next section (Sect. 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. ( [10]), Dávila ( [15]) or Ponce ( [33]) devoted to the design of nonlocal counterparts of Sobolev and BV seminorms 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. ( [28])) 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 ( [12]), 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 Sect. 4 that intends to qualitatively assess the proposed algorithm on several experiments in comparison with prior works.

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 For a positive real number 0 < λ < 1, the subspace C 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 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) 2 . Finally, we need the operator called normalized infinity Laplacian and defined as Δ ∞ (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 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 modeling 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 [23,24,27,[29][30][31][32][35][36][37] or [41]). In this work, to model the texture, we propose to use the space G(R 2 ) introduced by Meyer ( [30]), 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 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 [40]: 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 modeling 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. [32] and introduce an auxiliary variable naturally stemming from Helmholtz-Hodge decomposition as follows: g = ∇Q + P, with P a divergence-free vector disregarded afterward. 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: dx, κ ε , ξ ε , ζ suitable 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 modeled 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 minimizers are searched in the respective following functional spaces: (Ω) enabling u to be extended by 0 outside the domain Ω without being too restrictive.
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. The existence of minimizers is guaranteed as stated in [16,Theorem 11].

Toward a Nonlocal Related Model
Inspired by the strength and robustness of nonlocal methods exemplified in many image processing tasks in [26] (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 [10] introducing a new nonlocal characterization of the Sobolev spaces W 1, p , 1 < p < +∞, Aubert and Kornprobst [5] questioning the utility of this characterization to actually solve variational problems, Dàvila [15] and Ponce [33] 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 [28]. This new operator is derivativefree 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, then u ∈ W 2,2 (R 2 ) . We thus introduce a sequence of functionals F ε,n depending on n and such that the component is approximated by an integral operator involving a differential quotient and the radial mollifiers. The resulting model inherits fine analytical properties and has the advantage of being numerically more tractable than [28] 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 [16] 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 restricted to a small area around the point of interest (like the NL-means algorithm [12]) 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 ( [16, Theorem 10]). Theorem 1 Let u ∈ W 2,2 (R 2 ). Then Equipped with this theoretical result and characterization, we are henceforth able to derive a nonlocal counterpart of our local problem (1), phrased in terms of the functional F ε,n defined below and for which an existence result is provided (see [16,Theorem 12]).

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 +∞.
[ be a sequence of minimizers of (2) for each n ∈ N * . Let us assume additionally for technical purposes that v 2,n ∈ {v 2 ∈ Then there exist a subsequence still denoted by (u n , g n , Q n , v 1,n , v 2,n ) and a minimizer (ū,ḡ, Proof The proof is subdivided into three main steps that we sketch hereafter. The details are given in Appendix A.
Step 1. We have proved in what precedes that for any n ∈ N * , there exists a solution to (2). We consider a sequence of such minimizers (u n , g n , Step 2. In a second stage, we demonstrate that Step 3. The last step consists in showing that being a minimizer of (1).

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 x 1 , n x 2 ), the unit outward normal to the boundary. Making use of the absolutely minimizing Lipschitz extensions ( [3]) for the equation in Q, we get: combined with suitable boundary conditions. Focusing now on the variable u, denoting by η a test function, by a real number and setting we obtain, with suitable boundary conditions: 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 [14] for a general introduction including the main material and definitions.

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 with Q 0 ∈ W 1,∞ (R 2 ) and B 0 its Lipschitz constant. (To remove the problem of boundary conditions, we work on R 2 for the spatial domain.) We now give a thorough analysis of the PDE in Q in the viscosity solution theory framework. Let us first state the additional following assumption: div g is bounded and is Lipschitz continuous uniformly in time with κ g its Lipschitz constant independent of time. (H) For the sake of conciseness and using the normalized version of the infinity Laplacian, the evolution equation is now written in the form being the set of symmetric 2 × 2 matrices equipped with its Loewner partial order) defined by ⊗ representing the Kronecker product, and with the following properties : 1. The operators G, E : X → −γ trace(X ) and 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 We begin our analysis by establishing a comparison principle stating that if a subsolution and a super-solution are ordered at the initial time, then they are ordered at any time.

Theorem 4 (Comparison principle, adapted from [21]) Assume (H) and let q
Proof We refer the reader to Appendix B.
We now construct barriers in order to use Perron's method.

respectively, super-and subsolution 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.
Hence, q − is a subsolution and q + a super-solution of (EE).
We are now able to use Perron's method and to give the following existence result. Then function Q = sup{w, w subsolution of (EE) and 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 subsolution. Applying the comparison principle yields Q * ≤ Q * . But by definition, Q * ≤ Q * , so we finally get We now analyze the regularity of the solution. Let us first consider the regularity in space.

Then the solution of (EE) is Lipschitz continuous in space and
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 [21]) 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 well adapted for imaging problems. We recall that they satisfy ∀n ∈ N, ∀x ∈ R, ρ n (x) = ρ n (|x|); ∀n ∈ N, 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, 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 [10,15] and [33]). We thus redefine the weights ρ n (h) |h| 4 involved in the theoretical expression, inspired by the NL-means algorithm ( [12]). 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 behavior as seen in Lemma 1 and that are geographically close. In that purpose, we consider the following nonlocal weights ρ n (h) , 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 time-dependent scheme for Q = Q(x 1 , x 2 , t) (nonlinear over-relaxation method, see [13,Sect. 4]), and a stationary semi-implicit fixed-point scheme in . 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 Algorithm 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.
Compute the distances (we do not make use in practice of the Gaussian kernel G a , yielding equal weights for the contributions) 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 between all patches centered at x + ye 2 of size p with 0 ≤ y ≤ 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 minimizers 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 nt x interior points in the row direction (we removed the first and last rows of points) and nt y 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 Fig. 1 for an example), each subdomain comprising w −1 (w is the window size) rows of the above ghost cells, 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 [1]=ex, tab_bounds [2]=sy and tab_bounds [3]=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 . Created function: void voisinage(MPI_Comm comm2d, int * voisin, int * voisin_diagonale, int * coords, int * dims) (iv) 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 two-dimensional 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, columns and 2 × 2 diagonal arrays involved in the 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 allows replication of a datatype into locations that consist of equally spaced blocks, here a single column, a group of (w-1) adjacent columns and a group of (w-1) adjacent rows-are used. For the communications 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. (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 Fig. 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 : 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); (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).
For subdomains with at least one edge included in the image domain boundary, diagonal communications do not occur and are replaced by communications with the ghost cells of the involved contiguous subdomain. 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.  The numerical applications presented below have been obtained in this MPI setting, with the same orders of reduction factor when passing from a single to 224 tasks. 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 accuracy, 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, in both Fig. 3 (size 501 × 501) and Fig. 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 computed 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 refer 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,18]). 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 (Fig. 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 Fig. 4 is properly returned. We conclude the paper with two applications dedicated to crack detection on bituminous surfacing Figs 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 preprocessing step is applied. It consists in apprehending the problem first as an inpainting one ( [1]) and by considering these white spots as missing parts of the image that need to be filled. This is achieved with the MATLAB ® 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 5g, h, i and 6g, h, i are the results obtained without considering the decomposition part. Thanks to Figs. 5j, k, l and 6j, k, l showing the absolute difference between both results, we observe that u is less noisy with our method, and 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.

Conclusion
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 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.
. 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
We set E n (h) = R 2 |u n,e (x + he 1 ) − 2u n,e (x) + u n,e (x − he 1 )| 2 dx. It satisfies E n (2h) ≤ 16E n (h), ∀h ∈ R. By using Fubini-Tonelli theorem, we have 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) (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 . 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 1 t q+1 dt, C being a constant and the last integral being convergent since q = 3, resulting in the uniform boundedness of |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) |u| q+1 du dξ ≤ C . (The con-stant 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) D (1,0) u n on Ω, and D (0,1) u n,e = (D (0,1) dx dy, with C independent of n. From [17,Lemma 4.33,p. 200], we know that  1 2 , then (u n ) strongly converges toū in C 0,λ b (Ω) and so pointwise everywhere on Ω. Thenũ =ū on Ω andũ = 0 on ∂Ω, by uniqueness of the weak limit.
By definition of the sequence (u n , g n , Q n , v 1,n , v 2,n ), we have ∀n ∈ N * , F n,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤F n,ε (ū,ḡ,Q,v 1 ,v 2 ), and by taking the lim sup when n → +∞, it yields lim sup n→+∞F n,ε (u n , g n , Q n , v 1,n , v 2,n ) ≤F ε (ū,ḡ,Q,v 1 ,v 2 ). Indeed, thanks to [16,Theorem 9], we know that (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 and the conclusion follows. As We know that (v 2 (B(0, L R)). Using the dominated convergence theorem, we get that (2,0) u(x)| 2 dx −→ n→+∞ 0, and by letting ε tend to 0, we conclude that yielding, thanks to the second triangle inequality, from Cauchy-Schwarz inequality. It follows that It results from this that: Exchanging the role of v ε and u, we obtain: Combining these different elements in (A.1) yields the desired result.

−1 β
The last inequality implies that X Y . We then use the following lemma.