Christian Wolf
email: christian.wolf@liris.cnrs.fr

Liris

An iterative graph cut optimization algorithm for a double MRF prior

Keywords: Markov Random Fields, Bayesian estimation, Graph cuts, Document Image Restoration

published or not. The documents may come

An iterative graph cut optimization algorithm for a double MRF prior

Christian Wolf

Introduction

In a previous paper [START_REF] Wolf | Document ink bleed-through removal with two hidden Markov random fields and a single observation field[END_REF] we presented a double Markov random field (MRF) model designed for document restoration, more specifically document ink bleedthrough removal, i.e the extraction and replacement of the verso ink traversing the paper and showing on the recto side. The original method used simulated anneling for optimization, which is theoretically known to converge to global solution but which is painfully slow in practice. In this paper we propose an efficient optimization technique based on graph cuts.

The paper is organized as follows: section 2 very briefly describes the graphical model, details on its motivation and its derivation can be found in [START_REF] Wolf | Document ink bleed-through removal with two hidden Markov random fields and a single observation field[END_REF]. Section 3 formulates the iterative optimization algorithm based on graph cuts. Section 4 illustrates the experiments performed to evaluate the results on a dataset of scanned documents, and section 5 finally concludes.

The graphical model

The model introduced in [START_REF] Wolf | Document ink bleed-through removal with two hidden Markov random fields and a single observation field[END_REF] consists of two hidden label fields, one for the recto side and one for the verso side, as well as a single observation field (the scanned image). We therefore consider a segmentation problem where each pixel corresponds to two different hidden labels, one for each field, and where each label is chosen from a space of two labels: text and back ground. Formally, we have a graph G = {V, E} with a set of nodes V and a set of edges E. V is partitioned into three subsets: the two fields of hidden variables F 1 and F 2 and the field of observed variables D. The three fields are indexed by the same indices corresponding to the pixels of the image, i.e. F 1 s , F 2 s and D s denote, respectively, the hidden recto label, the hidden verso label and the observation for the same pixel s. The hidden variables F 1 s and F 2 s may take values from the set Λ = {0, 1}, where 0 corresponds to background and 1 corresponds to text. The set of edges E defines the neighborhood on the graph, i.e. there is an edge between to nodes r and s if and only if r ∈ N s and s ∈ N r .

The dependency graph (see figure 1) contains the following cliques types: first order and second order "intra-field" cliques in the subgraph F 1 , first order and second order "intra-field" cliques in the subgraph F 2 (we will assume the 3-node clique potentials to be zero) and finally the "inter-field" cliques between F 1 , F 2 and D.

The joint probability distribution of the whole graph can therefore be given as follows:

P (f 1 , f 2 , d) = 1 Z exp -U (f 1 ) + U (f 2 ) + U (f 1 , f 2 , d) /T (1)
The terms U (f 1 ) and U (f 2 ) correspond to two Potts models, one prior for each field:

U (f 1 ) + U (f 2 ) = U (f 1 , f 2 ) = {s}∈C1 α 1 f 1 s + {s,s ′ }∈C2 β 1 s,s ′ δ f 1 s ,f 1 s ′ + {s}∈C1 α 2 f 2 s + {s,s ′ }∈C2 β 2 s,s ′ δ f 2 s ,f 2 s ′ (2)
where C 1 is the set of single site cliques, C 2 is the set of pair site cliques and δ is the Kronecker delta defined as δ i,j = 1 if i = j and 0 else. The term U (f 1 , f 2 , d) corresponds to the data likelihood, which factorizes as follows:

P (d|f 1 , f 2 ) = s N (d s ; µ s , Σ s ) (3)
where µ s is the mean for class f s (in the degraded image) and Σ s is the covariance matrix for class f s given as follows (note that µ r = µ r etc.):

µ s =    µ r if f 1 s = text µ v if f 1 s = background and f 2 s = text µ bg else Σ s =    Σ r if f 1 s = text Σ v if f 1 s = background and f 2 s = text Σ bg else (4
) where µ r , µ v and µ bg are, respectively, and in the degraded image, the means for the recto class, the verso class and the background class, and the covariances are denoted equivalently.

The posterior probability and its maximization with graph cuts

Applying the Bayes rule, the posterior probability of the two label fields can be given as follows:

P (f 1 , f 2 |d) = 1 P (d) P (f 1 , f 2 )P (d|f 1 , f 2 ) ∝ P (f 1 )P (f 2 )P (d|f 1 , f 2 ) (5)
As usual, we can ignore the factor 1 P (d) not depending on the hidden variables and maximize the joint probability, or minimize its energy. Combining it with the two Potts models and the data likelihood, we get the following energy potential function:

U (f 1 , f 2 , d) = {s}∈C1 α 1 f 1 s + {s,s ′ }∈C2 β 1 s,s ′ δ f 1 s ,f 1 s ′ + {s}∈C1 α 2 f 2 s + {s,s ′ }∈C2 β 2 s,s ′ δf 2 s , f 2 s ′ + {s}∈C1 1 2 (d s -µ s ) T Σ -1 s (d s -µ s ) (6)
where and µ s and Σ s are the sufficient statistics for the observation model given by the labels f s and f s ′ . To estimate the binary images, equation ( 5) must be maximized. Unfortunately, the function is not convex and standard gradient descent methods will most likely return a non global solution. Simulated Annealing has been proven to return the global optimum under certain conditions [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], but is painfully slow in practice. Loopy belief propagation is another option, giving an approximative solution by iteratively applying Pearl's belief propagation algorithm originally designed for belief networks [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems[END_REF]. In this work we will take advantage of the nature of the dependency graph (binary labels and cliques with not more than 2 hidden labels) in order to derive an optimization algorithm based on the calculation of the minimum cut/maximum flow in a graph [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF][2][3] [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF].

For convenience we will rewrite the energy function for the whole graph in terms of unary functions U 1 and two types of binary functions U 2 and U ′ 2 as follows:

U (f 1 , f 2 , d) = {s}∈C1 α 1 U 1 (f 1 s ) + α 2 U 1 (f 2 s ) + {s,s ′ }∈C2 β 1 s,s ′ U 2 (f 1 s , f 1 s ′ ) + β 2 s,s ′ U 2 (f 2 s , f 2 s ′ ) + {s}∈C1 U ′ 2 (f 1 s , f 2 s ; d s ) (7)
where

U 1 (f s ) = f s , U 2 (f s , f s ′ ) = δ fs,f s ′ and U ′ 2 (f 1 s , f 2 s ; d s ) = 1 2 (d s -µ s ) T Σ -1 s (d s -µ s )
. We consider U ′ 2 (., .; .) a binary function since we do not maximize over the third argument, which is an observed variable.

Although the problem involves two possible labels for each hidden variable (|Λ| = 2), the exact solution for equation [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems[END_REF] cannot be found using algorithms based on graph cuts. As shown by Kolmogorov et al. [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF], a function of binary variables composed of unary terms and binary terms is graph-representable, i.e. it can be minimized with algorithms based on the calculation of the maximum flow in a graph, if and only if each binary term E(., .) is regular, i.e. it satisfies the following equation:

E(0, 0) + E(1, 1) ≤ E(0, 1) + E(1, 0) (8)
It can be easily seen that this is the case of the terms U 2 (., .) in equation ( 7), but not necessarily for all terms U ′ 2 (., .; .). According to the value of the observation d s at site s, U ′ 2 (f 1 s , f 2 s ; d s ) may be regular or not. In other words, if the observation likelihood for equal labels f 1 s and f 2 s is higher than the observation likelihood for different labels, then the term is regular for site s.

We therefore propose an adaptation and extension of the iterative α-expansion move algorithm proposed by Boykov et al. [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] for labeling problems with multiple labels (|Λ| > 2) and improved by Kolmogorov et al. [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. In the original formulation for multi label problems, each subproblem is a binary problem where each hidden variable may take two virtual labels: x s and α, where x s is the original (current) label, and α is a new label, whose value is changed at each iteration.

In our case, the iteratively solved binary labeled and regular subproblems arise fixing the hidden label of one of the two fields F 1 and F 2 and estimating the labels of the other one. Completely fixing a whole set of variables corresponds to running an α-expansion move algorithm on a single field dependency graph where each single hidden variable f s my take 4 values (background, recto, verso, recto-verso) and the pairwise clique potentials are adapted accordingly. This optimization schedule may be improved by fixing only the variables whose site s is not regular, and jointly estimate the variables f 1 s and f 2 s for the regular sites s. For convenience we introduce a binary matrix H indicating for each site s whether it is regular or not, i.e. whether the associated function Figure 2: The inference algorithm iteratively optimizing two different binary subproblems. H is a matrix storing for each pixel whether its likelihood term is regular or not. U →2 (f 1 , f 2 , d, H) and U →1 (f 1 , f 2 , d, H) correspond to the posterior energy with different hidden variables clamped.

Input: d (a realization of the observed field) Output: f 1 , f 2 (estimated label fields) F 1 , F 2 ← Initialize the label fields (e.g. with k-means) H ← determine the regular sites s repeat -Fix f 1 s for all s and f 2 s for H s = 0, estimate optimal f 2 for all s and f 2 s for

H s = 1 maximizing U →2 (f 1 , f 2 , d, H) -Fix f 2
s for all s and f 1 s for H s = 0, estimate optimal f 1 for all s and f 1 s for

H s = 1 maximizing U →1 (f 1 , f 2 , d, H) until Convergence U ′ 2 (f 1 s , f 2 s ; d s ) is regular or not: H s =    1 if U ′ 2 (0, 0, d s ) + U ′ 2 (1, 1, d s ) ≤ U ′ 2 (0, 1, d s ) + U ′ 2 (1, 0, d s ) 0 else (9)
Figure 2 outlines the inference algorithm, which iteratively calculates the exact solution of two different binary subproblems, maximizing, respectively,

U →2 (f 1 , f 2 , d, H) and U →1 (f 1 , f 2 , d, H
). These two energy functions are actually equivalent, however, the set of fixed variables and the set of estimated variables being different, they lead to two different cut graphs. In order to show the derivation of the cut graph, we will rewrite the two functions by reordering some terms. W.l.o.g., in the rest of this section we describe U →2 (f 1 , f 2 , d, H), i.e. the subproblem where a subset of the variables in F 1 is fixed, whereas the variables of F 2 and the complementary subset of F 1 are estimated. The function U →1 (f 1 , f 2 , d, H) corresponding to the complementary subproblem can be derived in similar way

After separating terms according to the contents of H, the corresponding energy function can be given as follows:

U →2 (f 1 , f 2 , d, H) 1 = {s}∈C1:Hs=0 α 1 U 1 (f 1 s ) 2 + {s}∈C1:Hs=1 α 1 U 1 (f 1 s ) 3 + {s}∈C1 α 2 U 1 (f 2 s ) 4 + {s,s ′ }∈C2:Hs=0∧H s ′ =0 β 1 s,s ′ U 2 (f 1 s , f 1 s ′ ) 5 + {s,s ′ }∈C2:Hs=1∧H s ′ =1 β 1 s,s ′ U 2 (f 1 s , f 1 s ′ ) 6 + {s,s ′ }∈C2:Hs =H s ′ β 1 s,s ′ U 2 (f 1 s , f 1 s ′ ) 7 + {s,s ′ }∈C2 β 2 s,s ′ U 2 (f 2 s , f 2 s ′ ) 8 + {s}∈C1:Hs=0 U ′ 2 (f 1 s , f 2 s ; d s ) 9 + {s}∈C1:Hs=1 U ′ 2 (f 1 s , f 2 s ; d s ) (10
) Written in this notation, The energy functions can be directly translated into a cut graph using the method introduced by Kolmogorov et al. [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF]. The cut graph then contains, besides the terminal nodes source and sink, one node for each variable F2 s as well as one node for each variable F1 s satisfying H s = 1. Each unary term is translated into a tedge, and each binary term is translated into an n-edge as well as two t-edges.

The terms in lines 1 and 4 of equation [START_REF] Tonazzini | A markov model for blind image separation by a mean-field em algorithm[END_REF] do not depend on estimated variables and therefore can be omitted during the minimization. The terms in lines 2 and 3 contain standard unary functions and will be represented by t-edges. The terms in lines 5 and 7 contain standard binary functions (pairwise cliques of the Potts model) and will be represented by n-edges. The terms in line 6 are binary functions (also pairwise cliques of the Potts model) in the full original expression (equation ( 7)), but one of the two arguments is fixed in equation [START_REF] Tonazzini | A markov model for blind image separation by a mean-field em algorithm[END_REF] describing the sub problem. They can therefore be represented as t-edges in the cut graph. Similarly, the terms in line 8 are non-regular pairwise functions of the observation model, which can be represented as t-edges. The terms in line 9, finally, correspond to the regular pairwise function of the observation model, which can be represented as nedges.

Table 1 gives a full description of the different edges of the cut graph and their weights. Figure 3 shows an example of a dependency graph for a toy problem, a 3 × 1 image, and two different cut graphs. Figure 3b shows the cut graph for the αexpansion move like algorithm, i.e. all sites s are considered as non-regular. The cut graph is shown for the case where the complete set of variables F 1 is fixed whereas the complete set of variables F 2 is estimated.

Figure 3c shows the extended algorithm, where the middle and the right site are considered regular, whereas the left site is considered non-regular. For the middle and the right site, the variables F 1 s and F 2 s are jointly estimated, whereas for the left site only F 2 s is estimated whereas F 1 s is fixed.

Experimental results

Evaluating document restoration algorithms is a non trivial task since ground truth is very hard to come by. Short of manually classifying each pixel in a scanned image, the only way to get reliable ground truth data on pixel level is to test the algorithm on synthetic data. These tests, on the other hand, may not be realistic enough to capture all the subtleties of a real environment. To evaluate our algorithm we therefore decided to test its ability to improve the performance of an OCR algorithm when applied to real scanned documents. We chose a dataset consisting of 104 pages of low quality printed text from the 18 th century, the Gazettes de Leyde. This journal in French language was printed from 1679 to 1798 in the Netherlands in order to escape the censorship in France at the 18 th century and relates news of the world. The Gazettes are currently used by several research projects in social and political sciences, some of which are currently collaborating with our team in the framework of digitization projects.

From an image processing point the view, the data situates itself between the difficulty of manuscripts and the regularity of printed documents. The images of sizes around 1030×1550 pixels are of very low quality compared to modern printed text. Recognition is possible, although the performance on the non-restored images is not very high. We chose the open source OCR software Tesseract published by Google 1 for our experiences, mainly because it is easily scriptable. We performed some selected experiments with the product of the market leader, Abby Finereader 8 2 , which performs slightly better without changing the ranking of the restoration methods.

We compared the proposed method with several competing methods. One group of algorithms purely exploits the fact that, according to the hypothesis that set recto pixels completely cover verso pixels, without taking into account interactions between neighboring pixels. Examples are the k-n-edges for node pairs:

Weight

Line in eq. ( 10)

F 2 s , F 2 s ′ : (s, s ′ ) ∈ C 2 -β 2 s,s ′ 7 F 1 s , F 1 s ′ : (s, s ′ ) ∈ C 2 ∧ H s = 1 ∧ H s ′ = 1 -β 1 s,s ′ 5 F 1 s , F 2 s : H s = 1 U ′ 2 (0, 1, d s ) + U ′ 2 (1, 0, d s ) -U ′ 2 (0, 0, d s ) -U ′ 2 (1, 1, d s ) 9 (a)
t-edges (to source if weight> 0) for nodes:

Weight

Line in eq. ( 10)

F 2 s α 2 3 F 2 s : H s = 0 U ′ 2 (f 1 s , 1, d s ) -U ′ 2 (f 1 s , 0, d s ) 8 F 1 s : H s = 1 α 1 2 F 1 s : H s = 1 s ′ :H s ′ =0∧f 1 s ′ =1 β 1 s,s ′ - s ′ :H s ′ =0∧f 1 s ′ =0 β 1 s,s ′ 6 F 1 s : H s = 1 U ′ 2 (1, 0, d s ) -U ′ 2 (0, 0, d s ) 9 F 2 s : H s = 1 U ′ 2 (1, 1, d s ) -U ′ 2 (1, 0, d s ) 9 (b)
Table 1: The edges added to the cut graph for the proposed inference algorithm: each edge corresponds to a term in eq. ( 10). Each t-edge is connected to the source if the weight is positive, or connected to the sink if the weight is negative, in which case the absolute value of the weight is used. Multiple edges between same nodes (taking into account the orientation) are replaced by a single edge, its weight being the sum of the individual weights. means clustering algorithm with k=3 clusters (followed by our restoration algorithm replacing verso pixels, explained in [START_REF] Wolf | Document ink bleed-through removal with two hidden Markov random fields and a single observation field[END_REF], as well as two thresholding algorithms. We chose two methods which represent the state of the art in adaptive thresholding: Niblack's algorithm [START_REF] Niblack | An Introduction to Digital Image Processing[END_REF] which performed best in a widely cited evaluation paper [START_REF] Trier | Goal-Directed Evaluation of Binarization Methods[END_REF] as well as an improvement of Niblack's algorithm by Sauvola et al. [START_REF] Sauvola | Adaptive Document Binarization[END_REF]. Since a restoration is not straightforward from a binary output, we directly fed the binary images to the OCR in the case of the two thresholding algorithms.

As mentioned in section 1, statistical source separation is one of the most active areas in bleedthrough removal with several works published by Tonazzini et al. on this subject [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF][10] [START_REF] Tonazzini | Bayesian MRFbased blind source separation of convolutive mixtures of images[END_REF] [START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF]. We therefore decided to compare the proposed method with two of them: since the scans of the Gazettes de Leyde are in color, the color model introduced in [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF] and which we described in section 1 is applicable. The second method, introduced in [START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF] and based on orthogonalization, is non-blind and therefore requires the presence of the verso side of the image. In a personal communication sent for the experiments in this paper, Prof. Tonazzini recommended the use of two different planes of the color image as recto and verso observations, which we did in our experiments. The source codes have kindly been provided by the author, Prof. Tonazzini itself.

The tested source separation methods are not automatic, they need user interaction in order to chose the correct output source plane. While the number of the correct recto plane may be different between different images, tests showed that for all 104 images of the Gazettes de Leyd, the order of the source images was the same. The source planes resembling the most to the assumed recto plane where, for both methods, source #1 and source #2, which we both included into the experiments. This was not the case for other images, as for instance the manuscripts shown in figures 7 and8.

The last method compared to the proposed algorithm is a standard single MRF with a Potts model and three labels (recto, verso and background ), optimized using Kolmogorov et al's version of the αexpansion move algorithm [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF] and combined with the same parameter estimation and pre-and postprocessing as our proposed method.

Figures 4 to 6 illustrate the OCR results on a small image taken from the Gazettes dataset. As we can see, being based on segmentation, the results for k-means and the two MRF methods are similar. The k-means result (Figure 4b) is noisy as opposed to the MRF results, the double MRF (Fig- ure 4d) improves the regularity of the single MRF (Figure 4c). The OCR output is a little bit cleaner for the double MRF case.

The results of Niblack's algorithm and Sauvola et al.'s algorithm show the typical weaknesses of these approaches: Niblack (Figure 5a) produces spurious components, especially in areas with few text, and Sauvola (Figure 5b) tends to cutting characters into

D F 2 F 1 F 2 source(0/non -text) sink(1/text) F 2 F 1 source(0/non -text) sink(1/text) Figure 3: (a)
The dependency graph of a simple model containing three pixels in a single row; (b) the cut graph for an α-expansion move like inference algorithm: inference of the verso pixels; (c) the cut graph for the proposed inference algorithm: joint inference of the verso pixels and of a subset of the recto pixels. In this exemple, the potential functions related to the observation model are regular for the middle and for the right pixel (H s = 1), but not for the left one (H s = 0). 83.23 74.85 32,537 100 † Not available: lack of OCR performance makes a correct evaluation impossible ‡ results obtained with a subset of the images only (absolute cost is not comparable).

Method

Table 2: OCR results on a database of 104 scanned document images: non-restored input images and different restoration methods. several parts due to its assumptions on the grayvalue distribution in the image.

Figures 5c and5d show the first two source components of the non-blind source separation method [START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF] applied to the color components red and green of the color input image. All source separation results are shown without the post processing recommended by the authors (see below).

The second, blind method [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF], shown in Figures 6a-c, delievers similar results: although we can identify a source component which does not include the verso text, the response itself is quite noisy and faint. Post-processing the image slightly improves the latter but tends to increase the noise. Figure 6d shows an image which corresponds to a grayscale conversion of a color image composed of the three different source components obtained with the color based method [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF]. Although this result was not intended, as the verso component is still part of the image, the result seems to be better than the ones consisting of a single source component only. Surprisingly, this result is the only one which produces at least limited OCR output, whereas the other images do not produce anything meaningful.

In order to evaluate the amount of recognition improvement of the restoration method, we manually created groundtruth for the 104 images, and calculated the Levenstein edit distance between two strings [START_REF] Wagner | The string to string correction problem[END_REF], which finds the optimal transformation from on string into another with elementary operations (insertion, deletion, substitution) minimizing the global cost of these operations. Additionally, we calculated character recall and character precision derived from the transformation operation of this distance. Table 2 compares the measures for the different methods described above, as well as the recognition performance on not restored images. Note, that precision and recall are indepen- dant of the dataset size, whereas the total transformation cost is not. We can see that all methods based on identifying the verso component (k-means and the two MRF methods, including the proposed one) are capable of significantly improving the recognition results compared to no restoration at all. Not surprisingly, regularizing the segmentation with a priori knowledge boosts the performance. Separating the regularization of the recto and verso side further improves recognition, gaining 1.2 percent points in recall compared to the single MRF and 2.7 percentage points in precision. Totally, compared to no restoration at all, the proposed method improves recognition at about 17 percentage points in terms of recall and around 25 percentage points in terms of precision.

Recognition on the results of Niblack's method produces only gibberish, probably because of the small ghost objects it creates. Sauvola et al.'s method overcomes this problem and the recognition performance almost attains the quality of the three class segmentation of performed by the k-means algorithm.

Surprisingly, the recognition performance on the results of the two source separation results was very disappointing. We performed recognition experiments for both planes of the first method [START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF] and all three planes of the second method [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF], respecting the author's recommendations to darken the images after applying the inverted mixture matrix. In a personal communication for the experiments in this paper, Prof. Tonazzini recommended subtracting the K component of the CMYK color decomposition. However, we obtained better results with a histogram stretch instead of the proposed method.

Unfortunately, the recognition performance on these results was not good enough to include it in the table. Most of the output was blank or gibberish, making an evaluation impossible. We managed to get some statistics on the first source plane of the first method, as well as on output images combining all three source planes of the second method. However, this was only possible when a subset of the dataset was removed. Even then, the results where not competitive.

Figures 7 and8 show restoration results on two different manuscript images. The source separation methods remove more of the verso text in Figure 7, but unfortunately the contrast is very low and they are significantly disturbed by the JPEG artifacts in the input image. The performance shown in figure 8 by the source separation methods are more spread out across the image and seem to touch more of the low frequency components.

Computational complexity

The computational complexity of the proposed method is dominated by the inference part with minimum cut/maximum flow whose complexity is bounded by O(|E| * f ), where |E| is the number of edges in the graph and f is the maximum flow. We use the graph cut implementation by Boykov and Kolmogorov [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF] which has been optimized for typical graph structures encountered in computer vision and whose running time is nearly linear in running time in practice [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. Table 3 gives effective run times measured on a computer equipped with an Intel Core 2 processor running at 2.5Ghz and 4GB of RAM (only one core was used). The running time of the proposed method is comparable to the running time of a single MRF with graph cut optimization and quite competitive given its restoration performance.

Conclusion and Outlook

In this paper we presented a method to separate the verso side from the recto side of a single scan of document images. The novelty of the method is the separation of the MRF prior into two different label fields, each of which regularizes one of the two sides of the document. This separation allows to estimate the verso pixels of the document which are covered by the recto pixels, which, again through the MRF prior, improves the estimation of the verso pixels not covered by recto pixels, thus increasing the performance of the regularization. We showed that this formulation leads to an efficient algorithm based on graph cuts.

The performance of the method has been evaluated on scanned document images from the 18 th century, showing that the restoration is able to improve the recognition performance of an OCR significantly, compared to non restored images but also compared to competing methods.

F 1 F 2 DFigure 1 :

 121 Figure 1: The dependency graph for a 2 × 2 pixel image. The model consists of the two label fields F 1 and F 2 ("empty" nodes) and the single observation field D (shaded nodes).

  bras 'd5:glgrComprzgnîc': |•Y• ;]L\ët4ôi;n àccoâniàagnë Lié ibn _ Uncle } ]F1'èré'dc"i'a""M ÈrèQ_ Lei j,cxn,c> Pxjnèci ziypir . four motif dc fa fqizc lc déiïr de/trouver ies mbicns i d'é"rë;ifèi' Qin Përç dg triiïcfgitar,. 0ùlî1.é;0it Yë(gliE1_'• bras de lx =>, Compsignielz émir accompagné de ibn _ Oncle ,'• Frère dela Mere, Le j,eune• Prince avoit . pour motif de la faire le delîr de trouver les moyens irer fh;] Père du t.tiPce état,. où. il, étoit réduit; bras de l:îC0mP;rg1rrE: Il étoit accompagné dé (on _ Oncle ,' Frère de fu Mere. Le ierzne Prince avoit . pour motif de la fuite le delîr de trouver les mbvens 'dë'r•• iter *}>;> Père du ttiftc état,. où il. étoit réduitr; ' bras de l:1Comp;1gni•:: Il étoit accompagne dé Ton _ ' Oncle, Frère de fa Mere. Le j.e;:ne• Prince avoit pour motif dc fa fuite le defir de trouver les moyens 'àië'r<<• iter '\}>>:• Père du triûe état,. ou il étoit réduit.:

Figure 4 :

 4 Figure 4: Small extracts of the OCR results obtained on scanned document images: (a) input image (no restoration); (b) k-means segmentation + restoration; (c) single MRF segmentation [5] + restoration; (d) double MRF (proposed method).

Figure 5 :

 5 Figure 5: Small extracts of the OCR results obtained on scanned document images: (a) segmentation with Niblack [6]; (b) segmentation with Sauvola et al. [8]; (c) Tonazzini et al. [12] source #1, Tonazzini et al. [12] source #2

Figure 6 :

 6 Figure 6: Small extracts of the OCR results obtained on scanned document images: (a) Tonazzini et al. [9] source #1 (b) Tonazzini et al. [9] source #2 (c) Tonazzini et al. [9] source #3 (d) Tonazzini et al. [9] all three sources combined.

Figure 7 :

 7 Figure 7: Restoration results on real data. From left to right, top to bottom: input image, k-means, single MRF & α-exp. move [5], double MRF (proposed method), Tonazzini et al. [12] source #1, Tonazzini et al. [12] source #2, Tonazzini et al. [9] source #1, Tonazzini et al. [9] source #2 (source #3 not displayed).

Figure 8 :

 8 Figure 8: Restoration results on real data. From left to right, top to bottom: input image, k-means, single MRF & α-exp. move [5], double MRF (proposed method), Tonazzini et al. [12] source #2, Tonazzini et al. [12] source #1, Tonazzini et al. [9] source #3, Tonazzini et al. [9] source #1 (source #2 not displayed).

Table 3 :

 3 reveils similar strengths and weaknesses, typical to the two types of approaches: the regularized segmentation approaches create crisp images but show localized artificacts, whereas the artificats created Execution times in seconds for various methods. The replacement of the background pixels is included if the method is based on segmentation.

	Type	Method	1026	2436
			×1587 ×3320
	Context-	Niblack [6]	0.6	2.9 †
	free	Sauvola et al. [8]	0.6	2.9 †
		K-Means (k=3)	1.9	10.5 †
	Source-	Tonazzini et al. [9]	36.9	134.6 ‡
	separation Tonazzini et al. [12]	17.0	74.5 ‡
	MRF	Single MRF [5]	7.2	34.9 †
		Double MRF	7.4	36.4 †
		† Code in C++
		‡ Code in Matlab/GNU Octave

http://code.google.com/p/tesseract-ocr

http://finereader.abbyy.com

Acknowledgements

We thank Anna Tonazzini for providing us with the source code of the two source separation methods and her kind help in setting up the corresponding experiments as well as for the interesting discussions.


Ben Hachmi
email: hachmi.ben-dhia@ecp.fr

Olivier Dhia
email: olivier.jamond@ecp.fr

Jamond

Ben Dhia

On the use of XFEM within the Arlequin Framework for the Simulation of Crack Propagation

Keywords: multimodel, Arlequin method, XFEM, crack propagation

The Arlequin method is a generic numerical method that allows, by local superposition and coupling of models, to address multimodel and multiscale mechanical problems. In particular, this method has already been used to super-impose cracked patches on sound structures, reducing this way the global simulation ressources a classical finite element approach would have required.

In this paper, one of the key features of the Extended-Finite Element Method, namely the Heaviside enrichment function, is used within the Arlequin framework to further reduce the costs of crack propagation simulations. The main goal of the paper is to describe the proposed methodology and to assess its performance through numerical experiments.

Introduction

Generation of classical regular finite element meshes that are capable of properly representing material discontinuities or accurately simulating the propagation of cracks in complex structures with high strain and stress localization is a tedious and time consuming endeavor. This explains the search over the last two decades for flexible numerical methodologies in order to efficiently tackle this difficult issue. Among others, methods that have been proposed are the so-called meshless methods, see e.g. [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF][START_REF] Belytschko | Element-free Galerkin methods[END_REF], the partition-of-unity method (PUM) [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF][START_REF] Babuška | The partition of unity method[END_REF] and its well-known instances, namely the generalized finite element method (GFEM) (e.g. [START_REF] Babuška | The design and analysis of the Generalized Finite Element Method[END_REF]) and the extended finite element method (XFEM), see e.g. [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Dolbow | Discontinuous enrichment in finite element with a partition of unity method[END_REF]). In these methods, somewhat reminiscent of the Treftz approach, the finite element space associated with a given coarse mesh is enriched by analytical or semi-analytical functions based on a priori knowledge of the unknown solution. In the XFEM approach, for example, one uses level set functions and asymptotic local modal functions around the cracks. XFEM/GFEM methods are particularly relevant and effective whenever the local behavior of the solution is known. Another class of approaches is given by the s-method, which is in many ways similar to the Fast Adaptive Composite-grid (FAC) method [START_REF] Mccormick | The fast adaptive composite-grid (FAC) method for elliptic equations[END_REF]. This approach can be considered as a genuine finite element method in which the coarse finite element space is locally enriched by fine scale basis functions associated with a fine mesh. Although effective in many regards, the method may suffer from redundancy. Moreover, the resulting problem in the s-method is solved globally. An approach similar to the s-method is one that involves the use of local meshes (also called patches), but whose resulting system of equations is solved by a relaxed Chimera algorithm [START_REF] Glowinski | Approximation of multi-scale elliptic problems using patches of finite elements[END_REF][START_REF] Glowinski | Finite element approximation of multi-scale elliptic problems using patches of elements[END_REF] (see also the Schwarz algorithm (see e.g [START_REF] Steger | On the use of composite grid schemes in computational aerodynamics[END_REF])). Our ultimate goal being the simulation of physically-based damage, localisations and crack propagation, we propose a novel and alternative approach, experienced in this paper on basic tests. The approach is based on the multimodel and multiscale Arlequin framework [START_REF] Dhia | Multiscale mechanical problems: the Arlequin method[END_REF][START_REF] Ben Dhia | Application of the Arlequin method to some structures with defects[END_REF][START_REF] Rateau | Méthode Arlequin pour les Problèmes Mécaniques Multi-échelles. Application à des Problèmes de Jonction et de Fissuration de Structures Élancées[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Global-local approaches: the Arlequin framework[END_REF][START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF]. The method essentially consists in superimposing a fine scale model (in a patch) to a large scale model combined together using a partition of energy. The construction of the resulting coupled model allows for great flexibility; indeed, one may independently refine the mesh or use different constitutive laws, governing equations, or geometries within the patch. For instance, the Arlequin framework has been recently employed for successfully coupling atomistic models with continuum models [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF][START_REF] Ben Dhia | Coupling of atomistic and continuum models in the Arlequin framework[END_REF][START_REF] Bauman | On the application of the Arlequin method to the coupling of particle and continuum models[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF]. The idea of using the Arlequin framework for the simulation of crack propagation, fatigue, or damage phenomena in a structure, has already been advanced a few years ago [START_REF] Rateau | Méthode Arlequin pour la propagation de défauts dans des pneus : Etude de Faisabilité[END_REF][START_REF] Dhia | The Arlequin method: a partition of models for concurrent multiscale analyses[END_REF][START_REF] Romdhane | [END_REF]. In order to expose our motivation and put in evidence our new numerical achievements, we report here some of our earlier results. In the first example (from [START_REF] Rateau | Méthode Arlequin pour la propagation de défauts dans des pneus : Etude de Faisabilité[END_REF]), we considered a solid composed of two materials, one soft (rubber) and one hard (steel). The structure was submitted to a cyclic loading on the top surface, to symmetry conditions along the bottom surface, while the other two boundaries were kept free. Figure 1 shows the evolution of damage as well as nucleation and propagation of macroscopic cracks whenever damage exceeded a prescribed threshold. The fine scale phenomena near the crack tip were taken into account using an Arlequin patch which was free to move and grow according to the evolution of the crack tip. In this case, the simulation was still relatively simple due to the fact that the direction of propagation was a priori known to be the straight bi-material interface. The second example, taken from [START_REF] Romdhane | [END_REF], illustrates the crack propagation in a linearly elastic structure as shown in Figure 2. The patch consisted of a very fine mesh at the tip and was progressively adapted following the evolution of the crack. We note that the underlying coarse mesh was kept unchanged during the whole process. Data for this test case (referred to as the double cantilever beam test) was taken from Belytschko and Black [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF]. The obvious drawback of the method was that the number of finite elements along the lips of the crack was far from being optimal, which unnecessarily increased the cost of the computation. Moreover, the treatment of more than one crack would quickly become an intricate issue.

In this work, we propose to circumvent the aforementioned issues by intro-ducing Heaviside enrichment functions, as in XFEM, within our Arlequin framework, with the objective of further reducing computational costs and enhancing the flexibility of crack propagation simulations in this framework.

In particular, we will show that the combined approach allows us to consider smaller patches to capture the local features near the tip since Heaviside functions are used to enrich the coarse model. The main goal of the present paper is to describe the proposed methodology and to assess its performance through numerical experiments.

The paper is organized as follows. Following the introduction, we provide in Section 2 the continuous formulation of the crack problem in a linearly elastic structure. In Section 3, we recall the main features of the extended finite element method. We then show in Section 4 how enrichment functions a la XFEM can be incorporated into our Arlequin framework. In particular, we describe the continuous Arlequin formulation for the crack problem in Subsection 4.1. The discrete counterpart is given in Subsection 4.2 and practical aspects related to numerical integration are detailed in Subsection 4.3. Numerical results in order to assess the performance of our approach are then presented in Section 5. We use to that effect two examples, one dealing with the simulation of a kinked crack and the double-cantilever beam test case. In these two examples, the patch is strictly contained in the interior of the global domain (no intersection between the boundary of the patch and the global domain boundary). A third numerical example shows the feasibility of the simulation of the propagation of a crack from a side of a square plate to the opposite one. This example shows, in particular, how the patch is adapted easily when it hurts the boundary of the global domain. The paper closes with concluding remarks and some ongoing works in Section 6.

Formulation of a cracked elastic body problem

Let us consider a linear elastic cracked solid, occupying the closure of a bounded domain Ω 0 . For the sake of clarity, we will consider the bidimensional case. The boundary of the solid is partitioned into Γ u , Γ g , Γ + and Γ -, where Γ u is the part where the solid is clamped and Γ g is the part on which an applied density of loads g is given. The other two parts of the boundary Γ + and Γ -define the faces of the crack whose tip is denoted by o (see figure 3). It is assumed that the crack is closed in the undeformed configuration of the solid. Thus, from a geometrical point of view, we have Γ + ≡ Γ -≡ Γ. But free boundary conditions are applied on Γ + and Γ -in the remainder of the article. That is to mean that no unilateral contact loads are considered between the faces of the crack. Moreover, for the sake of conciseness, we neglect external body forces. The classical weak form of the linearized primal elastic problem reads:

Ω 0 Γ o Γ u Γ + Γ -
Fig. 3. problem of a cracked linear elastic solid

find u ∈ V 0 ; ∀v ∈ V 0 a(u, v) = l(v) (1)
where:

V 0 = v ∈ H 1 (Ω 0 ); v = 0 on Γ u (2) a(u, v) = Ω 0 σ(u) : (v) dx (3)
l(v) = Γg g • v dγ (4)
with σ and standing for the stress and linearized strain tensors, the stess tensor being related to (u) = ∇ s (u) by the Hooke's law:

σ ij (u) = R ijkl kl (u) in Ω 0 (5)
where the Einstein's convention of summation on repeated indices is used and where the elasticity tensor R is assumed to verify:

R ijkl = R jikl = R klij (symmetry) (6) ∃c 0 ; R ijkl τ ij τ kl ≥ c 0 τ ij τ ij , ∀τ ij = τ ji (coercivity) (7) R ijkl ∈ L ∞ (Ω 0 ) (regularity) (8)
For a classical conforming finite element approximation of this problem, the mesh has to match the crack. But since the latter may propagate into the solid without necessarily fiting the existing mesh edges, one has to generate repetively new meshes compatible with the evolution of the crack in order to achieve an accurate numerical approximation of the problem solution. This procedure is clearly unflexible and costly. The Extended Finite Element Method is one of the most efficient methods that has been designed to address this issue. Since the XFEM is partly used in our approach, its main features are recalled in the next section where some notations are introduced.

XFEM discretisation

The XFEM discretization of the problem (1)-( 5) consits in constructing a finite dimensional subspace of V 0 by enrichment of a classical finite element subspace of

V Γ 0 = v ∈ H 1 (Ω Γ 0 ); v = 0 in Γ u (9)
where Ω Γ 0 = Ω 0 ∪ Γ denotes the uncracked domain. This is done as follows. Let T h be a mesh of the uncracked solid whose edges do not necessarily match Γ. Let (ϕ 1 , ϕ 2 , ..., ϕ n ) denote a basis spanning a classical Lagrange finite element space V Γ 0h , associated to T h (see e.g. [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF][START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]) and assumed to be a subspace of V Γ 0 . The compact support of a shape function ϕ i is denoted by ω i and its interior by

• ω i . The set of indices k such that o ∈ ω k is denoted by K s (K s = {7, 8, 9
, 10} in the scalar case depicted by figure 4). The set of indices k such that k / ∈ K s and 4). The union of ω k , k ∈ K d , defines a domain we denote by Ω d . The latter is divided by the crack Γ into two domains, denoted by Ω + d and Ω - d . A precise definition of the domain Ω + d (and thus of Ω - d ) relies on the signed distance function: by defining at each point p ∈ Γ a classical positively oriented frame (τ p , ν p ) where τ p and ν p refer to unit tangent and normal vectors to the crack at point p (see figure 5), one can associate to each point x ∈ Ω d (one of) its nearest point(s) x * on Γ (as in contact mechanics, e.g. [START_REF] Ben Dhia | Level-Sets, placement and velocity-based lagrangian formulation of impact problems[END_REF]). The domain Ω + d is then defined by: (see figure 5)

• ω k ∩ Γ = ∅ is denoted by K d (K d = {1, 2, 3, 4, 5, 6} in figure
Ω + d = {x ∈ Ω d ; (x -x * ).ν x * ≥ 0} (10)
The XFEM finite dimensional space is then defined by: zone enriched with Heaviside function zone enriched with singular modes

ϕ 1 ϕ 2 ϕ 3 ϕ 7 ϕ 8 ϕ 4 ϕ 5 ϕ 6 ϕ 9 ϕ 10 Γ o ϕ 12 ϕ 11 ϕ ... Fig. 4. XFEM discretization V 0h = span (ϕ i ) i=1,n , (Hϕ j ) j∈K d , (F k ϕ l ) k=1,4; l∈K s (11) ϕ 1 ϕ 2 ϕ 3 ϕ 7 ϕ 8 ϕ 4 ϕ 5 ϕ 6 ϕ 9 ϕ 10 r Ω + d Ω - d o p τ p ν p θ Fig. 5. local frame
where H refers to the sign-like field, defined by:

H(x) =    +1 in Ω + d -1 elsewhere (12)
and where the functions F j are given by:(see e.g. [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Dolbow | Discontinuous enrichment in finite element with a partition of unity method[END_REF][START_REF] Laborde | High order extended finite element method for cracked domains[END_REF][START_REF] Chahine | Crack tip enrichment in the XFEM method using a cut-off function[END_REF], for details and more general considerations)

{F j (r, θ), j = 1, 4} = √ r sin θ 2 , √ r cos θ 2 , √ r sin θ 2 sin θ, √ r cos θ 2 sin θ (13
) with (r, θ) referring to the polar coordinates in the local frame (o, τ o , ν o ) (see figure 5). Observe here that the enriched XFEM solution is defined by:

u h0 = i=1,n u i φ i + j∈K d a j Hφ j + k=1,4 l∈K s b k l F k φ l (14)
where (u i ) i=1,n , (a j ) j∈K d and (b k l ) k=1,4; l∈K s are the unknown degrees of freedom.

To summarize, for the simulation of a cracked solid, XFEM is based on a double enrichment of classical finite element spaces by:

• a localized sign-like field that reproduces the material discontinuity due to the crack; • localized asymptotic modes in the near crack tip zone.

The first enrichment could be considered as universal. In the contrary, for general mechanical problems, the asymptotic knowledge of the behaviour of the solution in critical zones is unfortunately unavailable and one has to resort to numerical approximations in these zones. For the numerical treatment of these general problems, we suggest here to use only the first enrichment within the Arlequin framework.

Use of XFEM within the Arlequin framework

In this section, by using the Arlequin framework, a local model is superimposed to the global one in the neighbourhood of the crack tip (see figure 6). The two models are combined together using a partition of energy and a coupling operator, leading to a continuous Arlequin formulation of the problem, detailed in Subsection 4.1. A first key (modeling) feature of this continuous formulation is that it al-

Γ o Ω 1 Ω 0
Fig. 6. Arlequin domains lows us, by means of different approximations of the super-imposed models, to elaborate multimodel (and also multiscale) discrete formulations. The one elaborated in Subsection 4.2 uses: (see figure 7)

• the XFEM to enrich a coarse global finite element approximation of the global model; • a fine finite element space for the approximation of the local model. (Observe here that discrete models such as atomistic ones can also be used locally instead of the continuum elasticity model [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF][START_REF] Ben Dhia | Coupling of atomistic and continuum models in the Arlequin framework[END_REF][START_REF] Bauman | On the application of the Arlequin method to the coupling of particle and continuum models[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF], but these topics are not considered in this paper).

Fig. 7. Examples of meshes for the two arlequin models

A second key (computational) feature of the Arlequin formulation is related to its substructuring character which allows for the use of domain decompositionlike solvers (see e.g. [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Ben Dhia | Multi-modeling of multi-alterated structures in the Arlequin framework. Solution with a domain decomposition solver[END_REF]).

Continuous Arlequin Formulation

Let Ω 1 be a (copy of a) non-zero measured given regular subdomain of Ω 0 , located in the near crack tip zone. We assume for clarity that Ω 1 is strictly embedded in Ω 0 (see figure 6). This subdomain is partitioned into two regular non overlapping domains Ω c and Ω f such that (see figure 8):

o ∈ Ω f (15) ∂Ω 1 ⊂ ∂Ω c (16) meas(Ω i ) > 0, i = c, f (17)
Ω 0 Γ Ω 1 Ω f Ω c

Fig. 8. Arlequin domains

The domains Ω c and Ω f are the models coupling and free zones. In the former, the models are coupled; so that they have to be quite similar. In the latter, the models are super-imposed but free from each other; so that they can potentially be very different.

A Lagrangian Arlequin formulation of our problem, using an energy scalar product as coupling operator, reads as following: (see e.g. [START_REF] Ben Dhia | Application of the Arlequin method to some structures with defects[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF])

Find (u 0 , u 1 , Φ) ∈ V 0 × V 1 × M ; (18) ∀ v 0 ∈ V 0 , a 0 (u 0 , v 0 ) + c(Φ, v 0 ) = l 0 (v 0 ) (19) ∀ v 1 ∈ V 1 , a 1 (u 1 , v 1 ) -c(Φ, v 1 ) = 0 (20) ∀ Ψ ∈ M , c(Ψ, u 0 -u 1 ) = 0 ( 21
)
where V 0 is defined by ( 2) and where V 1 = H 1 (Ω 1 ) and the mediator space M = H 1 (Ω c ). Furthermore, for i = 0, 1, the weighted internal and external virtual works are defined by:

∀ u i , v i ∈ V i , a i (u i , v i ) = Ω i α i σ(u i ) : ε(v i ) dx (22) ∀ v 0 ∈ V 0 , l 0 (v 0 ) = Γn g.v 0 dγ (23)
The internal energy weight parameter functions (α i ), i = 0, 1, are defined in th whole domain Ω 0 and satisfy: ( [START_REF] Dhia | Global-local approaches: the Arlequin framework[END_REF][START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF])

α i ∈ [0, 1] in Ω 0 (24) α 0 + α 1 = 1 in Ω 0 ( 25
) α 0 = 1 in Ω 0 \Ω 1 (26) ∃ α 0 > 0; α i ≥ α 0 in Ω f (27)
The coupling operator c is defined by:

∀(ψ , v) ∈ M × V 1 , c(ψ, v) = Ωc κ 0 {ψ • v} + κ 1 { (ψ) : (v)} dx (28)
The positive parameters κ 0 and κ 1 are chosen such that the two terms defining c are homogeneous. In the numerical section, we have taken, as in [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF], κ 0 = 1 and κ 1 = l 2 c , where l c refers to the thickness of the coupling zone Ω c . Before closing this subsection, let us here recall two mathematical results concerning the Arlequin problem, defined above:

• Stability: the mixed Arlequin problem, defined by ( 19)-( 28), has been analyzed mathematically and proved to be well-posed in [START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF] and previously in [START_REF] Ben Dhia | Mathematical analysis of the mixed Arlequin method[END_REF], but under more stringent conditions on the α i (not allowing for regular α i ). • Relevance of the Arlequin formulation: if one modifies, for Ω 0 , the part of the crack contained in Ω f , without modifying the same part of the crack for the local domain Ω 1 , then, under some appropriate hypotheses on the data [START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF]; mainly when, the restriction to Ω f of α 1 is a constant that goes to 1, the solution of the problem ( 19)-( 28) converges to the solution of the (mono-)model problem, defined in Section 2. This result, proved in [START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF] and assessed by previously published numerical tests (e.g. [START_REF] Ben Dhia | Application of the Arlequin method to some structures with defects[END_REF][START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF]), is used in the sequel for a flexible construction of the discrete Arlequin problems.

Discretisation of the Arlequin problem using XFEM

To solve efficiently the Arlequin problem defined in the previous subsection, we suggest the following scheme: let T 0 h be a mesh of the unfractured domain Ω Γ 0 . Let T 1 h be a mesh of Ω 1 , significantly finer than T 0 h , especially in the very near crack tip zone. Let Ω s o denote the set of elements of T 0 h containing the crack tip o. Let Γ 0 denote the part of the crack Γ exterior to Ω s o . Finally, let B 0 = (ϕ 0i , i = 1, n 0 ) be a classical vector-valued finite element basis associated to the mesh T 0 h and denote by V r 0h the finite element space spanned by B 0 . Referring to XFEM recalled in Section 3, the coarse space V r 0h is enriched as following: let K d 0 be the set of indices k such that

• ω k ∩Γ 0 = ∅ (recall that
• ω k stands for the interior of the support of the basis function ϕ 0k ). Let Ω d denote the union of these supports and define Ω + d as in [START_REF] Ben Dhia | Application of the Arlequin method to some structures with defects[END_REF]. The space V r 0h is enriched by the space V d 0h = span {Hϕ 0j , j ∈ K d 0 } where H is the sign field, defined by [START_REF] Dhia | Global-local approaches: the Arlequin framework[END_REF]. The following finite dimensional space:

V 0h = V r 0h + V d 0h = span{(ϕ 0i ) i=1,n 0 , (Hϕ 0j ) j∈K d 0 } (29)
is then taken as an approximation of V 0 . Now, based on the fine mesh T 1 h of the cracked subdomain Ω 1 , a classical finite element subspace of V 1 , denoted V 1h , is constructed. The discrete mediator space, denoted M h , is taken as the restriction of V 1h to the coupling zone Ω c (though other choices are possible, as reported in [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insigths by theoretical investigations of the multiscale Arlequin method[END_REF]). Finally, when using low order finite elements to generate the coarse spaces V 0h , we assume that:

Ω s o ⊂ Ω f (30)
The condition ( 30) is a kind of relaxation of the condition [START_REF] Ben Dhia | Coupling of atomistic and continuum models in the Arlequin framework[END_REF]. It takes into account the finite scale introduced by the coarse discretization and means that the free zone Ω f has at least to contain the coarse elements that contain the singular crack tip. Its optimality in the discrete finite element framework as a low bound for the size of Ω f is assessed in the numerical Section (see also [START_REF] Ben Dhia | Multi-modeling of multi-alterated structures in the Arlequin framework. Solution with a domain decomposition solver[END_REF] for other investigations on this topic). With these finite dimensional spaces, our discrete enriched finite element approximation of the problem defined by ( 19)-( 28), using XFEM, reads as follows.

Find(u r 0h , u d 0h , u 1h , φ h ) ∈ V r 0h × V d 0h × V 1h × M h ; ∀(v r 0h , v d 0h , v 1h , ψ h ) ∈ V r 0h × V d 0h × V 1h × M h a 0 (u r 0h , v r 0h ) + a d 0 (u d 0h , v r 0h ) + c(φ h , v r 0h ) = l 0 (v r 0h ) (31) a d 0 (u r 0h , v d 0h ) + a 0 (u d 0h , v d 0h ) + c d (φ h , v d 0h ) = l 0 (v d 0h ) (32) a 1 (u 1h , v 1h ) -c(φ h , v 1h ) = 0 (33) c(ψ h , u r 0h -u 1h ) + c d (ψ h , u d 0h ) = 0 (34)
where:

a d 0 (u r 0h , v d 0h ) = Ω 0 ∩Ω + d α 0 σ(u r 0h ) : ε(v d 0h ) dx + Ω 0 ∩Ω - d α 0 σ(u r 0h ) : ε(v d 0h ) dx (35) c d (ψ h , u d 0h ) = Ωc∩Ω + d κ 0 ψ h • u d 0h ) + κ 1 (ψ h ) : (u d 0h )) dx (36)
+

Ωc∩Ω - d κ 0 ψ h • u d 0h ) + κ 1 (ψ h ) : (u d 0h )) dx
As a matter of fact, one can observe that, unless the crack tip o belongs to the boundary ∂Ω d , the space V 0h is not a subspace of V 0 . Indeed, assume that the crack tip o / ∈ ∂Ω d (see figure 9). Under this hypothesis, the set Ω s o is necessarily composed of only one element of T 0 h . Consider now the subset of basis functions Hϕ 0k such that Ω s o is included in the support ω k of ϕ 0k . These functions are clearly discontinuous, not only on the part of the crack contained in Ω s o , but also on its fictive straight extension along the tangent vector τ o from the crack tip o till the point o where the extension crosses the boundary of Ω s o (see figure 9). Thus, the discrete problem, defined by ( 31)-( 33), is a nonconforming approximation of the continuous problem ( 19)- [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]. However, and this is one of the main features of the Arlequin method recalled at the end of the previous subsection, this non-conformity is with no significant impact on the solution when the Arlequin parameters are chosen appropriately: in the considered scenario, with a sufficienly large free zone, one has simply to take the weight parameter α 1 constant and approximately equal to one in Ω f to stress the patch model and compute a good Arlequin solution.

Some computational aspects

The discrete problem (31)-( 34) is equivalent to the following linear mixed system:

          K rr 0 K rd 0 0 C r 0 K rd 0 K dd 0 0 C d 0 0 0 K 1 C 1 C r 0 C d 0 C 1 0                     U r 0 U d 0 U 1 Φ           =           F r 0 F d 0 F 1 0           (37)
The stiffness matrix K rd 0 involves terms that require specific numerical integration, when compared to the ones developed for the standard Arlequin formulation, given in [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF] and for the XFEM, given in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. Indeed, the generic terms of K rd 0 read:

(K rd 0 ) ij = a d 0 (ϕ 0i , ϕ d 0j ) ( 38
)
where ϕ 0i is an element of the basis of V r 0h and where ϕ d 0j is an element of the basis of V d 0h . These terms are classically evaluated by summing elements contributions. However, one may have particular situations where coarse elements are crossed by the coupling zone boundary and the crack. For each of these particular elements, one has to create an integration mesh that fits both the parts of the coupling zone boundary and the crack, contained in this element, as shown in figure 10, in order to use classical quatrature formulae in this integration mesh. Before giving some numerical results, let us observe that system 37 is similar to systems derived by domain decomposition methods. It can then be solved by Krylov solvers. Very recently, the FETI method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] has been adapted to the solution of such systems [START_REF] Ben Dhia | Multi-modeling of multi-alterated structures in the Arlequin framework. Solution with a domain decomposition solver[END_REF] and proved to be very efficient. This efficiency is due in particular to the fact that thanks to the energy coupling operator involved by the Arlequin formulation, it is possible de derive "natural" local pre-conditioners (see [START_REF] Ben Dhia | Multi-modeling of multi-alterated structures in the Arlequin framework. Solution with a domain decomposition solver[END_REF], for details).

Numerical results

Three tests are considered in this Section to assess the methodology developed in this paper. For the three of them, we have taken α 0 = 0 in the coupling zone and α 0 = 10 -3 , in the free zone. The first involves a 2D-structure with a kinked crack; the kink being very close to the crak tip. Our results are compared to the ones obtained with a classical globally refined finite element model. The propagation of a crack in

K i K i Γ Ω f Ω c Γ ∂Ω c K i
Fig. 10. Mesh subdivision for the integration issue an elastic beam is studied as a second test and our results are compared to those obtained in [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF]. The third result shows the capability of the method to propagate a crack from a side to the opposite one of a square plate. Let us observe here that for all of the tests, there is no compatibility requirement between the coarse and fine meshes. We consider a homogeneous and isotropic cracked elastic plate in plane strain conditions. The dimensions of the plate (in meters) are given in figure 11. The Young's modulus is equal to 2 10 5 MP a and the coefficient of Poisson is equal to 0.3. The plate is clamped on its lower edge. A given normal displacement is prescribed on the upper edge. The two other edges are free. A straight crack is deviated sharply with a relatively small increment as shown by figure 11. This situation may represent the deviations of a crack during its propagation with small and unknown increments. For the numerical simulation of this test, the two interacting singularities have to be considered, either by introducing the appropiate singular modes by XFEM, for instance, or, as done here, by refining sufficiently the local domain mesh around the two singular points. For comparison purpose, we have first computed a classical globally refined monomodel finite element solution. The mesh size in the near crack tip is taken equal to about 0.02m in the very near crack tip zone. The deformed mesh is shown in figure 12. Three cases are computed with our methodology. In the first, a uniform quadrangular 9 × 9 elements fixed mesh is used for the global domain. The edges of this mesh do not match the crack. The patch is a disk of radius r = 1.2m, centered at the crack tip. The coarse element, of diameter H c = 1.73m, containing the two hot points (the kink and the crack tip) is included in the free zone. The thickness of the coupling zone is taken equal to r/8 and the mesh size of the patch in the very near crack tip zone is of about 0.02m. The deformed meshes are given in the left part of figure 13. A zoom of the critical zone is shown in the right part of this figure. The red lines represent the position of the deformed faces of the crack in the global model. One can observe that while the coarse deformed crack does fit the fine crack in the coupling zone, the two surfaces do not match in the free zone. But, as mentioned at the end of Subsection 4.2, this is with no impact on the Arlequin solution, dominated by the fine local solution in this zone thanks to the fact that the fine model is stressed in this zone by choosing its weight parameter α 1 nearly equal to 1. Figure 14 shows the Von Mises stresses at the Gauss quadrature points of the patch and one can notice that the two singularities are captured (up to the precision of the used numerical models: no h-adaptivity is herein carried out). The similarity of the reference solution given in figure 12 with the one shown in figure 13 is clear. For further comparison, stress intensity factors KI and KII and energy release rates G are computed for the considered reference monomodel and for Arlequin-XFEM models with patches having radii varying from 0.2m to 3m, with 0.125 ratio between the thickness of the coupling zone and the radius of the patches. The right part of figure 15 shows a reasonable decrease, followed by a stabilisation of the error with the increase of the ratio between the the radius of the free zone, denoted by R f , and the coarse mesh size, denoted by H c . Notice that the relative error is already less than 5% when R f Hc is greater than 0.8 and of order 1% when this ratio is of order 1. The left part of figure 15 × 100, at the two critical points o and k, shown in figure 16. One can also notice significant errors when the radius of the free zone is too small with respect to H c . In figure 17, we show the deformed meshes for such an inappropriate choice of the free zone dimension. These results are to be linked to hypothesis [START_REF] Nayroles | Generalizing the finite element method: diffuse approximation and diffuse elements[END_REF].

In a second case, the same procedure has been used by changing only the coarse mesh size: a 19 × 19 coarse structured mesh is used instead of the 9 × 9 mesh. In this case, the diameter of the coarse elements is H c = 0.74m. The obtained results are shown in figure 18. A remarkable feature is the similarity of the profile of error between the two considered numerical models. For instance, for both models, the error in terms of stress intensity factors and energy release rate, is of order 1% when the ratio

R f
Hc is of order 1. e K2 e G Fig. 18. Evolution of the error on chosen mechanical quantities of interest with the increase of the radius of the free zone and with a fixed ratio of patch radius vs. free zone radius

In a third case, the uniform quadrangular 19 × 19 elements fixed mesh is used for the global domain, as in the second case. But in this third case, the patch is a disk of a fixed radius r = 1.5m, centered at the crack tip. It is meshed finely as in the first test. We have considered various thiknesses of the coupling zone which induce various radii R f for the free zones, while keeping fixed the patch mesh. The same comparisons in terms of the same quantities of interest as for the first and second tests are reported in figure 19. Hc . This can be explained by the fact that in the considered case, the fine local mesh is kept fixed, while it is not in the first and second ones.

The double-cantilever beam test

The test studied in [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF] and shown in the introduction of the present paper is analyzed here by our new methodology. It consists in a rectangular L × h 2D homegeneous and isotropic elastic structure in plane stress conditions, with L = 0.3m and h = 0.1m. The Young modulus is equal to 2 10 5 MP a and the Poisson coefficient is equal to 0.3. The structure contains an initial straight crack of lenght a = 0.1m, aligned with the x axis. A small increment of 0.002m with an angle θ of 5.71 o is added to the initial crack and the structure is submitted to a symmetric pointwise opening loading P = 877N, as shown in figure 20. We study the quasi-static growth of the crack. The used procedure is the following: the stress intensity factors using the domain form of the interaction integrals are computed. With these quantities, the propagation's angle, corresponding to the maximum orthoradial stress σ θθ is evaluated and the crack is advanced 0.002m in this direction. Figure 21 shows the propagation of the crack which is qualitatively similar to the one obtained in [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF]. Figure 22 shows zooms of the crack tip zone at the last step of the propagation. The deformed meshes are given by the left part of the figure and the Von Mises stresses in the quadrature points of the patch's mesh by the right part. the approximation of the coarse Arlequin solution. When completely included in the global domain, the radius of the considered (whole) patch is equal to 1.5m and the thickness of the coupling zone is equal to 0.2m. The mesh size of the patch, in the near crack tip zone, is taken equal to 0.01m. The patch geometry is adapted whenever it touches the boundary of the plate. This is simply realized by considering the restriction of the whole patch to the global domain. Observe then that the coupling zone is no more the crown satifying the condition ( 16), but the portion of the crown, satisfying the following adapted condition:

∂Ω 1 ∩ Ω 0 ⊂ ∂Ω c (39)
The same procedure as the one used in the previous subsection is used for the quasi-static propagation of the crack. Four snapshots of the deformed meshes during the propagation of the crack are represented in figure 24. A zoom of the crack tip zone in the first opening step is given in figure 25. This test could possibly be delicate to achieve with classical enrichment methods.

Conclusion

A new methodology using the Heaviside enrichment key feature of the extendend finite element method within the Arlequin framework has been developed When compared to our previous works (e.g. [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF]), the new methodology reduces significantly the global computational costs since a small cracked patch is used in the working zone; the part of the crack located outside of the critical zones being taken into account by means of a Heaviside function.

For the propagation of a crack, the patch, finely meshed, is attached to the crack tip and moves with it over the underlying global domain, approximated by a fixed coarse finite element enriched by a Heaviside function. This methodology is assessed in this paper by some simple but relevant mecanical tests. A kinked crack is calculated by following two different scenarii. For both of these scenarii, reasonable low errors, in terms of quantities of interst, are reached when the size of the free zone is of the same order than the size of the coarse elements containing the singularities. The potential of the methodology to propagate cracks with significantly enhanced flexibility is shown by firstly investigating a double cantilever test without using a priori knowledge of the local singular modes. Our results are shown to compare qualitatively well with the ones given by Belytschko and Black [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF]. Secondly, a final feasibility test is carried out to show the capapility of the method to propagate a crack from a side to another of a structure by adaptig the patch geometry. The extension of this methodology to the simulation of the propagation of

Fig. 1 .

 1 Fig. 1. Propagation of damage and creation of a crack within the Arlequin framework (three selected snapshots)

Fig. 2 .

 2 Fig. 2. Propagation of a crack for the double-cantilever elastic beam (three selected snapshots)

Fig. 9 .

 9 Fig. 9. Element containing the crack tip

Fig. 11 .

 11 Fig. 11. kinked crack problem

Fig. 12 .

 12 Fig. 12. Classical globally refined monomodel finite element solution

Fig. 13 .Fig. 15 .Fig. 16 .

 131516 Fig. 13. deformed meshes of the substrate and the patch

1 =

 1 shows the same tendancy for local relative displacement errors e

Fig. 17 .

 17 Fig. 17. Deformed meshes for inappropriate choice of the free zone dimension

Fig. 19 .

 19 Fig.[START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. Evolution of the error on chosen mechanical quantities of interest with the increase of the radius of the free zone and fixed patch radius

Fig. 20 .

 20 Fig. 20. Geometry of the double cantilever beam test

Fig. 21 .Fig. 22 .Fig. 23 .

 212223 Fig. 21. Deformed meshes of the substrate and the patch

Fig. 24 .

 24 Fig. 24. Deformed meshes of the substrate and the patch

Fig. 25 .

 25 Fig.[START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]. zoom on the initial crack at the beginning of the propagation cracks, while taking into account more complex physics in the working zones, is an ongoing work.