Christian Wolf

Document Ink bleed-through removal with two hidden Markov random fields and a single observation field

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

We present a new method for blind document bleed through removal based on separate Markov Random Field (MRF) regularization for the recto and for the verso side, where separate priors are derived from the full graph. The segmentation algorithm is based on Bayesian Maximum a Posteriori (MAP) estimation. The advantages of this separate approach are the adaptation of the prior to the contents creation process (e.g. superimposing two hand written pages), and the improvement of the estimation of the recto pixels through an estimation of the verso pixels covered by recto pixels; Moreover, the formulation as a binary labeling problem with two hidden labels per pixels naturally leads to an efficient optimization method based on the minimum cut/maximum flow in a graph. The proposed method is evaluated on scanned document images from the 18 th century, showing an improvement of character recognition results compared to other restoration methods.

I. INTRODUCTION

General image restoration methods which do not deal with document image analysis have mostly been designed to cope with sensor noise, quantization noise and optical degradations as blur, defocussing etc. (see [START_REF] Sezan | Survey of recent developments in digital image restoration[END_REF] for a survey). Document images, however, are often additionally subject to further and stronger degradations:

1) non stationary noise due to illumination changes.

2) curvature of the document.

3) ink and coffee stains and holes in the document. 4) ink bleed through : the appearance of the verso side text or graphics on the scanned image of the Christian Wolf is with Université de Lyon, CNRS; and with INSA-Lyon, LIRIS, UMR5205, F-69621, France. email: christian.wolf@liris.cnrs.fr recto side. This is an important problem when very old historical documents are processed. 5) low print contrast. 6) errors in the alignment of multiple printing or imaging stages.

In this paper we concentrate on ink bleed through removal, i.e. the separation of a single scanned document image into a recto side and a verso side. We assume that a scan of the verso side is not available (blind separation). In this case, the task is equivalent to a segmentation problem: classify each pixel as either recto, verso, background, or eventually recto-and-verso (simultaneously), making immediately available the vast collection of widely known segmentation techniques. However, document images are a specific type of images with their own properties and their own specific problems.

At first thought it might be a good idea to interpret the task as a blind source separation problem similar to the "cocktail party" problems successfully dealt with by the (audio) signal processing community. The widely used technique independent components analysis (ICA) has been applied to document bleed through removal, mainly by Tonazzini et al. [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF]. However, ICA assumes a linear model which makes this formulation questionable: ds = Af s where ds is the observation vector, f s is the source vector and A is the mixing matrix. The source vectors, corresponding to the pixels at sites s, are mostly chosen to be three dimensional: the recto signal, the verso signal and an additional signal adding the background color [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF]. In this case, the column vectors of the mixing matrix become the color vectors for, respectively, recto pixels, verso pixels and background pixels, as can be seen setting f s = [ 1 0 0 ] T , [ 0 1 0 ] T and [ 0 0 1 ] T and ds to the respective color vector and solving for A. We can easily verify that the linear hypothesis cannot be justified for ink bleed through by calculating the color of an observed pixel created by a source pixel which contains overlapping recto and verso pixels (f s = [ 1 1 0 ] T ), thus the sum of the color vectors for the recto and the verso pixel, which is unlikely.

The same authors present a non-blind technique also applicable to grayscale images [START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF], the different components corresponding to the recto scan and the verso scan. The inverse of the mixing matrix is calculated using orthogonalization justified by several assumptions on the degradation process. In [START_REF] Tonazzini | A markov model for blind image separation by a mean-field em algorithm[END_REF] the same authors introduce a double MRF model similar to our proposition combined with a likelihood term consisting of a linear mixing model. However, whereas our graphical model is directly employed for classification, the MRF in [START_REF] Tonazzini | A markov model for blind image separation by a mean-field em algorithm[END_REF] guides an EM algorithm estimating the inverse of the mixing matrix. As with the other algorithms based on mixing, the biggest weakness is the linearity of the model. In [START_REF] Tonazzini | Bayesian MRF-based blind source separation of convolutive mixtures of images[END_REF], the model is extended to convolutive mixtures. Sharma presents a non-blind restoration algorithm, i.e. a method which requires a scan of the recto as well as the verso side of the document [START_REF] Sharma | Show-through cancellation in scans of duplex printed documents[END_REF]. The two images are aligned using image registration techniques. A reflectance model taking into account a bleedthrough spread function is created, approximated and corrected with an adaptive linear filter. Another nonblind method is proposed by Dubois and Pathak [START_REF] Dubois | Reduction of bleed-through in scanned manuscript documents[END_REF], where the emphasis is set to the image registration part, the restoration itself is performed using a thresholdinglike heuristic.

Tan et al. propose a non-blind method where the alignment is done manually [START_REF] Tan | Restoration of archival documents using a wavelet technique[END_REF]. Foreground strokes are detected using both images and enhanced using a wavelet decomposition and restoration. The same authors also present a directional blind wavelet method exploiting the hypothesis that handwriting is (very) slanted, and therefore that the strokes of the recto and the verso side are oriented differently [START_REF] Wang | Directional wavelet approach to remove document image interference[END_REF].

Nishida and Suzuki describe a method based on the assumption that high frequency components of the verso side are cut off in the bleeding process [START_REF] Nishida | Correcting show-through effects on document images by multiscale analysis[END_REF]. Their restoration process uses a multi-scale analysis and edge magnitude thresholding. Leydier et al. propose a serialized (i.e. adaptive) version of the k-means algorithm [START_REF] Leydier | Serialized Unsupervised Classifier for Adaptative Color Image Segmentation: Application to Digitized Ancient Manuscripts[END_REF]. Drira et al. propose an iterative algorithm which recursively applies the k-means algorithm to the image reduced with principal components analysis [START_REF] Drira | Restoring ink bleed-through degraded document images using a recursive unsupervised classification technique[END_REF].

The method presented by Don [START_REF] Don | A noise attribute thresholding method for document image binarization[END_REF] is justified by a noise spot model with Gaussian spatial distributions and Gaussian gray value distributions. Under this assumption, thresholds near the class means produce spurious connected components. A histogram of connected component counts is created and thresholded using standard techniques. MRF regularization has already been used for this kind of problem. For instance, Tonazzini et al. present a document recognition system which restores selected difficult parts of the document image with MRF regularization [START_REF] Tonazzini | Analysis and recognition of highly degraded printed characters[END_REF]. As prior model they chose an Ising model with non-zero clique types of 1, 2, 3, and 9 pixels. The observation model contains a convolution term with an unknown parameter. Donaldson and Myers apply MRF regularization with two priors to the problem of estimating a super-resolution image from several low-resolution observations [START_REF] Donaldson | Bayesian super-resolution of text in video with a text-specific bimodal prior[END_REF]. A first prior measures smoothness, whereas a second prior measures a bi-modality criterion of the histogram.

In this approach we ignore degradations no. 1 (non stationarity) and 3 (stains) mentioned at the beginning of the paper and propose an approach based on a stationary model. Non homogeneous observation models will be developed in forthcoming publications. The geometry changes in curvature degradations can be treated with preprocessing, e.g. with the method developed in our team [START_REF] Lebourgeois | Document images analysis solutions for digital libraries[END_REF] or other recent work [START_REF] Brown | Image restoration of arbitrarily warped documents[END_REF], [START_REF] Zhang | An improved physicallybased method for geometrical restoration of distorted document images[END_REF]. The illumination changes inherent in strong curvature degradation can be tackled by non homogeneous observation models.

We formulate our method as Bayesian MAP (maximum a posteriori) estimation problem in terms of two different models:

• the a priori knowledge on the segmented document is captured by the prior model. In our case, the prior model consists of two MRFs, one for each side of the document. • the observation model captures the knowledge on the document degradation process. In a previous paper, we described a MRF model for document image segmentation [START_REF] Wolf | Binarization of Low Quality Text using a Markov Random Field Model[END_REF]. The goal, however, was to learn the spatial properties of text in document images in order to improve the binarization performance. In this paper the emphasis is set to regularization. Therefore, a parametric prior model has been chosen.

The contribution of this paper is threefold: 1) Creation of a double MRF model with a single observation field. 2) Formulation of an iterative optimization algorithm based on the minimum cut/maximum flow in a graph. The proposed inference algorithm is an extension of the widely used α-expansion move algorithm [4][19]. 3) A complete restoration process beginning with an algorithm for the initial recognition of recto and verso labels without using any color or gray value information and finishing with a hierarchical algorithm for the calculation of the background gray value replacing the verso pixel. This paper is organized as follows. Section II proposes a dependency graph for the joint probability density of the full set of variables (the hidden recto and verso variables as well as the observed variables) and derives the prior probability. Section III proposes the observation model. Section IV describes the posterior probability and formulates an iterative inference algorithm based on graph cuts. Section V outlines the estimation procedure for the prior parameters and the parameters of the observation model. Section VI illustrates the preand post processing steps of the method and section VII presents the experiments we performed on scanned document images in order to evaluate the performance of the proposed method. Section VIII finally concludes.

II. THE PRIOR MODEL

MRFs capture the spatial distribution of the pixels of an image by assigning a probability (or an energy) to a given configuration, i.e. a given segmentation result. This is normally used to regularize the segmentation process, i.e. to favor certain configurations which are considered more probable. One of the most widely used assumptions is the smoothness criterion -homogeneous areas are considered more probable then frequent label changes.

This assumption is normally justified1 considering that very often high frequencies in the image content correspond to noise and assuming that the MRF model has been adapted to the prior knowledge on the image content. However, this changes when the observed image is the result of the superposition of two or more "source" images, which is the situation dealt with in this work. In this case, a priori knowledge may be available for each of the source images, but not for the mixture of these images. Applying a simple regularization on the combined image may over-smooth areas which should actually contain high frequency edges due to the superposition process.

We therefore propose to create a prior model with two different label fields : one for the recto side (F 1 ) and one for the verso side (F 2 ). Instead of a segmentation problem with a configuration space of 3 or eventually 4 labels for each site (recto, verso, background, and eventually recto-and-verso), we get 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 background. The advantages of this formulation are three-fold:

• the separation into two different label fields creates a situation where the priors regularize fields which directly correspond to the natural process "creating" the contents (e.g. hand writing letters), as opposed to the single field case, where the prior tries to regularize a field which is the result of overlapping two content fields. • Correctly estimating verso pixels which are shadowed by recto pixels, which is only possible with two separate fields, is not just desirable in the case where the verso field is needed. More so, a correct estimation of the covered verso pixels, through the spatial interactions encoded in the MRF, helps to correctly estimate verso pixels which are not covered by a recto pixel, increasing the performance of the algorithm. • As we will see in section IV, the formulation with separate labels leads to a simple optimization routine based on graph cuts.

Note, that a similar result could be achieved with a single hidden label field and by adapting the prior model such that its regularization handles different label interactions differently. In general this produces rather complicated energy functions equivalent to rather simple interactions in the respective fields. Moreover, the formulation of the inference algorithm would have been more complex.

In the following and as usual, uppercase letters denote random variables or fields of random variables and lower case letters denote realizations of values of random variables or of fields of random values. In particular, P (F =f ) will be abbreviated as P (f ) when it is convenient.

MRFs [15][25] are non causal models on undirected graphs which treat images as stochastic processes. A field F of random variables Fs, s=1 . . . N , where N is the number of variables (pixels in our case), is a MRF if and only if

P (F =f ) > 0 ∀f ∈ Ω and P (Fs=fs|Fr=fr, r = s) = P (Fs=fs|Fr=fr, r ∈ Ns) (1)
where f is a configuration of the random field, Ω is the space of all possible configurations and Ns is the neighborhood of the site s. In other words, the conditional probability for a pixel depends only on the pixels of a pre-defined neighborhood around it.

On a graph, each neighborhood defines a set of cliques, where a clique is fully connected sub graph. According to the Hammersley-Cifford theorem [START_REF] Hammersley | Markov fields on finite graphs and lattices[END_REF] [2], the joint probability density functions of MRFs are equivalent to Gibbs distributions defined on the maxima cliques, i.e. are of the form

P (f ) = 1 Z exp {-U (f )/T } (2)
where Z = P f e -U (f )/T is a normalization constant, T is a temperature factor which can be assumed to be 1 for simplicity, U (f ) = P c∈C Vc(f ) is a user defined energy function, C is the set of all possible cliques of the field and Vc(f ) is the energy potential for the realization f defined on the single clique c.

Given the nature of the problem, the three different label fields (two hidden and one observed) should be considered in a holistic way in order to precisely describe the interactions between the two fields and to define a joint probability distribution on the full set of labels. In the rest of this paper we therefore consider a full 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 F2 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 Ds 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 two nodes r and s if and only if r ∈ Ns and s ∈ Nr.

The model described in this work is generative, i.e. it tries to explain the process of creating the observed variables from the hidden ones. Considering the relationships between the observed variables and the hidden variables, i.e. the degradation process (see section III), we assume a first-order MRF, which means that the following two conditions hold (a common assumption in the MAP-MRF framework):

1) The observations Ds are independent conditional to the hidden labels

F 1 and F 2 . 2) P (Ds|F 1 , F 2 ) = P (Ds|F 1 s , F 2 s )
As a consequence, the dependency graph (see figure 1) contains the following clique types: first order and second order "intra-field" 2 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. For reasons which will become clear in section III, we will set the potentials for the pairwise inter-field cliques to zero, i.e. the second order cliques with one node ∈ F 1 and one node ∈ D as well as the second order cliques with one node in ∈ F 2 and one node ∈ D. The only contributing inter-field cliques are therefore three-node cliques with one node of each respective field (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 n - " U (f 1 ) + U (f 2 ) + U (f 1 , f 2 , d) " /T o (3) Splitting the partition function 1 Z into two factors 1 Z1
and 1 Z2 and uniting the cliques involving hidden labels only in a single function U (f 1 , f 2 ), which is a change of notation only, we get:

P (f 1 , f 2 , d) = 1 Z 1 exp n -U (f 1 , f 2 )/T o • 1 Z 2 exp n -U (f 1 , f 2 , d)/T o (4
) Using the Hammersely-Clifford theorem [START_REF] Baird | Document image defect models and their uses[END_REF] and Bayes rule, we can interpret equation ( 4) as a Bayesian problem, which leads us to:

P (f 1 , f 2 , d) = P (f 1 , f 2 )P (d|f 1 , f 2 ) (5)
The first factor on the right hand side corresponds to the prior knowledge and the second factor corresponds to the data likelihood determined by the observation/degradation model. Inferring the set of hidden labels from the observed labels corresponds to a maximization of the posterior probability (see section IV).

Let us now direct our attention to the prior probability P (f 1 , f 2 ). In the derivation above we saw that the prior was composed of cliques involving hidden variables only, and that there are no cliques containing variables from both fields F 1 and F 2 , which can also directly be seen in the dependency graph: the two hidden label fields F 1 and F 2 do not share any common nodes nor edges. Therefore,

P (f 1 , f 2 ) = 1 Z exp n -U (f 1 , f 2 )/T o = 1 Z exp n - " U (f 1 ) + U (f 2 ) " /T o = 1 Z 1 exp n -U (f 1 )/T o • 1 Z 2 exp n -U (f 2 )/T o = P (f 1 )P (f 2 ) (6)
We can see that the prior probability is actually the product of the two probabilities of the two fields F 1 and F 2 . In other words, the writing on the recto is independent of the writing on the verso page, which makes sense since the two different pages do not necessarily influence each other -they may even have been created by different authors. However, this independence only concerns the situation where no observation has been made. In the presence of observations (the scanned image), the two hidden fields are not independent anymore due to the cliques involving pairs of hidden variables and one observed variable. Intuitively speaking this can be illustrated by the following example: if the observation of a given pixel suggests that at least one of the document sides contains text on this spot (e.g. the gray value is rather low for a white document with dark text), then the knowledge that the recto label is background will increase the probability that the verso pixel will be text.

For a single hidden field, we adopted the widely used Potts model:

U (f ) = X {s}∈C1 αfs + X {s,s }∈C2 β s,s δ fs,f s (7)
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. We chose a stationary and anisotropic model, therefore the single site parameter α depends on the label fs of the corresponding site s whereas the pair site parameters β s,s depend on the direction of the clique (horizontal or vertical). Combining ( 6) and ( 7), we can see that the whole prior energy defined on both hidden fields is given as the sum of two Potts models:

U (f 1 , f 2 ) = X {s}∈C1 α 1 f 1 s + X {s,s }∈C2 β 1 s,s δ f 1 s ,f 1 s + X {s}∈C1 α 2 f 2 s + X {s,s }∈C2 β 2 s,s δ f 2 s ,f 2 s (8)
Note that only the intra-field cliques from the sets C 1 and C 2 are used in the prior model, the clique potentials from the set C 3 are part of the observation model and will be defined in the next section. This choice results in a prior parameter vector θp which consists of 6 parameters (3 for the recto field and 3 for the verso field):

θp = [α 1 , β 1 h , β 1 v α 2 , β 2 h , β 2 v ] T (9)
where superscripts denote the chosen field (1 corresponds to the recto field whereas 2 corresponds to the verso field) and subscripts indicate the direction of the pairwise clique (h denotes horizontal and v denotes vertical).

III. THE OBSERVATION MODEL

We propose a degradation model which allows for several variations and degradations:

• There is no assumption whatsoever on the specific gray level or color of the recto and verso text. However, we assume constant or Gaussian variability of the colors of recto text, verso text and background. • Eventual linear attenuation of the verso colors by the bleed-through process. • 100% opaque ink, i.e. that in the observation field a recto text pixel totally covers the corresponding verso pixel, whereas a recto background pixel does not. • Additional Gaussian noise from the scanning device. The Gaussian assumption may seem to be an oversimplification of the complex process involved in the degradation of historic documents which very often have been stored for centuries in not optimal conditions. However, the choice is motivated by several reasons : the simplicity of the Gaussian function makes the mathematical formulation of the model easy and very often the oversimplifications of the observation model are compensated by the regularizing effect of the prior.

Degradation models designed for document images do exist and are widely used in the document image community. Unfortunately most of them have been developed for the evaluation of document analysis algorithms and therefore have been designed as binary operations, e.g. a series of morphological operations [START_REF] Baird | Document image defect models and their uses[END_REF] [START_REF] Zheng | Morphological degradation models and their use in document image restoration[END_REF]. In [START_REF] Kanungo | Global and local document degradation models[END_REF], Kanungo et al. propose a degradation model which takes into a account the page bending process as well as the perspective distortion and the illumination change which results from it. These formulations are hard, up to impossible to use in a probabilistic estimation framework.

The assumption of 100% opaque ink could theoretically lead to problems in a situation where the verso strokes are darker than the recto strokes, requiring that a lighter recto stroke splits a darker verso strokes into two parts. However, this situation is somewhat unlikely due to the fact that the verso strokes are normally brightened by the ink bleeding process. In our large archive of manuscripts and early printed documents, we did not a find a single occurrence of a situation violating this assumption.

The given assumption has not been chosen to simplify the problem. Partial transparent ink would allow us to use a mixture model and therefore infer information on the value of a verso pixel even if the recto pixel is "text". In this case, the mixture model can sometimes be inverted using ICA or related methods, or it can be included into the likelihood term of our proposed model. This does not question neither the framework itself nor the inference algorithm.

The assumptions described above can be expressed as follows:

D = φ(F 1 , F 2 , µ r , Σ r , µ v , Σ v , µ bg , Σ bg ) + N (0, Σ N ) (10)
where N is the normal law, µ r , µ v , µ bg are, respectively, the mean of the recto, verso and background and the covariances Σ * are defined similarly. Σ N denotes the covariance of the noise process and φ is given as follows:

φ(. . .) = 8 < : N (µ r , Σ r ) if F 1 = 1 A(N (µ v , Σ v )) -b if F 1 = 0 ∧ F 2 = 1 N (µ bg , Σ bg )

else (11)

The matrix A and the vector b correspond to a possible dampening of the colors by the bleed-through process. For completely transparent paper, A corresponds to the identity matrix and b to the null vector.

As can be seen from the form of the combining function φ, the degraded process D consists of three different classes, each following a normal law. The zero mean Gaussian noise of the scanning device as well as the dampening of the verso color are not explicitly modeled, as their effect can be integrated into the parameters of the normal laws of the degraded image. Note that the conditional independence assumptions stated on page 4 are justified by the given model.

As a consequence, the likelihood factorizes as follows:

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

µ s = 8 < : µ r if f 1 s = text µ v if f 1 s = background and f 2 s = text µ bg else Σs = 8 < : Σr if f 1 s = text Σv if f 1 s = background and f 2 s = text Σ bg else (13
) 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.

IV. THE POSTERIOR PROBABILITY AND ITS MAXIMIZATION WITH GRAPH CUTS

Applying Bayes rule to equation ( 5) and combining the result with equation ( 6), we get the posterior probability of the two label fields:

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 ) (14)
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 [START_REF] Brown | Image restoration of arbitrarily warped documents[END_REF], [START_REF] Donaldson | Bayesian super-resolution of text in video with a text-specific bimodal prior[END_REF] and [START_REF] Dubois | Reduction of bleed-through in scanned manuscript documents[END_REF] we get the following energy potential function:

U (f 1 , f 2 , d) = X {s}∈C1 α 1 f 1 s + X {s,s }∈C2 β 1 s,s δ f 1 s ,f 1 s + X {s}∈C1 α 2 f 2 s + X {s,s }∈C2 β 2 s,s δf 2 s , f 2 s + X {s}∈C1 1 2 (ds -µ s ) T Σ -1 s (ds -µ s ) (15)
where and µs and Σs are the sufficient statistics for the observation model given by the labels fs and f s .

To estimate the binary images, equation ( 14) 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][5][10] [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) = X {s}∈C1 h α 1 U 1 (f 1 s ) + α 2 U 1 (f 2 s ) i + X {s,s }∈C2 h β 1 s,s U 2 (f 1 s , f 1 s ) + β 2 s,s U 2 (f 2 s , f 2 s ) i + X {s}∈C1 U 2 (f 1 s , f 2 s ; ds) (16
) where U 1 (fs) = fs, U 2 (fs, f s ) = δ fs,f s and

U 2 (f 1 s , f 2 s ; ds) = 1 2 (ds -µ s ) T Σ -1
s (ds -µ s ). We consider U 2 (., .; .) as 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 ( 16) 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) (17)
It can easily be seen that this is the case of the terms U 2 (., .) in equation ( 16), but not necessarily for all terms U 2 (., .; .). According to the value of the observation ds at site s, U 2 (f 1 s , f 2 s ; ds) may be regular or not. In other words, only 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 iterative formulation for multi label problems, each subproblem is a binary problem where each hidden variable may take two virtual labels: xs and α, where xs 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 labels 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 fs may take 4 Fig. 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) for all s and f 1 for Hs = 1 maximizing

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

U →1 (f 1 , f 2 , d, H)
until Convergence 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 sites s are not regular, and jointly estimating 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

U 2 (f 1 s , f 2
s ; ds) is regular or not:

Hs = 8 < : 1 if U 2 (0, 0, ds) + U 2 (1, 1, ds) ≤ U 2 (0, 1, ds) + U 2 (1, 0, ds)
0 else (18) Figure 2 outlines the inference algorithm, which iteratively calculates the exact solution of two different binary subproblems, maximizing, respectively, 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 graphs, we will rewrite the two functions by reordering some terms. Without loss of generality, 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 variables in 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) and U →1 (f 1 , f 2 ,
U →2 (f 1 , f 2 , d, H) 1 = X {s}∈C1:Hs=0 α 1 U 1 (f 1 s ) 2 + X {s}∈C1:Hs=1 α 1 U 1 (f 1 s ) 3 + X {s}∈C1 α 2 U 1 (f 2 s ) 4 + X {s,s }∈C2:Hs=0∧H s =0 β 1 s,s U 2 (f 1 s , f 1 s ) 5 + X {s,s }∈C2:Hs=1∧H s =1 β 1 s,s U 2 (f 1 s , f 1 s ) 6 + X {s,s }∈C2:Hs =H s β 1 s,s U 2 (f 1 s , f 1 s ) 7 + X {s,s }∈C2 β 2 s,s U 2 (f 2 s , f 2 s ) 8 + X {s}∈C1:Hs=0 U 2 (f 1 s , f 2 s ; ds) 9 + X {s}∈C1:Hs=1 U 2 (f 1 s , f 2 s ; ds) (19)
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 F 2 s as well as one node for each variable F 1 s satisfying Hs = 1. Each unary term is translated into a t-edge, 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 ( 19) 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 ( 16)), but one of the two arguments is fixed in equation ( 19) describing the sub problem. They can therefore be represented as t-edges in the cut graph. Similarly, the terms in line 8 are nonregular 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 n-edges.

Table I 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.

V. PARAMETER ESTIMATION

Since realizations of the label fields F 1 and F 2 are not available, the parameters of the prior model and the observation model must be estimated from the observed data or from intermediate estimations of the label fields. In this work we chose to estimate the parameters in a unsupervised manner, i.e. different parameters are estimated specifically from and for each input image.

To this end, we create initial label fields with the kmeans method and median filter them before applying the supervised estimation technique described below.

Alternatives to this unsupervised approach would be, for instance, iterated conditional estimation [START_REF] Braathen | Global and Local Methods of Unsupervised Bayesian Segmentaiton of Images[END_REF] or the mean field theory [START_REF] Zhang | The mean field theory in em procedures for Markov random fields[END_REF].

The parameters of the observation model are estimated using the classical maximum likelihood estimators, i.e. the empirical means and covariances.

A. The MRF hyper-parameters

For the supervised estimation of the MRF parameters we use least squares estimation, which was first proposed by Derin et al. [START_REF] Derin | Modeling and Segmentation of Noisy and Textured Images Using Gibbs Random Fields[END_REF]. For a single MRF the estimation procedure may be described as follows.

The potential function for a single site s may be given as

U (fs, f Ns , θp) = θ T p N (fs, f Ns ) ( 20
)
where Ns are the intra-field neighbors of s: Ns = {fwe, fea, fno, fso}, θp is the prior parameter vector and N (fs, f Ns ) can be derived from [START_REF] Brown | Image restoration of arbitrarily warped documents[END_REF] as follows:

N (fs, f Ns ) = [ δ fs,1 , δ fs,fwe + δ fs,fea δ fs,fno + δ fs,fso ] T (21)
From [START_REF] Kumar | Discriminative random fields[END_REF] and the basic definition of conditional probabilities on MRFs:

P (fs|Ns) = e -U (fs,f Ns ,θp) P fs∈L e -U (fs,f Ns ,θp) (22)
n-edges for node pairs: Weight Line in eq. ( 19)

F 2 s , F 2 s : (s, s ) ∈ C 2 -β 2 s,s 7
F 1 s , F 1 s : (s, s ) ∈ C 2 ∧ Hs = 1 ∧ H s = 1 -β 1 s,s 5 F 1 s , F 2 s : Hs = 1 U 2 (0, 1, ds) + U 2 (1, 0, ds) -U 2 (0, 0, ds) -U 2 (1, 1, ds)
9 (a) t-edges (to source if weight> 0) for nodes: Weight Line in eq. ( 19)

F 2 s α 2 3 F 2 s : Hs = 0 U 2 (f 1 s , 1, ds) -U 2 (f 1 s , 0, ds) 8 F 1 s : Hs = 1 α 1 2 F 1 s : Hs = 1 X s :H s =0∧f 1 s =1 β 1 s,s - X s :H s =0∧f 1 s =0 β 1 s,s 6
F 1 s : Hs = 1 U 2 (1, 0, ds) -U 2 (0, 0, ds) 9 F 2 s : Hs = 1 U 2 (1, 1, ds) -U 2 (1, 0, ds)

(b) TABLE I THE EDGES ADDED TO THE CUT GRAPH FOR THE PROPOSED INFERENCE ALGORITHM: EACH EDGE CORRESPONDS TO A TERM IN EQ. (19). 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.

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) Fig. 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 example, the potential functions related to the observation model are regular for the middle and for the right pixel (Hs = 1), but not for the left one (Hs = 0).

the following relationship can be derived [START_REF] Derin | Modeling and Segmentation of Noisy and Textured Images Using Gibbs Random Fields[END_REF]:

θ T p [N (f s , f Ns ) -N (fs, f Ns )] = ln " P (fs, f Ns ) P (f s , f Ns ) « ( 23
) where f s is a label different of fs. The RHS of ( 23) can be estimated using histogram techniques [START_REF] Derin | Modeling and Segmentation of Noisy and Textured Images Using Gibbs Random Fields[END_REF], counting the number of occurrences of the clique labellings in the label field. Considering all possible combinations of fs, f s and f Ns , (23) represents an over determined system of linear equations which can be rewritten in matrix form as follows:

N θp = p (24)
where N is a M × 6 matrix, M being the number of data points, i.e. the number of different combinations of label pairs fs and f s having the same neighborhood labels f Ns . The rows of N contain the transposed vectors [N (f s , f Ns ) -N (fs, f Ns )] T . The rows of the vector p contain the corresponding values from the RHS of [START_REF] Lebourgeois | Document images analysis solutions for digital libraries[END_REF]. The system (24) can be solved using standard least squares techniques, as for instance the pseudo inverse.

For practical purposes, note that labeling pairs with one or both of the probabilities P (fs, f Ns ) and P (f s , f Ns ) equal to zero cannot be used. Furthermore, Derin et al. suggest to discard equations with low labeling counts in order to make the estimation more robust.

Adapting the estimation procedure for a double MRF is straight forward. We estimate the parameters on the recto field only, since this field is more stable -all its labels are directly related to the observation field. We get the parameters of the verso field from the assumption that, statistically speaking, the verso field is a flipped version of the recto field, which does not change the parameters.

VI. PRE-AND POSTPROCESSING A. Initialisation of the label fields

An initial estimation of the two label fields f 1 and f 2 is needed for the iterative algorithm described in the previous section. A natural choice is to apply a segmentation technique without regularization, e.g. a k-means segmentation, in order to classify the pixels into three clusters. However, we need to determine for each cluster whether it is background, recto or verso. For most images that could be done using gray level information only, background being the lightest cluster and recto being the darkest one. In order to make this choice more robust, we developed a cluster labeling method which does not use the gray level of the pixels. Instead, it is based on the following two assumptions:

Assumption 1: Most space on the document page is occupied by background.

Assumption 2: The ink is 100% opaque and therefore a recto text pixel completely covers a verso pixel. The first assumption is used to determine the background cluster as the one having most pixels, which is rather straightforward and very efficient. The second assumption is used to determine which one of the two remaining cluster labels -henceforth denoted label a and label b -is the recto label. The basic idea is the following: since recto pixels cover verso pixels, connected components in the (unobservable) verso label field are often cut into several connected components in the observation field when they interact with connected components from the (unobservable) recto label field.

Since we do not have the unobservable label fields -which would make the task trivial -we use histogram statistics on the initial segmentation to exploit this fact. We search for all the places in the image where the two labels interact, i.e. where there are two neighboring pixels, one having label a and the other having label b. The recto label is obtained counting the different connected components touching each label and choosing the label having a minimum number of different components.

As an example, consider the two transitions indicated by the two points A and B in figure 4. For these two transitions, 2 connected components are counted for the gray-ish label, while only one connected component is counted for the darker label.

B. Restoration

The principle of the restoration algorithm is simple: replace the color or gray value of the pixels classified as verso by the color or gray value of the background. Directly using the mean of the background class will produce visible artifacts due to the noise in the image. A better solution is to use the mean of the neighboring pixels classified as background.

Searching these pixels might be laborious in cases where we need to fill larger areas of verso pixels. We therefore resort to a hierarchical pyramidal structure for the calculation of the replacement values. The pyramid is characterized by a 2 × 2 reduction function and a receptive field of 3 × 3 children for each parent site. Each sites holds the mean value of the children's grayvalues as well as the number of children which have been labeled as background. In order to replace a verso pixel, we traverse the pyramid from bottom to the top and stop at a level where enough pixels have been found contributing to a meaningful background value.

VII. 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 nonrestored images is not very high. We chose the open source OCR software "Tesseract" published by Google, mainly because it is easily scriptable 3 , but we also performed some selected experiments with the product of the market leader, Abby Finereader 8 4 , which performs slightly better without changing the ranking of the restoration methods.

As mentioned in section V, the parameters of the method have been estimated in an unsupervised manner, i.e. for each image we estimate specific parameters.

We compared the proposed method with several competing methods. One group of algorithms purely exploits the fact that, according to the hypothesis stated in section III, set recto pixels completely cover verso pixels, without taking into account interactions between neighboring pixels. Examples are the k-means clustering algorithm with k=3 clusters (followed by our restoration algorithm replacing verso pixels, explained in section VI-B), 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 I, statistical source separation is one of the most active areas in bleed-through removal with several works published by Tonazzini et al. on this subject [34][35][36] [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] which we described in section I 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 herself.

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 8 and 9.

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 post-processing as our proposed method.

The k-means method has been tested with different color spaces: grayscale, RGB, L*a*b*, L*u*v*, using the Euclidean distance for each space. Additionally, the CIE94 color metric [START_REF] Melgosa | Testing cielab-based color-difference formulas[END_REF] has been tested for the L*a*b* space (see table III). The best results have been obtained with the L*a*b* and the Euclidean distance.

Figures 5 to 7 illustrate the OCR results on a small image taken from the Gazettes dataset. As we can see, being based on segmentation, the results for kmeans and the two MRF methods are similar. The kmeans result (Figure 5b) is noisy as opposed to the MRF results, the double MRF (Figure 5d) improves the regularity of the single MRF (Figure 5c). The OCR 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 8 and9 show restoration results on two different manuscript images. The source separation methods remove more of the verso text in Figure 8, 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 9 reveils similar strengths and weaknesses, typical to the two types of approaches: the regularized segmentation approaches create crisp images but show localized artifacts, whereas the artifacts created by the source separation methods are more spread out across the image and seem to touch more of the low frequency components.

A. Computational complexity

The computational complexity of the proposed method is dominated by the inference part based on the minimum cut/maximum flow in a graph 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 IV gives effective run times measured on a laptop 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.

VIII. CONCLUSION AND OUTLOOK

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.

Involved in several digitization projects around the world, our team is currently looking into the following perspectives of this work:

• Creation of a homogeneous (adaptive) observation model, which increases the performance on larger images. This model needs to take into account several text colors, as well as other kinds of degradation (see section I).

F 1 F 2 DFig. 1 .

 121 Fig. 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).

F 1 , F 2 ←

 12 Initialize the label fields (e.g. with k-means) H ← Determine the regular sites s repeat -Fix f 1 s for Hs = 0, optimally estimate f 2 s

Fig. 4 .

 4 Fig. 4. (a) an input image; (b) the result of the initial kmeans clustering. As described by assumption 2, the word in the upper part of the image, which belongs to the verso side of the document, has been cut into several connected components by a letter written on the recto side.

  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ërc ¸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.:

Fig. 5 .

 5 Fig. 5. 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 [19] + restoration; (d) double MRF (proposed method).

Fig. 6 .

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

Fig. 7 .

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

Fig. 8 .

 8 Fig. 8. Restoration results on manuscripts. From left to right, top to bottom: input image, k-means, single MRF & α-exp. move [19], double MRF (proposed method), Tonazzini et al. [37] source #1, Tonazzini et al. [37] source #2, Tonazzini et al. [34] source #1, Tonazzini et al. [34] source #2 (source #3 not displayed).

Fig. 9 .

 9 Fig. 9. Restoration results on manuscripts. From left to right, top to bottom: input image, k-means, single MRF & α-exp. move [19], double MRF (proposed method), Tonazzini et al. [37] source #2, Tonazzini et al. [37] source #1, Tonazzini et al. [34] source #3, Tonazzini et al. [34] source #1 (source #2 not displayed).

  ] and

	Color space Distance	Recall	Prec.	Cost
	RGB	Euclidean	78.91	68.23	42,835
	Grayvalue	Euclidean	79.82	68.43	42,675
	L*u*v*	Euclidean	78.30	68.50	41,800
	L*a*b*	Euclidean	78.57	69.43	40,375
	L*a*b*	CIE94 [26]	78.50	68.95	41,142

•

  Creation of a hierarchical Markov model in the lines of[3][18][START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF] which is able to take into account larger neighborhood structures.• Creation of a discriminative model, as for instance a CRF[START_REF] Lafferty | Conditional random fields: Probabilistic models for segmenting and labeling data[END_REF][20] adapted to the nature of the problem, allowing us to model dependencies between the observations. Code in C++ ‡ Code in Matlab/GNU Octave

	Type	Method	1026×1587 2436×3320
	Context-	Niblack [27]	0.6	2.9 †
	free	Sauvola et al. [30]	0.6	2.9 †
		K-Means (k=3)	1.9	10.5 †
	Source-	Tonazzini et al. [34]	36.9	134.6 ‡
	separation Tonazzini et al. [37]	17.0	74.5 ‡
	MRF	Single MRF [19]	7.2	34.9 †
		Double MRF	7.4	36.4 †
	†

TABLE IV EXECUTION

 IV TIMES IN SECONDS FOR VARIOUS METHODS. THE REPLACEMENT OF THE BACKGROUND PIXELS IS INCLUDED IF THE METHOD IS BASED ON SEGMENTATION.

Label changes on the borders of regions can be ignored, dealt with by a separate line processes[START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF] or directly in the main process[START_REF] Charbonnier | Deterministic edge-preserving regularization in computated imaging[END_REF].

The reader may have noticed that we frequently denote the subsets of sites F 1 , F 2 and D as "fields" and will excuse the slight ambiguity with the "full" Markov random field which consists of all three fields.

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

http://finereader.abbyy.com

IX. 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. We thank Madiha Nadri for the help in groundtruthing the document image dataset.

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 6a) produces spurious components, especially in areas with few text, and Sauvola (Figure 6b) tends to cutting characters into several parts due to its assumptions on the grayvalue distribution in the image.

Figures 6c and6d 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 7a-c, delivers similar results: although we can identify a source component which does not include the verso text, the response itself is quite noisy and faint. Postprocessing the image slightly improves the latter but tends to increase the noise. Figure 7d 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 im-provement 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 II 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 independant 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


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

Gérald Gavin
email: gerald.gavin@univ-lyon1.fr

Inference and parameter estimation on hierarchical belief networks for image segmentation

Keywords: Belief networks, image segmentation, graph cuts

We introduce a new causal hierarchical belief network for image segmentation. Contrary to classical tree structured (or pyramidal) models, the factor graph of the network contains cycles. Each level of the hierarchical structure features the same number of sites as the base level and each site on a given level has several neighbors on the parent level. Compared to tree structured models, the (spatial) random process on the base level of the model is stationary which avoids known drawbacks, namely visual artifacts in the segmented image. We propose different parameterizations of the conditional probability distributions governing the transitions between the image levels. A parametric distribution depending on a single parameter allows the design of a fast inference algorithm on graph cuts, whereas for arbitrary distributions, we propose inference with loopy belief propagation. The method is evaluated on scanned documents, showing an improvement of character recognition results compared to other methods.

Introduction

Image segmentation techniques aim at partitioning images into a set of non overlapping and homogeneous regions taking into account prior knowledge on the results as well as a probabilistic model of the observation (degradation) process. Belief networks, but also undirected probabilistic graphical models, are widely used to incorporate spatial dependencies between the pixels into the classification process, often formulating the problem as Bayesian estimation.

In their seminal paper [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF], Geman and Geman introduced a maximum a posteriori (MAP) estimation technique for Markov random fields (MRF). An alternative to the two-dimensional MRFs are hidden Markov chains (MC) on one-dimensional traversals (Hilbert-Peano scans) of an image [START_REF] Abend | Classification of binary random patterns[END_REF] or hybrid MC/MRF techniques [START_REF] Fjortoft | Unsupervised classification of radar images using hidden Markov chains and hidden Markov random fields[END_REF]. The Markov chain models have been extended to belief networks with a pseudo 2D structure [START_REF] Kuo | Keyword spotting in poorly printed documents using pseudo 2-d hidden Markov models[END_REF] and to full 2D connectivity [START_REF] Levin | Dynamic planar warping for optical character recognition[END_REF].

Hierarchical models introduce a scale dependent component into the classification algorithm, which allows the algorithm to better adapt itself to the image characteristics. The nodes of the graph are partitioned into different scales, where lower scale levels correspond to finer versions of the image and higher scale levels correspond to coarser versions of the image. Lower scales manage interactions on pixel level, whereas higher scales manage interactions of groups of pixels (regions). Examples are stacks of flat MRFs [START_REF] Bello | A combined Markov random field and wave-packet transform-based approach for image segmentation[END_REF], pyramidal graph structures [START_REF] Kato | A hierarchical Markov random field model and multitemperature annealing for parallel image classification[END_REF] and the scale causal multi-grid [START_REF] Mignotte | Sonar image segmentation using an unsupervised hierarchical mrf model[END_REF]. Bouman and Shapiro were among the first to propose a hierarchical belief network for image segmentation [START_REF] Bouman | A Multiscale Random Field Model for Bayesian Image Segmentation[END_REF] (refined by Laferte et al. [START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF]). A quad tree models the spatial interactions between the leaf pixel sites through their interactions with neighbors in scale. The main problem of the quad tree structure is the non stationarity it induces into the random process of the leaf sites, which results in "blocky" artifacts in the segmented image.

In the same paper a second model is proposed, where each node has three parents. At first sight, the structure of the dependency graph is similar to our solution (which features four parents for each site), however, the model proposed by Bouman is a pyra-midal model in that the number of nodes decreases at each level. Moreover, as an approximation the inference algorithm proposes a change of the graph after each inference step, whereas in our work the whole graph keeps its full connectivity.

The work described in this paper concentrates on the solution to the lack of shift invariance of the quad tree structured network. Our new model combines several advantages:

-Adaptation to the image characteristics with a hierarchical graph structure (similar to the quad tree) -A stationary random process at the base level (where each site corresponds to one pixel of the input image). -Fast inference using minimum cut/maximum flow algorithms for a subclass of transition probability distributions. The paper is organized as follows: section 2 describes the quad tree structured network and section 3 extends it to the cube. Section 4 presents an inference algorithm using loopy belief propagation and section 5 outlines an interpretation of the hidden variables of the model. Section 6 presents a fast inference algorithm for a parametric class of transition probability distributions. Section 7 describes parameter estimation for the latter class of distributions and section 8 deals with the computational complexity and memory requirements. Section 9 experimentally validates the method. Finally, section 10 concludes.

Quad tree structured models

In the following we describe belief networks, thus graphical models defined on directed acyclic graphs G = {G, E}, where G is a set of nodes (sites) and E is a set of edges. The edges of the graph assign, to each site s, a set of parent sites (written as s -) and a set of children sites (written as s -). The hierarchical nature of the graph partitions the set of nodes into levels G (i) , i ∈ 0..L -1, G (0) being the base level corresponding to finest resolution.

Each site s is assigned a discrete random variable X s taking values from the label set Λ = {0, . . . , C -1} where C is the number of classes. X G , or short X denotes the field of random variables of the graph, whereas X G (l) denotes the field of random variables at level l. The space of all possible configurations of the field X is denoted as Ω = Λ |G| .

In the case of the quad tree structured model [START_REF] Bouman | A Multiscale Random Field Model for Bayesian Image Segmentation[END_REF][START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF], the graph G forms a tree structure with a single root node r ∈ G (L-1) , four children nodes for each node and a single parent node for each node except

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨© " " " " " " " " " " " " " " " " " # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § § § § § § § § § § § § § § § § § § § § §
§ § § § ¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨© " " " " " " " " " " " " " " " " " the root node (see Fig. 1a). Each hidden variable X s is related to an observed variable Y s which is conditionally independent of the other observed variables given a realization of X s : P (y s |x) = P (y s |x s ) and P (y|x) = s∈G P (y s |x s ) (see Fig. 1b).

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $
The objective is to estimate the hidden variables x given the observed variables y, in the case of the MAP estimator this is done maximizing the posterior distribution: x = arg max x∈Ω p(x|y). The absence of cycles in the dependency graph allows the application of an extension of the Viterbi algorithm [START_REF] Viterbi | Error bounds for convolutional codes and an asymptotically optimal decoding algorithm[END_REF][START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF].

The proposed cube model

The main disadvantage of the Markov quad tree is the non stationarity introduced into the random process of the leaf sites G (0) due to the fact that, at any given level, two neighboring sites may share a common parent or not depending on their position on the grid. We therefore propose an extension shown in figures 2a-d (for easier representation the one dimensional case -a dyadic tree -is shown)

First, a second dyadic tree is added to the graph, which adds a common parent to all neighboring leaf sites which did not yet share a common parent. In the full 2D case, three new quad trees are added. The problem is solved for the first level, where the number of parents increased to 4 (for the full 2D model). The result of this step is seen in figure 2b. We repeat the process for each level. New trees connect sites of the original quad tree, but also sites of the trees added at the lower levels. The final result can be seen in figure 2d. Note, that the final graph is not a pyramid anymore, since each level contains the same number of nodes. In general, each node has 4 parents (2 in the 1D representation) except border nodes.

The whole graph can be efficiently implemented by a cube of dimensions N ×M × log 2 max(N, M ) , N × M being the size of the image. In practice, the full height of the cube is not always required, althought a certain number of levels is needed to ensure a proper amount of high-level interaction. The parents and

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ (a) (b) ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § § § § § § § § § § § § § § § § § § § ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § § § § § § § § § § § § § § § § § § § § § § § ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " " # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) (c) (d)
Figure 2: A one dimensional representation of the stepwise extension of the quad tree [START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF] (shown as a dyadic tree) to the proposed cube model children of site s with coordinates x and y on level l are given as follows:

s -=        x + ∆ n , y + ∆ n x + ∆ n , y + ∆ p x + ∆ p , y + ∆ n x + ∆ p , y + ∆ p s -=        x + ∆ n , y + ∆ n x + ∆ n , y + ∆ p x + ∆ p , y + ∆ n x + ∆ p , y + ∆ p where ∆ n = -1 if l = 0 -2 l-1 else ∆ p = 0 if l = 0 2 l-1 else ∆ n = 0 if l = 1 -2 l-2 else ∆ p = 1 if l = 1 2 l-2 else
Sampling from the joint probability distribution represented by the graph can be done in a single top down sweep, since the directed dependency graph does not contain cycles (taking into account the direction of the edges). The graph as it is described in figure 2d (in a 1D representation) corresponds to the hidden part, i.e. the prior model p(x) in the Bayesian sense. As for the quad tree, we consider one observed node related to each hidden node and independence of the observed nodes conditional to the hidden nodes, i.e. the full joint probability distribution factorizes as follows:

p(x, y) = s∈G (0) p(x s ) s∈G (l) , l>0 p(x s |x s -) s∈G p(y s |x s ) (1)
The full Markov cube including observed nodes is parametrized through three probability distributions: the discrete prior distribution of the top level labels p(x), the transition probabilities p(x s |x s -) and the likelihood of the observed nodes given the corresponding hidden nodes p(y s |x s ).

For the inference algorithm, observations at different cube levels are needed. These observations may directly be taken from multi-resolution data, as long as the specific reduction function of the graphical structure is taken into account. In most cases this will require resampling the data in all levels except the finest one. In image segmentation applications, the only available observations are at the base level (y G (0) ). The higher levels can be calculated recursively, e.g. through a mean filter.

Inference with loopy belief propagation

The MAP estimator maximizes equation ( 1), which is hard if done directly. Indeed, the graph of the model contains cycles when the direction of the edges is not taken into account. In other words, the factor graph of the model contains cycles, which makes an exact calculation using message passing algorithms intractable. The junction tree algorithm is capable of calculating the exact solution on general graphs, however, its complexity is exponential in the maximal clique size of the graph after moralization and triangulation. Calculating the optimal triangulation of a general graph being NP-complete. We applied heuristics to find good triangulations of the graph, which did not manage to provide us with tractable clique sizes, i.e. clique sizes which are independent of the size of the graph. Looking at the energy function derived from the (1), i.e. its negative logarithm, Kolmogorov and Zabih showed [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF] that the minimization of the general form of energy functions involving terms of second order or higher is NP-complete, unless the energy potentials are submodular and the number of labels is equal to two1 . Loopy belief propagation (LBP) is an approximative inference technique for general graphs with cycles [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems[END_REF]. In practice, convergence does occur for many types of graph structures. Murphy et al. present an empirical study [START_REF] Murphy | Loopy beliefpropagation for approximate inference: An empirical study[END_REF] which indicates that with LBP the marginals often converge to good approximations of the posterior marginals.

Loopy belief propagation is equivalent to the sumproduct (or max-product) algorithm proposed for

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ ¢ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ £ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¥ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ ¦ § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § § (a) (b) (c) (d)
Figure 3: A cube (left) and its factor graph (right) in one dimension (top) and two dimensions (bottom)

factor graphs [START_REF] Kschischang | Factor graphs and the sum-product algorithm[END_REF]. Any directed or undirected dependency graph can be transformed into a factor graph, and Figure 3 shows the 1D representation of a Markov cube without observed nodes as well as a small part of the full 2D Markov cube with observed nodes and their corresponding factor graphs. The sum-product algorithm operates by passing messages between the nodes of the graph, each message being a function of the corresponding variable node. The message passing schedule for the cube alternates between bottom up and top down passes.

Interpretation of the hidden variables

In this section, we propose a methodology to estimate the conditional probability distributions p(x s |x s -) by taking into account statistical invariance of images belonging to the same corpus. We propose to give an interpretation of the hidden variables x s (i.e. the variables belonging to level l>0) such that :

1. the independence model given by the structure is satisfied. Given a topological enumeration of vertices, a variable x s should be independent of all smaller index variables given its parents x - s . 2. the conditional probabilities are significantly different of conditional probabilities obtained on randomly binary images.

the conditional probabilities are close for all images of the corpus

For simplicity reasons, in the following we describe the binary case (C = 2), the adaptation to multiple labels is straightforward. Let x s be a vertex of the Markov cube and l its level. We call U x the set of vertices of level 0 reachable by a directed path from x. U xs is a 2 l * 2 l square on the image. Then, we naturally define the class of x s as the class with the maximum number of variables U xs (in case of equality, we choose the class randomly). In order to achieve estimation, we just have to compute the frequency of label 0 (resp 1) for each parent configuration. In our corpus, the 3 issues claimed above were verified. This interpretation allows several strategies for estimation of the conditional probabilities:

-Nonparametric definition of the conditional probabilities. Given initial labels at the base level, the labels at the upper levels are computed as described above and the probabilities are estimated using histograming. -Parametric functions are fitted to the histograms. This strategy is pursued in the next section.

Inference with graph cuts

Algorithms calculating the minimum cut/maximum flow in a graph are a powerful tool able to calculate the exact MAP solution on a number of binary labeling problems [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF][START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF] with low order polynomial complexity. It has been shown recently, that energy functions for which the optimal solution is equivalent to the minimum cut in an appropriate graph contain only submodular terms on binary labels [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF],

where submodular means that any projection of a term E(x i , x j , x k , . . . ) onto any subset of two of its arguments satisfies the following condition2 :

E(0, 0) + E(1, 1) ≤ E(0, 1) + E(1, 0) (2)
In the case of the proposed model, not all energy terms are submodular, especially the terms corresponding to the logarithm of the transition probabilities ln p(x s |x s -), so the general model cannot be solved with graph cuts. However, for a large sub class with interesting properties, graph cut solutions can be found. We propose a regularizing term based on the number of parent labels which are equal to the child label: where α l is a parameter depending on the level l, ξ(x s , x s -) is the number of labels in x s -equal to x s and Z is a normalization constant. The such defined transition probabilities favor homogeneous regions, which corresponds to the objective of an image segmentation algorithm. We then decompose it into a sum of binary terms:

p(x s |x s -) = 1 Z α ξ(xs,x s -) l (3)
ln p(x s |x s -) = s ∈s - [(ln α) δ xs,x s ] -Z (4)
where δ a,b is the Kronecker delta defined as 1 if a = b and 0 else. It should be noted that each binary term is submodular since the δ a,b is submodular for all a, b. Fig. 4 shows a cut graph constructed for the dependency graph of Fig. 3b: the cut with minimum cost separating source S from sink T corresponds to the exact MAP estimate for a Markovcube with binary labels (C = 2). Each non terminal node is connected to one of the terminal nodes with weight | ln p(y s |x s = 1)/p(y s |x s = 0)|, according to the sign inside the absolute value. The weights of top level nodes s contain an additional term ln p(x s = 1)/p(x s = 0). Additionally, each non terminal node is connected to each of its parents with an undirected edge and weight ln α.

Minimum cut algorithms are restricted to binary labeling problems (C = 2). Discontinuity preserving energy minimization with multiple labels is NP-hard [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], but the α-expansion move algorithm introduced in [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF] allows to find a local minimum with guaranteed maximum distance to the global minimum. It consists of iteratively applying the minimum cut algorithm to the sub problem of labeling each node of the whole graph between two labels: keeping the current label and changing the a new label α, which is changed at each iteration.

Parameter estimation

We chose the unsupervised technique Iterative Conditional Estimation (ICE) [START_REF] Pieczynski | Convergence of the iterative conditional estimation and application to mixture proportion identification[END_REF] for parameter identi-fication. Given supervised estimators of the parameters from a realization of the full set of variables (X, Y ), an iterative procedure estimates a new set of parameters as the conditional expectation of the parameters conditioned on the observed field Y and the current parameters. In pratice, the expectations are hard to calculate but can be approximated by simulations of the label field based on the current parameters. The initial set of parameters can be obtained from an initial segmentation of the input image.

The prior probabilities of the top level labels β i can be estimated using histogram techniques. Similarly, for most common observation models, maximum likelihood estimators of the sufficient statistics of the conditional distributions are readily available. In this paper, we work with a simple observation model assuming Gaussian noise, requiring as parameters means and (co)-variances for each class. Arbitrary complex likelihood functions are possible using Gaussian mixtures.

For the parameters α l of the transition probabilities, we propose a solution based on least squares estimation similar to the works proposed by Derin et al. for the estimation of Markov random field parameters [START_REF] Derin | Modeling and Segmentation of Noisy and Textured Images Using Gibbs Random Fields[END_REF]. For each level l, we consider pairs of different site labels x s and x s (s ∈ G (l) ) with equal parent labels x s -= x s -. Note that the parent sites are different, whereas their labels are equal. From (3) the following relationship can be derived:

P (x s |x s -) P (x s |x s -) = α ξ(xs,x s -) l α ξ(x s ,x s -) l (5)
Expressing the conditional probabilities through absolute probabilities and taking the logarithm we get:

ln α l [ ξ(x s , x s -) -ξ(x s , x s -) ] = ln P (x s , x s -) P (x s , x s -) (6)
The right hand side of the equation can be estimated from the label process, e.g. by histogramming, whereas the factor in the left hand side can be calculated directly. Considering a set of different label pairs, we can augment this to b T [ln α l ] = a where b is a vector where each element corresponds to the value in the left hand side of equation ( 6) for a given label pair and each value in the vector a corresponds to the right hand side of equation ( 6) for a given label pair. The solution of the over determined linear system can be found using standard least squares techniques. Measuring OCR improvement -To further evaluate our algorithm we tested it on a real application, namely the restoration of images degraded with ink bleedthrough. The goal is to remove the verso component from the recto scan, which makes it a three class segmentation problem. We chose a dataset consisting of 6 images of pages scanned with 600dpi containing low quality printed text from the 18 th century, the Gazettes de Leyde, a journal in French language printed from 1679 to 1798. We tested our method's ability to improve the performance of an optical character recognition (OCR) algorithm and compared it to several widely cited algorithms: k-means clustering, a flat Markov random field (MRF) with graph cuts optimization [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF], as well as two well known methods3 based on source separation [START_REF] Tonazzini | Independent component analysis for document restoration[END_REF][START_REF] Tonazzini | Fast correction of bleed-through distortion in grayscale documents by a blind source separation technique[END_REF].

Figure 6 shows parts of the images together with the OCR results. As can be seen, the cube tends to regularize stronger, which smoothes the region boundaries. We manually created groundtruth and calculated the recall and precision measures on character level, which are given in table 3.

As we can see, our general purpose model outperforms all other segmentation algorithms. The flat MRF directly models the interactions between pixels, which in theory is more powerful than the scale interactions of the markov cube. However, this is only interesting in cases where no long run interactions are needed, e.g. in images with small structures. In images with larger and, more importantly, scale varying content, the hierarchical nature of the markov cube manages to better model the image contents, which directly translates into a better restoration segmentation and restoration performance.

Surprisingly, the two source separation results performed so poorly that we were unable include their recognition performance in the table, since most of the output was blank or gibberish.

Figure 5 shows a zoom into the results comparing the flat MRF and the markov cube and shows that the hierarchical nature of the cube results outperforms the MRF, removes artifacts and fills holes.

Conclusion and discussion

We presented a causal model featuring the advantages of hierarchical models, i.e. scale dependent behavior and the resulting adaptivity to the image characteristics, without the main disadvantage of the quad tree model, i.e. the lack of shift invariance. Bayesian MAP estimation on this model has been tested on binarization and ink bleed through removal tasks and compared to widely used methods. Segmentation quality is better or equal to the results of a MRF model, the difference depending on the scale characteristics of the input image. LBP and a graph cut based algorithm are proposed for inference. Finally, let's note that the model can be extended to 3D data trivially, albeit with an increase in complexity.

Figure 1 :

 1 Figure 1: The quad tree structured model with and without observed nodes.

-Figure 4 :

 4 Figure 4: The cut graph constructed for the binary problem from the dependency graph shown in Fig. 3b, including the two terminal nodes S and T .

-Figure 6 :

 6 Figure 6: Restoration and OCR results on real data, from left to right, top to bottom: input image, k-means, MRF[10], markovcube+GC, 4× Tonazzini et. al [19] (plane #1, plane #2, plane #3, all 3 planes combined), 2× Tonazzini et al. [20] (plane #1, plane #2).

Figure 7 :

 7 Figure 7: Zoom on binarization (C=2) results on synthetic images. From top to bottom, left to right: input image, kmeans+filtering, quad tree, MRF-GC, Cube-LBP, Cube-GC.

& & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & & '

Note that the mode of posterior marginals (MPM) estimator[START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF] is equally difficult to calculate due to the presence of cycles in the dependency graph. Although the partial posterior marginals p(xs|y d (s)), where d(s) denotes the set of all descendants of s, are computable in O(C 8 ), i.e. with high complexity but still in polynomial time, the computation of the posterior marginals p(xs|y) is difficult due to the lack of a common root node which could start a recursion.

In[START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF], the term regular is used instead of submodular

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.

Complexity and storage

Inference complexity for LBP can be given as

where I is the number of iterations and H is the height of the cube and bounded by log 2 max(N, M ) . Storage requires N • M • (H -1) • 15C variables. In practice, LBP in its original form is applicable for low numbers of classes (2, 3 or maximum 4), which is enough for a large number of problems. For higher numbers of classes, the classes may be quantized and the message passing equations slightly changed.

Inference with minimum cut/maximum flow is considerably faster with a complexity 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 1 gives effective run times and memory requirements measured on a computer equipped with a single core Pentium-M processor running at 1.86Ghz.

Experimental results

We evaluated the model on synthetic data as well as real scanned images. In all experiments, we initialized the label field with k-means clustering after low pass filtering.

Pixel level evaluation -To be able to evaluate the model's segmentation performance quantitatively, we applied it to 30 synthetic images of size 512×512 (60 images total) and very low quality subject to multiple degradations: low pass filtering, amplification of ring shaped frequency bands causing ringing artifacts, low quality JPEG artifacts and additional Gaussian noice in various stages (with variances between σ=20 and σ=40). We compared the cube model with different methods of the state of the art: flat MRF segmentation with a Potts model and graph cut optimization [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF], a quad tree [START_REF] Laferte | Discrete Markov image modelling and inference on the quad tree[END_REF] and k-means clustering. The k-means algorithm is only method whose performance is improved when the image is low pass filtered before the segmentation. Table 2 shows the error rates on the different sets, and figure 7 shows a zoom on an example image.