Alexander Serrano

Benjamin Girault

Antonio Ortega

Graph Variogram: A novel tool to measure spatial stationarity

Keywords: Graph Stationarity, Intrinsic Stationarity, Variogram, Empirical Variogram

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Stationarity is a key characteristic used to describe random signals in both signal processing and geostatistical literature. In practice, the most common flavor of stationarity is second order stationarity, which assumes both first and second order moments statistics exist and are invariant under arbitrary spatial displacements or time shifts [START_REF] Hayes | Statistical digital signal processing and modeling[END_REF]. Stationary signals are characterized in the spectral domain by uncorrelated frequencies, so that they can be analyzed and processed using the Power Spectrum Density (PSD) [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF]. These statistical assumptions allow for a wide range of applications where the signal measured is assumed stationary, such as ARMA models [START_REF] Box | Time series analysis: forecasting and control[END_REF] or optimal filtering of noise (Wiener filters) [START_REF] Papoulis | Probability, random variables, and stochastic processes[END_REF].

Geology, atmospheric science, ecology, among other disciplines that require data collection from different spatial locations, often use stationarity [START_REF] Cressie | Statistics for Spatial Data[END_REF]. In particular, spatial stationarity is appealing here as it conveys information about the random signal defined over a Euclidean space.

Our objective is to build a definition of stationarity for signals defined on graphs, that features consistent descriptions of the random signal we observe, irrespective of how we sample continuous space.

Recent contributions to extending stationarity to graph signals have considered both global [START_REF] Girault | Stationary graph signals using an isometric graph translation[END_REF][START_REF] Girault | Signal processing on graphs -Contributions to an emerging field[END_REF][START_REF] Perraudin | Stationary signal processing on graphs[END_REF][START_REF] Marques | Stationary graph processes and spectral estimation[END_REF] and local [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF][START_REF] Girault | Local stationarity of graph signals: insights and experiments[END_REF] definitions of graph stationarity. The former is summarized as second order graph stationarity characterized by uncorrelated spectral components leading to a straightforward definition of graph PSD [START_REF] Girault | Signal processing on graphs -Contributions to an emerging field[END_REF]. In the vertex domain, these signals are invariant through graph translation [START_REF] Girault | Signal processing on graphs -Contributions to an emerging field[END_REF][START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF], hence directly extending the framework of second order stationarity for Euclidean domains using the graph translation instead of a spatial shift. Unfortunately, in this context of sampled Euclidean domains, relating the graph Fourier transform to the continuous Fourier transform is not an easy task, with the closest result being an asymptotic relation between graph Fourier modes and continuous Fourier modes when the number of sampled vertices grows to infinity [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF]. An important consequence of this observation is that the graph translation does not easily relate to a spatial shift, or diffusion operator, in the continuous domain.

To lift the difficulty of extending the shift operator for Euclidean domains to a sensible operator on graph signals, we propose to leave the definition of second order stationarity and work with intrinsic stationarity of random fields defined over Euclidean domains.

Intrinsic stationarity is tightly linked to the variogram, and its empirical equivalent, which allows us to measure spatial correlations for different sized neighborhoods of nodes on both local and global scales. More precisely, instead of extending shifts and PSD, we propose a new measure of spatial stationarity for random sensor network signals defined over an irregular graph domain: The graph variogram.

Our proposed graph variogram extends the original variogram by combining graph structure with distance information between pairs of vertices. More generally, for a given graph structure, local graph variograms can be computed by selecting local subgraphs corresponding to local neighborhoods. Using these subgraphs, we are able to finely measure the local dependence of sample correlation with distance.

Finally, averaging the local dependence of sample correlation with distance over all these local subgraphs yields our graph variogram. A major benefit of this definition is its robustness to variation of the sampling used to measure the continuous signals of interest through the use of distances between samples: the graph variogram provides consistent descriptions of spatial stationary signals, even with nonuniform samples.

The rest of the paper is organized as follows. In Sec. 2, we review and compare definitions of stationarity from signal processing and geostatistics. We then describe the graph signal processing framework and graph stationarity in Sec. 3, and our definition of the graph variogram and its properties in Sec. 4. Finally, Sec. 5 illustrates our contribution on synthetic examples.

BACKGROUND ON STATIONARITY

Second Order Stationarity

For continuous spatial models, we denote the random field defined at any particular spatial location s ∈ R d as X(s). The collection of spatial data from observations of the random field at a subset of N locations D ⊂ R d can be represented as a vector x ∈ R N , where xi = X(s i ), ∀s i ∈ D. For sensor network applications we use either

D ⊂ R 2 (d = 2) or D ⊂ R 3 (d = 3).
In practice, second order or wide sense stationarity is an assumption used for modelling spatial datasets. Second order stationarity describes a random field whose first and second order moments exist and are both invariant under spatial displacements ∀s, h ∈ R d :

E X(s) = µ, (1)
C(s + h, s) = cov X(s + h), X(s) = C(h). (2)
C(•) in ( 2) is called the autocovariance function. A second order stationary random field whose autocovariance function C(h) only depends on the magnitude of the displacement vector, ||h|| = h is called isotropic.

Intrinsic Stationarity

The previous definition of stationarity is more commonly found in the signal processing literature. Instead, researchers in geostatistics focus on the first and second order moment of the difference X(s + h) -X(s), rather than X(•). Past work has shown that intrinsic stationarity is a less restrictive form of stationarity than second order stationarity [START_REF] Myers | To be or not to be... stationary? that is the question[END_REF]. An intrinsic stationary random field [START_REF] Cressie | Statistics for Spatial Data[END_REF] satisfies the following relations ∀s, h ∈ R d :

E X(s + h) -X(s) = 0 (3) 2γ(s + h, s) := var X(s + h) -X(s) = 2γ(h) (4)
where the quantity 2γ(•) is known as the variogram. 2γ(•) measures the variance of the difference between two random field values at two corresponding locations separated by lag displacement h. Moreover, if the variogram of the intrinsically stationary random field only depends on ||h|| = h, it is called isotropic.

An interesting connection can be made between the covariance function and the variogram for a second order stationary random field. Assuming ( 1) and ( 2), the variogram evaluated at spatial locations s + h and s verifies [START_REF] Cressie | Statistics for Spatial Data[END_REF]:

2γ(s+h, s) = E X(s+h) 2 + E X(s) 2 -2E X(s+h)X(s) = 2 C(0) -C(h) , (5)
which in fact proves that any second order stationary random field inherits intrinsic stationarity. The converse, however, is not verified, with Wiener-Levy processes as simple counter-examples [START_REF] Myers | To be or not to be... stationary? that is the question[END_REF].

Variogram & Empirical Variogram

The original variogram for signals defined over continuous space was first introduced in 1963 [START_REF] Matheron | Principles of geostatistics[END_REF]. Letting X(s) be the random field evaluated at a particular location s of a geometrical region V , one calculates the variogram for displacement h as follows:

2γ(h) = 1 |V | V X(s + h) -X(s) 2 ds. (6)
In order to calculate ( 6), all possible pairs of points differing by the exact displacement h would be sampled. Since it is impossible to sample at all pairs of spatial locations, empirical variograms are computed instead using a finite, discrete set of spatial positions in Euclidean space.

In practice, we are also limited to a finite number of realizations of the random field. An empirical measure of the variogram using data obtained at a discrete set of spatial positions {s i }i under the constant mean assumption is the following empirical variance: where

2γ(h) = 1 |N(h)| (s i ,s j )∈N(h) X(s i ) -X(s j ) 2 (7)
N(h) = {(s i , s j ) : s i -s j ≈ h, ∀s i , s j ∈ D}.
When data are irregularly sampled from Euclidean space, the variogram is usually smoothed by using pairs with displacement vectors within some tolerance region described by s i -s j ≈ h [START_REF] Cressie | Statistics for Spatial Data[END_REF].

Empirical Local Variogram

The empirical variogram shown in ( 7) is used to measure spatial variations on a global scale for displacement h. Two additional separate hypotheses can be incorporated into an empirical variogram: isotropy and approximation using local neighborhoods [START_REF] Johannesson | Finding large-scale spatial trends in massive, global, environmental datasets[END_REF]. Under the isotropic assumption, let N(h) denote the set of pairs of spatial locations (s i , s j ) ∈ D such that s i is approximately at distance h from s j . For most applications in geostatistics, the variogram defined in ( 6) is an increasing function of distance h, since for many applications, the farther both samples of the random field are from one another, the more they differ on average [START_REF] Matheron | Principles of geostatistics[END_REF].

More precisely, we can study heterogeneous spatial variations using local neighborhoods. We define Nr(h; Using these two hypotheses, an empirical isotropic local variogram at spatial location s k will have a similar form to that of shown in previous work [START_REF] Johannesson | Finding large-scale spatial trends in massive, global, environmental datasets[END_REF]:

s k ) = {(s i , s j ) | ||s i - s j || ≈ h, ||s i -s k || ≤ r, ||s j -s k || ≤ r, ∀s i , s j ∈ D}.
2γ(h; s k ) = 1 |Nr(h; s k )| (s i ,s j )∈Nr (h;s k ) X(s i ) -X(s j ) 2 (8)
The empirical isotropic local variogram in (8) can be interpreted as a local unweighted average of squared differences between values of X(•) at (s i , s j ) ∈ Nr(h; s k ). A generalization of (8) using weighting function W(s i , s j ; h; s k ) and a corresponding normalization term W(h; s k ) is shown in [START_REF] Johannesson | Finding large-scale spatial trends in massive, global, environmental datasets[END_REF]:

W(s i , s j ; h; s k ) = w h (||s i -s j || -h)ws k (||s i -s k ||) × ws k (||s j -s k ||) (9)
W(h; s k ) = i,j W(s i , s j ; h; s k ) (10)
where w h (•) is a window that describes the tolerance region and ws k (•) defines the neighborhood of s k . This yields the generalized variogram:

2γ(h; s k ) = i,j W(s i , s j ; h; s k ) W(h; s k ) X(s i ) -X(s j ) 2 (11)
The variogram formulation in [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF] is equivalent to that of in (8) when w h (•) and ws k (•) are binary valued (0-1) decreasing functions.

BACKGROUND ON GRAPH SIGNAL PROCESSING

Fundamentals & Notation

Following the GSP literature [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF], we model sensor networks as undirected graphs G(V, E) with a vertex set consisting of N nodes.

The adjacency matrix of the graph is denoted as A ∈ R N ×N . If nodes i and j are connected through an edge weighted by wij, then Aij = Aji = wij. In the context of sensor networks, i ∈ V represents the i th sensor, and xi denotes the value of the random field X(•) at the i th sensor's position s i . For all N sensors, we define our graph signal as x ∈ R N . The Euclidean distance between the i th and j th sensors is denoted as dij = ||s i -s j ||.

There exist multiple ways to construct edge weights wij to model local neighborhood relationships between data points [START_REF] Shuman | Vertexfrequency analysis on graphs[END_REF]. We use here a fully connected graph by linking all nodes with edges weighted by a Gaussian kernel of the Euclidean pairwise distances wij = exp(-d 2 ij /2σ 2 ). Alternatively, one can also reduce the number of edges by connecting the K-nearest neighbors to each vertex, with weights given by a Gaussian kernel. Given a fixed adjacency matrix A, the combinatorial or unnormalized graph Laplacian matrix is defined as L = D -A, where D is a diagonal degree matrix verifying Dii = j wij. Conventional methods in graph signal processing use the eigenvectors of L to define the graph Fourier transform (GFT) [START_REF] Shuman | The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains[END_REF][START_REF] Ortega | Graph signal processing: Overview, challenges, and applications[END_REF].

Graph PSD and Graph Stationarity

It is worthwhile to study how past contributions to graph stationarity using the graph PSD describe intrinsic stationary signals. Fully connected, uniformly and nonuniformly spaced graph structures are generated as shown in Fig. 2a Fig. 2b , respectively. For both sampling schemes, we compute spectral autocorrelation matrices and graph PSDs using 1000 graph signal realizations. Comparing Fig. 2c with Fig. 2d, we observe that even though the signals are spatially stationary, the graph PSD is far from being consistent across realizations.

SPATIAL STATIONARITY USING GSP

Relating the Variogram to the Graph Laplacian

In the context of sensor networks, we have direct access to pairwise distances between nodes, which allows us to introduce the concept of the variogram.

For our problem formulation we use the graph Laplacian quadratic form of a graph signal x:

x T Lx = 1 2 N i,j=1 wij xi -xj 2 = (i,j)∈E wij xi -xj 2 , ( 12
)
where x T Lx measures squared differences of random field values across all edges (i, j) ∈ E. For constant equal weights wij = 1, this quadratic form measures the variation of the graph signal on a global-scale.

Comparing [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF] and ( 12), we see that the empirical isotropic local variogram, 2γ(h; s k ), is a valid graph Laplacian quadratic form. More precisely, the normalized weights

W(s i ,s j ;h;s k ) W(h;s k )
play the same role as that of wij in [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF]. Computing the generalized empirical isotropic local variogram is therefore equivalent to computing the quadratic form [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF] for a graph Laplacian matrix that best describes Nr(h; s k ). Note that for both sampling schemes high correlations for short distances (γ(0) = 0) lead to low energy in the higher spectrum since higher graph frequencies correspond to large variations on edges with large weights.

Defining Nr(h; s k ) on a Graph

Though they are similar, the key difference between ( 11) and ( 12) is that edge weights wij in [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF] correspond to arbitrary distances. This motivates the need to bin distances in order to make [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF] closely match the definition of the empirical isotropic variogram in [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF]. This amounts to creating a different Laplacian for each (h, s k ).

To make (12) match more closely with [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF], we first define a new adjacency matrix A ∆ h , that assigns zero weights to edges with corresponding distance length outside prespecified bin width tolerance range ∆ h := (h -δ 2 , h + δ 2 ). This reduces to applying an elementwise binary operator on the original adjacency matrix. For any pair of nodes (i, j):

[A ∆ h ]i,j = 1, dij ∈ ∆ h 0, dij / ∈ ∆ h (13)
The corresponding graph Laplacian L ∆ h is calculated from A ∆ h . Suppose out of all pairs of nodes, the maximum distance is dmax. ∆ h is therefore contained within (0, dmax). We propose here to break the interval (0, dmax) into H mutually disjoint intervals {∆ h } H h=1 to obtain corresponding Laplacians {L ∆ h } H h=1 . We will study alternative binning methods in a future communication. By this stage, using the isotropic assumption we have built the necessary tools to compute the global variogram shown in [START_REF] Perraudin | Stationary signal processing on graphs[END_REF].

Alternatively, we can compute the generalized empirical local variogram shown in [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF]. To do so, we use non-binary windows centered on vertex k decaying with distance to k. We propose using a graph signal g k ∈ R N localized at node k, to play the role of a local window in the vertex domain. Using g k to modify our previous adjacency matrix in [START_REF] Myers | To be or not to be... stationary? that is the question[END_REF], we write A (∆ h ,k) :

A (∆ h ,k) = G k A ∆ h G k with G k = diag g k (14)
The normalization constant that would be analogous to W (h; s k ) in ( 10) is accounted for by computing the quadratic form of the degree matrix D (∆ h ,k) :

W (∆ h ,k) = 1 T D (∆ h ,k) 1 (15)
Here D (∆ h ,k) accounts for degree variations as a function of both ∆ h and k.

Graph Variogram

We can now proceed with writing our final expression for calculating the GSP-based empirical isotropic local variogram. We call this quadratic form the graph local variogram 2γG, where 2γG(∆ h , k) defines an empirical local variogram relative to center node k ∈ V with respect to distance bin ∆ h :

2γG(∆ h , k) = 2 x T L (∆ h ,k) x 1 T D (∆ h ,k) 1 (16)
where

L (∆ h ,k) = D (∆ h ,k) -A (∆ h ,k) . 2γG(∆ h , k
) is undefined whenever the neighborhood described by (∆ h , s k ) is empty. The 2 in the r.h.s. of ( 16) is included to account for the quadratic Laplacian form in [START_REF] Belkin | Convergence of Laplacian eigenmaps[END_REF] using each edge only once, whereas the original variogram in [START_REF] Girault | Translation on graphs: An isometric shift operator[END_REF] counts each edge twice. We further define the graph global variogram as the average over vertices of the graph local variogram:

2γG(∆ h ) = 1 N N k=1 2γG(∆ h , k) (17)
2γG(∆ h ) measures spatial variations on a global scale for given ∆ h . Finally, we propose the following definition of intrinsic global graph stationarity: Definition: (Intrinsic Global Graph Stationarity) A signal x defined over G(V,E) is intrinsically stationary if and only if

E γG(∆ h , k) = E γG(∆ h ) , ∀k ∈ V (18)
Intuitively, this definition states that no matter the center vertex k, the local graph variogram is the same.

RESULTS

To validate the definition of the global variogram, we now perform experiments showing closeness to the theoretical true variogram. More precisely, we compute realizations of an intrinsic stationary field with an isotropic true variogram 2γ(h). Simultaneously, we generate an independent sampling of the field for each realization of the isotropic model with the goal of showing robustness to both sources of noise. Both uniform and non-uniform samplings (see Fig. 2a and Fig. 2b) have been used within a square grid defined by {∀(x, y) ∈ R 2 : 0 ≤ x ≤ 1, 0 ≤ y ≤ 1} as shown in the following set of figures, which were generated using GraSP [START_REF] Girault | Grasp: A matlab toolbox for graph signal processing[END_REF].

We build the original graphs using two schemes: fully connected and sparse using K-neareast neighbor (K =100). In both cases, the edge weights are chosen using a Gaussian kernel of the Euclidean distance with parameter σ = 0.05.

Each of the four cases allows us to compute covariance matrices Σx for our graph signals defined at N = 500 spatial locations. Using these covariance matrices, we generate 1000 realizations of isotropic intrinsically stationary graph signals as x ∼ N (0, Σx) and compute empirical averages and standard deviations of the empirical global graph variogram. In this section, we are interested in global graph variograms and therefore use g k = 1. We leave the study of local variograms and different g k for future work.

From Fig. 3 intervals ∆ h . Due to the fact that the spatial range of our nodes is limited, |N h; s k | decreases with increasing h. However, comparing and contrasting Fig. 3a with Fig. 3b as well as Fig. 3c with Fig. 3d, we obtain estimates of the true semivariogram model γ(h) with small bias, regardless of the node spatial distribution. Fig. 3e and Fig. 3f both show low variance and bias of global variograms across 100 realizations of graph structures.

CONCLUSIONS

This work extends the concept of the variogram to stochastic signals defined over Euclidean spaces. Experiments involving intrinsically globally stationary theoretical semivariogram models suggest that the global variogram renders consistent statistical measures irrespective of uniform/nonuniform spatial graph structure. Overall, this contribution best served to provide an alternative perspective on graph stationarity for a new class of stochastic graph signals and allows for a definition of shift invariance in the vertex domain that matches that of the Euclidean domain. Future work will include more variogram models, more complex windowing operations, and distance binnings. Finally, we are looking into defining a local stationarity test according to the local graph variogram definition as well as performing variogram analysis using graph filters.

Fig. 1 .

 1 Fig. 1. Illustration of Nr(h; s k ). Dashed lines indicate pairs of sensors that would be excluded in empirical isotropic local variogram calculations. Solid lines indicate pairs of sensors that would be included. Circle defines window's localization with respect to node k.

Fig. 1

 1 illustrates how Nr(h; s k ) would look like for a sensor network modeled using a graph structure.

Fig. 2 .

 2 Fig. 2. Left column uniform sampling. Right column nonuniform sampling. (a)-(b) Intrinsically stationary signal in Euclidean space with γ(h) = 1 -exp -h 0.2 and sampled on the graph vertex locations. (c)-(d) empirical graph PSD as a function of normalized graph frequency, λ.Note that for both sampling schemes high correlations for short distances (γ(0) = 0) lead to low energy in the higher spectrum since higher graph frequencies correspond to large variations on edges with large weights.

Fig. 3 .

 3 Fig. 3. Left column uniform sampling. Right column nonuniform sampling. (a)-(b) Statistics of global variogram using 1 fully connected graph and 1000 signal realizations according to γ(h) = 1 -exp -h 0.2 . (c)-(d) Statistics of global variogram using a single 100-nearest neighbor graph and 1000 signal realizations. There is a shorter support for the global variogram in the uniform sampling scheme since KNN removes most samples across longer distance ranges for uniform structures than it does for the nonuniform structure. (e)-(f) Statistics of global variograms over 100 fully connected graphs each with 1000 signal realizations.

  we observe higher variance in empirical global variogram measures for edge sets ∆j corresponding to longer distance

	1.5					1.5
	1					1
	0.5				γ µγ ∞	0.5					γ µγ ∞
					±σγ ∞						±σγ ∞
	0					0
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
			h						h
			(a)						(b)
	1.5					1.5
	1					1
	0.5				γ µγ 100	0.5					γ µγ 100
					±σγ 100						±σγ 100
	0					0
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
			h						h
			(c)						(d)
	1.5					1.5
	1					1
	0.5			γ µ γ∞(h) G ±σ γ∞(h) G	0.5				γ µ γ∞(h) G ±σ γ∞(h) G
	0					0
	0	0.2	0.4	0.6	0.8	1	0	0.2	0.4	0.6	0.8	1
			h						h
			(e)


Ait Bella

Fatim Zahra

Mohammed El Rhabi

Abdelilah Hakim

Amine Laghrib

Reflectance and illumination decomposition using nonlocal models

Keywords: Nonlocal operators, Partial differential equation, Text image, Viscosity solutions theory

In this paper, we investigate two new reflectance and illumination decomposition models based on a nonlocal partial differential equation (PDE) applied to text images. Taking into consideration the higher regularity level of the illumination compared to the reflectance, we propose a nonlocal PDE which deals with repetitive structures and textures that characterize the text image much better compared to the classical local PDEs. The aim of this approach is to use the repetitive features of the reflectance to efficiently extract it from the nonuniform illumination. This idea is motivated by extending the range of application of the nonlocal operators to such problem. Numerical experiments on both grayscale and color text images show the performance and strength of the proposed nonlocal PDE.

Introduction

In image analysis, reflectance and illumination decomposition is a basic process which aims to enhance the image by reducing or removing the effects of the nonuniform illumination. This technique, which is named the Retinex problem, represents a necessary pretreatment step in a wide range of image processing applications.

The image decomposition is considered as an ill-posed problem since it consists in extracting two unknown parts from an observed image. The Retinex theory was first proposed by Land and McCann in [START_REF] Edwin | Lightness and retinex theory[END_REF] to address and improve the quality of an image when the lighting conditions are not satisfactory [START_REF] Edwin | The retinex theory of color vision[END_REF][START_REF] Mccann | Lessons learned from mondrians applied to real images and color gamuts[END_REF]. The main idea of this theory is to consider the human visual system as a color perception where we can define the reflectance associated to a field in which both illumination and reflectance are unknown. Indeed, for a given scene, our visual system detects and sees the same color even if the illuminations conditions vary. This explains that the color of the objects is considered invariant in spite of illumination transformations.

The main goal of Retinex theory is to decompose a given image u into two components represented by the following model: u(x, y) = I(x, y)R(x, y).

(

) 1
The model assumes that the observed image u is the product of two main components [START_REF] Horn | Robot vision[END_REF], namely the illumination I and the reflectance image R.

According to the single-scale Retinex algorithm, the reflectance is estimated as described by the following equation: log(R i (x, y)) = log(u i (x, y))log(F (x, y) * u i (x, y)),

where * is the convolution operator, R i is the reflectance image on the i th channel and F is a Gaussian filter defined by

F (x, y) = C exp -( x 2 +y 2 σ 2 ) , (3)
where C is a normalization factor. This single-scale Retinex formula is extended to a multi-scale version in [START_REF] Daniel J Jobson | A multiscale retinex for bridging the gap between color images and the human observation of scenes[END_REF][START_REF] Petro | Multiscale retinex. Image Processing On Line[END_REF]. This approach is pursued and introduced by Brelstaff [START_REF] Brelstaff | Computing lightness[END_REF] and Horn [START_REF] Berthold | Determining lightness from an image[END_REF], who proposed a Retinex algorithm based on Poisson equation. In 2003, Kimmel and Elad proposed a variational Retinex approach [START_REF] Kimmel | A variational framework for retinex[END_REF]; to form a penalty functional, the authors call into question some hypotheses about the illumination image: the illumination is spatially smooth, it is close to the intensity image u and I ≥ u, moreover they assume that the illumination continues smoothly as a constant beyond the image boundaries. Other approaches are presented using the same principle. First, Ma and Osher have introduced a total variation regularized formulation and an efficient algorithm based on the split Bregman method, then, they took advantage of the nonlocal total variation to enhance the resulting images [START_REF] Ma | A tv bregman iterative model of retinex theory[END_REF]. Further, Chen et al have used a TV-L1 based variational Retinex approach [START_REF] Ma | A tv bregman iterative model of retinex theory[END_REF] . More recently, Liang and Zhang suggested a reflectance and illumination decomposition model for the Retinex problem using a convex variational model based on high-order total variation and L 1 decomposition [START_REF] Liang | Retinex by higher order total variation lˆ{1} decomposition[END_REF]. This method can effectively decompose the gradient component of images into salient edges and slightly smoother illumination field.

Another famous method to correct the non-uniform illumination is the Homomorphic filtering. This approach was introduced by Oppenheim et Al [START_REF] Van Oppenheim | Nonlinear filtering of multiplied and convolved signals[END_REF], the main idea of this technique is to decompose the image into high and low frequencies components, as a result the enhanced image is obtained through a high-pass filter in order to remove the illumination component. The improvement is then achieved by the simultaneous compression of the intensity range and the contrast enhancement. This method gave reasonable results and have been used in many applications [START_REF] Fries | Image enhancement by stochastic homomorphic filtering[END_REF][START_REF] Fan | Homomorphic filtering based illumination normalization method for face recognition[END_REF][START_REF] Gorgel | A wavelet-based mammographic image denoising and enhancement with homomorphic filtering[END_REF].

On the other hand and precisely in document-image processing, many approaches are proposed to enhance the quality of a document-image. This generates a text which is easier to read in order to extract it thereafter using an optical character recognition (OCR). In 1995, W. Jiang has investigated a thresholding and enhancement method [START_REF] Wei | Thresholding and enhancement of text images for character recognition[END_REF] which consists of converting text images of low spatial resolution to bi-level images of higher spatial resolution. In the same principle, a simple algorithm based on empirical mode decomposition was proposed in [START_REF] Pei | Enhancement of uneven lighting text image using line-based empirical mode decomposition[END_REF] to improve the illumination problem for document-image. However, there are several other methods specified for document-image enhancement with different degrees of success, see for examples [START_REF] Visvanathan | Enhancement of camera captured text images with specular reflection[END_REF][START_REF] Caner | Shape-dna: Effective character restoration and enhancement for arabic text documents[END_REF][START_REF] Wei | Thresholding and enhancement of text images for character recognition[END_REF][START_REF] Fan | Cartoontexture image decomposition via non-convex low-rank texture regularization[END_REF].

The main objective of this paper is to provide an improved nonlocal approach to remove the shadow from a document-image. Nonlocal operators have received much attention in the past decade. Such operators were first introduced by Gilboa and Osher [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. The purpose of their paper is the use of nonlocal operators to treat signal and image processing by adopting the so called nonlocal PDE. Convinced by the effectiveness and efficiency of the nonlocal operators in dealing better with image textures, we introduce two new nonlocal PDEs to enhance the text quality and make it more meaningful and expressive for the human eye. These two nonlocal PDEs take into consideration the fact that illumination varies more smoothly than reflectance and also consider the effect of the illumination on image features such as texture. This can help the model to be able to take into account the surroundings of the object much better compared to some classical approaches.

The paper is organized as follows. In section 2, we will focus on introducing a local model which we aim to extend to the nonlocal case, we recall some essential and basic definitions and properties of viscosity solutions theory and we prove the existence of a unique solution of this problem, section 3 deals with the formulation of the new nonlocal models, in section 4, two algorithms are presented to solve the proposed nonlocal models, and in order to assess the validity of the proposed approach we will present some numerical tests in section 5, finally we shall sum up by a conclusion where we will outline the big points covered in the paper.

Problem Formulation

We consider a document-image u acquired by a camera phone. In general, many degradation factors are involved by this acquisition, such as : noise, non-uniform illumination and motion blur. Note that in this paper, the blur is not considered. In fact, the studied model is based on the separation of the two components of the image: the illumination I and the reflectance R u(x, y) = I(x, y)R(x, y) + n(x, y),

where n is the noise. As discussed above, we firstly neglect the noise and we apply the logarithm function, we obtain then the additive model log(u(x, y)) = log(I(x, y)) + log(R(x, y)).

(

) 5
Our aim is to consider a decomposition approach to recover the illumination I by estimating the term " log I". In particular, we focus on extracting illumination I using a nonlocal PDE in a well posed functional framework. Our proposed model is considered as a nonlocal version of the PDE proposed in [START_REF] Mahani | Text enhancement by pde's based methods[END_REF]. The idea behind this choice is related to recently developed nonlocal operators [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. Therefore, in the following, we recall the studied PDE and also we present some theoretical results. The local PDE is defined by

         ∂u ∂t = max(0, Δu) in Ω × (0, T ) ∂u ∂n = 0 on ∂Ω × (0, T ) u(., 0) = u 0 in Ω (6)
where Ω is a bounded open subset of R 2 and Δ is the Laplacian operator defined by

Δu = ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 . ( 7
)
The method proposed in [START_REF] Mahani | Text enhancement by pde's based methods[END_REF] which aims to resolve the equation ( 6) supposes that the document-image u has been acquired by a cameraphone. The authors consider two problem types; namely: image distortions and noise. Especially, they address the so-called non uniform illumination or variations of brightness. In image processing applications, a shadow is considered as a region with low lightness and high gradient contours. There are in fact two kinds of shadows, the own shadow and the shadows arising from the acquisition process. Own shadow occurs when the light hits a surface with a slope change. Otherwise, the brightness of pixels decreases when the angle of incidence deviates from the normal of the surface. While the brightness reaches its minimum value in the case when the incident light and the surface normal are orthogonal. A shadow always occurs when the light source is obscured by an object before the light reflection on the surface. Taking into account that the image contains only text, the background is preponderant compared to the text in the image (or there is much more background than text in the image). Also knowing that the intensity of the background is superior to the intensity of the text (light background on dark writing), it is easy to approximate the background. Indeed, intuitively, in order to estimate the background, we replace each pixel by considering the maximum mean of its neighboring pixels and we iterate that is to say if u n ij represents the intensity of the image in the spatial location (i, j) at the n th iteration, we get:

u n+1 ij = max � u n ij , � u 1,n ij , � u 2,n ij , � u 3,n ij , � u 4,n ij � , (8)
where � u k is the mean calculated in the k th direction using the neighboring pixels (see Fig. 1). By subtracting u n ij from both sides of the equation ( 8), the discrete derivatives in all directions appear and we get the discretized version of the equation:

∂u ∂t = 1 2 max � 0, ∂ 2 u ∂x 2 , ∂ 2 u ∂y 2 , ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 + 2 ∂ 2 u ∂x∂y , ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 -2 ∂ 2 u ∂x∂y � ( 9
) in a similar way, considering only the horizontal and vertical directions, one can easily get the equation ( 6). Realeyes 3D company and M. El Rhabi have successfully implemented this equation as an embedded application for document scanning in more than 150 million smartphones in the world from 2006 to 2010. In this paper, as a first step, we choose a simplified form of this equation which concides with equation ( 6). Then, we rewrite (9) in its non local formulation which we analyse.

The equation ( 6) has demonstrated its robustness in the extraction of the variable luminance and it has given quite satisfactory results (see [START_REF] Mahani | Text enhancement by pde's based methods[END_REF]). However, the theoretical aspect of this equation has not been investigated, which left a lack of information on the behavior and the existence of the solution. As a first contribution in this work, we are interested in the theoretical side of this model. The complexity to show the existence of a solution to the equation ( 6) comes from the fact that we cannot say anything about the maximum of two operators. We intend here to find a continuous viscosity solution of [START_REF] Buades | Nonlocal image and movie denoising[END_REF]. We first give a brief overview of some essential notions of the viscosity solutions theory [START_REF] Michael G Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF][START_REF] Giga | Surface evolution equations: A level set approach[END_REF][START_REF] Bardi | Viscosity solutions and applications[END_REF]. The theory deals with some partial differential equations of the form :

F (x, u, ∇u, ∇ 2 u) = 0, ( 10
)
where

x is defined in an open subset Ω of R n . F is a real-valued function on D = Ω × R × R n × S n
, where S n stands for the set of symmetric n × n matrices, u denotes the real-valued unknown function on Ω, ∇u and ∇ 2 u designate respectively, the gradient and the Hessian matrix of u.

The function F is called degenerate elliptic if it verifies the following inequality:

F (x, r, p, X) ≤ F (x, r, p, Y ) for Y ≤ X, (11)
where Y ≤ X signifies:

< Xξ, ξ >≤< Y ξ, ξ > ∀ξ ∈ R n . ( 12
)
If F is degenerate elliptic, we can say that it is proper if:

F (x, r, p, X) ≤ F (x, s, p, Y ) for Y ≤ X, r ≤ s. ( 13
)
Let u : E → R ∪ {±∞}, where E is a set of a metric space. The upper semi-continuous envelope u * , and the lower semi-continuous envelope u * of u are given respectively by

u * (y) = lim r→0 sup{u(z)/ z ∈ B(y, r) ∩ E}, u * (y) = lim r→0 inf{u(z)/ z ∈ B(y, r) ∩ E},
where B(y, r) is a closed ball of radius r centred at y ∈ E.

Definition 1. [START_REF] Giga | Surface evolution equations: A level set approach[END_REF] Let Ω be an open set in R n and T > 0 a fixed time parameter.

Let [START_REF] Ghergu | Nonlinear PDEs: Mathematical models in biology, chemistry and population genetics[END_REF] in U , if it is both a viscosity subsolution and a viscosity supersolution of (14) in U .

F : Ω × [0, T ] × R × R n × S n → R be a continuous function. Let U be an open set in Ω×]0, T [. A function u : U → R ∪ {-∞} is a viscosity sub-solution of ∂u ∂t + F (x, t, u, ∇u, ∇ 2 u) = 0, ( 14
)
in U if: • u * (x, t) < ∞ for (x, t) ∈ U . • If (ϕ, x, t) ∈ C 2 (U ) × U satisfies max U (u * -ϕ) = (u * -ϕ)(x, t). ( 15
)
then ∂ϕ ∂t (x, t) + F (x, t, u * (x, t), ∇ϕ(x, t), ∇ 2 ϕ(x, t)) ≤ 0. ( 16
) Definition 2. [15] A function u : U → R∪{+∞} is a viscosity supersolution of (14) in U , if • u * (x, t) > -∞ for (x, t) ∈ U . • If (ϕ, x, t) ∈ C 2 (U ) × U satisfies min U (u * -ϕ) = (u * -ϕ)(x, t), ( 17
)
then ∂ϕ ∂t (x, t) + F (x, t, u * (x, t), ∇ϕ(x, t), ∇ 2 ϕ(x, t)) ≥ 0 (18) Definition 3. [15] A continuous function u : U → R ∪ {±∞} is a viscosity solution of
We apply theorem 3.1 in [START_REF] Ishii | Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain[END_REF] to get the following existence result:

Theorem 1. Let Ω be a regular bounded open subset of R 2 . Let u ∈ U SC([0, T [×Ω) (upper
-max(0, Δu) = 0 in Ω×]0, T [ (19)
and

∂u ∂n = 0 on ∂Ω×]0, T [ (20)
If u(x, 0) ≤ v(x, 0) f or x ∈ Ω (21)

Proof

First of all, we use the following decomposition of the max operator max(0, Δu) = H(Δu)Δu [START_REF] Wei | Thresholding and enhancement of text images for character recognition[END_REF] where H is the Heaviside step function, defined by

∀ x ∈ R, H(x) = � 0 if x < 0 1 if x ≥ 0 (23)
Afterwards, we consider the function F : S n →] -∞, 0], given by:

F (X) = -H(trace(X))trace(X). (24)
According to the theorem 2.1 in [START_REF] Ishii | Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain[END_REF], we have to prove that F is continuous and there exists a continuous function w

: [0, ∞[→ [0, ∞[ satisfying w(0) = 0, such that if X, Y ∈ S n and µ 1 , µ 2 ∈ [0, ∞[ satisfy � X 0 0 Y � ≤ µ 1 � I -I -I I � + µ 2 � I 0 0 I � , (25)
then

F (X) -F (-Y ) ≥ -w(µ 2 ). ( 26
)
We start then by the continuity which is evident, since the functions

f : S n → R X � → trace(X) and g : R → ] -∞, 0] x � → -H(x)x
are continuous. Now, we suppose that there exist µ 1 , µ 2 ∈ [0, ∞[ such that the inequality ( 25) is satisfied.

We multiply both sides of the inequality by the matrix

� I I I I � ,
We then apply the trace operator, we obtain

trace(X) + trace(Y ) ≤ 2nµ 2 . ( 27
)
Afterwards, taking the decreasing function g(x) = -H(x)x yields

-H(2nµ 2 )2nµ 2 = -2nµ 2 ≤ -H(trace(X) + trace(Y ))(trace(X) + trace(Y )) ≤ 1 2 (-trace(X) -trace(Y ) -|trace(X) + trace(Y )|) ≤ 1 2 (-trace(X) -trace(Y ) -|trace(X)| + |trace(Y )|) = F (X) -F (-Y ). ( 28
)
Finally, we take w(r) = 2nr which satisfies w(0) = 0. According to the theorem 2.1 and theorem 3.1 in [START_REF] Ishii | Nonlinear oblique derivative problems for singular degenerate parabolic equations on a general domain[END_REF], the results of comparison and existence of a unique viscosity solution is reached. � As demonstrated in [START_REF] Mahani | Text enhancement by pde's based methods[END_REF], the used PDE (6) has given good results, but does not preserve well the fine details of the text. This can be interpreted by the use of a local image reconstruction algorithm using observations in a neighborhood of a pixel of the interest. However, in some cases the selection of the neigborhood plays a major role in the quality of the image reconstitution. This led us to think to the nonlocal version of the equation ( 6). In the nonlocal techniques, algorithms analyse data in a larger neigborhood and collect the observations from the whole image, searching similar features, which is much better and accurate. Indeed, the nonlocal methods, compared to local methods, are more effective, in particular to restore fine and repetitive structures of the image, which is the case of document-images. These approaches which are called also patch methods, allow to characterize a pixel with a vector of attributes wider than the intensity or color value of the pixel only. A patch representing a pixel is considered as a square sub-image of the original image, centered at this pixel, of a selected width. The patch-based methods were first used in the texture synthesis [START_REF] Efros | Texture synthesis by nonparametric sampling[END_REF]. Also, they were then used for filtering and restoring degraded image [START_REF] Buades | A non-local algorithm for image denoising[END_REF]. More recently, it have also been extended to inpainting problems [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. In the following section, we introduce the proposed nonlocal model to extract the illumination I and we give a result of the existence and uniqueness of the solution.

The proposed nonlocal approaches

In this section, we introduce the two nonlocal PDEs and we give some theoretical results. We start by the definition of the first model.

The first model

As discussed above, we propose a new model for document-image enhancing. Therefore, we intend to obtain a new nonlocal generalization of the local model presented in [START_REF] Mahani | Text enhancement by pde's based methods[END_REF]. As a preliminary, we recall some definitions of nonlocal operators introduced by Gilboa and Osher [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF]. Let Ω ⊂ R 2 , and u : Ω → R denote a real function. The nonlocal gradient ∇ w u(x) at x ∈ Ω is defined as the vector of all partial derivatives ∇ w u(x, .):

(∇ w u)(x, y) := (u(y) -u(x)) � w(x, y), y ∈ Ω, (29)
where w stands for a symmetric and positive weight function.

For any pair of

v 1 , v 2 : Ω × Ω → R, one defines (v 1 .v 2 )(x) := � Ω v 1 (x, y)v 2 (x, y)dy, (30)
and the inner product

< v 1 , v 2 > := � Ω×Ω v 1 (x, y)v 2 (x, y)dxdy. (31)
With the above definition of inner product and nonlocal gradient, the nonlocal divergence of v : Ω × Ω → R at x is defined as the negative adjoint of the nonlocal gradient:

(div w v)(x) := � Ω (v(x, y) -v(y, x)) � w(x, y)dy. (32)
The nonlocal Laplacian of u : Ω → R is then defined by:

(Δ w u)(x) := 1 2 div w ((∇ w u)(x)) = � Ω (u(y) -u(x))w(x, y)dy. (33)
At this point, we can use these operators to introduce the proposed nonlocal PDE. Indeed, to estimate the variable illumination we introduce the following model:

� ∂u ∂t = max(0, Δ w u) in Ω × [0, ∞[ u(., 0) = u 0 in Ω (34)
where u 0 ∈ L 2 (Ω) is a given degraded image. This problem is seen as an evolution equation with homogeneous Neumann boundary condition, which mean that when integrating in Ω, we suppose that there is no flux across the boundary. Before resolving this PDE, we have to check the existence and uniqueness of the solution. We are now ready to state our main result Proof Basically, we use the following decomposition of the max operator max(0,

Δ w u) = 1 2 (Δ w u + | Δ w u |) (35)
According to the definition of Δ w u and (35), we consider the bounded linear operator A given by

Au = 1 2 u(.) � Ω w(., y)dy. (36)
The operator A is the infinitesimal generator of S(t) = e -tA , t ≥ 0. On the other hand, the function:

f (u) = 1 2 ( � Ω u(y)w(., y)dy + | � Ω (u(y) -u)w(., y)dy|) (37)
is an M-Lipschitz continuous function, indeed ∀t > 0, ∀u, v ∈ L 2 (Ω):

| f (u(t)) -f (v(t)) | L 2 (Ω) ≤ | Ω || w | L ∞ (Ω×Ω) (2 | u -v | L 1 (Ω) + | u -v | L 2 (Ω) ) ≤ | Ω || w | L ∞ (Ω×Ω) (2c + 1) | u -v | L 2 (Ω)
We introduce now the operator ψ given by:

ψ(u)(t) = S(t)u 0 + � t 0 S(t -s)f (u(s))ds. ( 38
)
We shall show that ψ maps

E α = {u ∈ C([0, ∞[; L 2 (Ω)); | u | Eα < +∞} equipped with the norm | u | Eα = sup t (e -αt | u(t) | L 2 (Ω) ) into E α ∀α > 0. E α is a complete space since it is a closed subset of C([0, ∞[; L 2 (Ω)).
we have

∀t > 0, | ψ(u)(t) | L 2 (Ω) ≤ | S(t)u 0 | L 2 (Ω) + � t 0 � S(t -s) � L(L 2 (Ω)) | f (u(s)) | L 2 (Ω) ds ≤ | u 0 | L 2 (Ω) + � t 0 | f (u(s)) | L 2 (Ω) ds ≤ | u 0 | L 2 (Ω) + � t 0 (| f (u(s)) -f (0) | L 2 (Ω) + | f (0) | L 2 (Ω) )ds ≤ | u 0 | L 2 (Ω) + � t 0 (M | u(s) | L 2 (Ω) +C)ds. ( 39
)
We multiply this last inequality by e -αt for t > 0. We obtain then :

∀t > 0; e -αt | ψ(u)(t) | L 2 (Ω) ≤ e -αt | u 0 | L 2 (Ω) + � t 0 e -α(t-s) e -αs (M | u(s) | L 2 (Ω) +C)ds ≤ e -αt | u 0 | L 2 (Ω) + M sup s (e -αs | u(s) | L 2 (Ω) ) � t 0 e -α(t-s) ds + tCe -αt ≤ | u 0 | L 2 (Ω) + 1 α M | u | Eα +C sup t (te -αt ). ( 40
)
Thus

sup t>0 e -αt | ψ(u)(t) | L 2 (Ω) < +∞. ( 41
)
Therefore, we can define the function ψ as follows, ∀α > 0:

ψ : E α → E α
Now, we shall show that ψ is a contraction on E α : let u, v ∈ E α and t ≥ 0, we have

| ψ(u)(t) -ψ(v)(t) | L 2 (Ω) = | � t 0 S(t -s)(f (u(s)) -f (v(s)))ds | L 2 (Ω) ≤ M � t 0 | u(s) -v(s) | L 2 (Ω) ds. ( 42
)
Hence

∀t > 0, e -αt | ψ(u)(t) -ψ(v)(t) | L 2 (Ω) ≤ M α | u -v | Eα . (43)
This shows that the operator ψ is a contraction on E α if α > M . Consequently, using the Banach fixed-point theorem [START_REF] Pazy | Semigroups of linear operators and applications to partial differential equations[END_REF] we deduce the existence of a unique fixed point of ψ in E α . This fixed point is the desired solution of the integral equation

u(t) = S(t)u 0 + � t 0 S(t -s)f (u(s))ds (44)
Which concludes the proof. �

The second model

We conclude this section by presenting another new nonlocal model using the nonlocal second derivatives. The main idea is to take into consideration the second order opposite derivatives instead of considering only the nonlocal Laplacian operator. In fact, we consider all the directions presented in Fig. 1. This technique is inspired by the benefits of considering the different directions in a document-image to handle the complexity of text texture in the degraded image. In [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF], Gilboa and Osher derived also the nonlocal second order derivatives. Let u : Ω ⊂ R 2 → R be a function, and w : Ω × Ω → R be a nonnegative and symmetric weight function, then for i = 1, 2:

� ∂ 2 u ∂x 2 i � w := ∇ w (∇ w (u).∇ w (x i )).∇ w (x i ) = ∇ w ( � Ω (u(y) -u(x))w(x, y)(y i -x i )dy).∇ w (x i ) (45)
Following the same idea, we establish the following text enhancement model to extract the variable illumination from an image:

           ∂u ∂t = max(0, ( ∂ 2 u ∂x 2 1 ) w , ( ∂ 2 u ∂x 2 2 ) w , ( ∂ 2 u ∂x 2 1 ) w + ( ∂ 2 u ∂x 2
2

) w + 2( ∂ 2 u ∂x 1 x 2 ) w , ( ∂ 2 u ∂x 2 1 ) w +( ∂ 2 u ∂x 2 2 ) w -2( ∂ 2 u ∂x 1 x 2 ) w ) in Ω × (0, T ) u(, 0) = u 0 in Ω (46)
In the next section, we describe briefly the proposed algorithms for solving the problems ( 34) and (46).

Numerical implementation

We begin by defining the weight function to measure the similarities between the pixels which is generally based on a gray level intensity vector comparison between two pixels. A pixel is characterized by its neighboring pixels belonging to a neighborhood window centered at its position, and the weight is calculated using the difference between these patches. Indeed, several forms of weight functions can be used [START_REF] Yaroslavsky | Digital picture processing: an introduction[END_REF][START_REF] Buades | Nonlocal image and movie denoising[END_REF]. Note that the choice of the weight function is crucial for the success of non-local methods. In this paper, we propose a non-negative real-valued and symmetric weight function depending on the spatial distance only:

w(x, y) = exp � -| x -y | 2 h � , (47)
where u :

Ω → R is a given image, Ω ∈ R 2 is
a bounded open domain and h is a positive parameter. Before introducing the two proposed algorithms, we consider the following discretization of the nonlocal Laplacian operator:

� w (u i ) = � j∈N i (u j -u i )w i,j j ∈ N i , (48)
where u i is the value of a pixel i (1 ≤ i ≤ N ), N i = {j : |i -j| ≤ r} denotes the neighbors set of the pixel i, and w i,j is the discretized form of the weight function :

w i,j = exp � -| i -j | 2 h � . (49)
To discretize the non-local second derivative operators, we give the discretized version of the non-local gradient:

∇ w (u i ) = (u j -u i ) √ w i,j j ∈ N i , (50)
and the discrete inner product for vectors is given by:

(v 1 .v 2 ) i := � j∈N i (v 1 i,j v 2 i,j ) j ∈ N i . (51)
The second order discretized nonlocal derivatives are given by:

(u i ) �� k,l = ∇ w (∇ w (u i ).∇ w (i k )).∇ w (i l ) k, l ∈ {1, 2}, (52)
where i k is the k th component of the pixel i and

∇ w (i k )(j) = (j k -i k ) √ w i,j j ∈ N i . (53)
Now we are ready to implement the two proposed algorithms based on the discretization form of non-local operators given above. We give firstly the algorithm 1 for solving the problem (34) as follows: Algorithm 1 for solving [START_REF] Morel | Screened poisson equation for image contrast enhancement[END_REF] Input : The acquired image u initialization : We set I 0 = log(u + 1) and choose dt > 0, h > 0 Compute :

I n+1 i = I n i + dt max(0, Δ w (I n i )) i = 1, ..., N Output : The illumination I.
The reflectance R = exp(log(u + 1) -I) In the same way, the algorithm 2 is given for resolving the second model in (46) using the nonlocal second order derivatives instead of the nonlocal Laplacian operator:

Algorithm 2 for solving (46) Input : The acquired image u initialization : We set I 0 = log(u + 1) and choose dt > 0, h > 0 Compute :

I n+1 i = I n i + dt max(0, (I n i ) �� 1,1 , (I n i ) �� 2,2 , (I n i ) �� 1,1 + (I n i ) �� 2,2 + 2 (I n i ) �� 1,2 , (I n i ) �� 1,1 + (I n i ) �� 2,2 -2 (I n i ) �� 1,2 ) i = 1, ..., N Output :
The reflectance image R = exp(log(u + 1) -I)

Numerical results

This section is devoted to the experimental part to test the two proposed nonlocal algorithms. Our aim is to confirm the performance and the robustness of the nonlocal models compared to some local methods in extracting the reflectance of a document-image. The considered images in the reflectance and illumination decomposition problem are assumed to be degraded images acquired under variable illumination. For that, we consider six document-images taken in low light conditions. The first four images are given in grayscale level while the two last ones are color document-images. We compare the two nonlocal algorithms with some available methods in the literature, namely the variational method proposed in [START_REF] Kimmel | A variational framework for retinex[END_REF] (called Variational Retinex), the multiscale Retinex approach [START_REF] Petro | Multiscale retinex. Image Processing On Line[END_REF] (MSRetinex) and two Retinex methods based on partial differential equations proposed in [START_REF] Morel | Screened poisson equation for image contrast enhancement[END_REF][START_REF] Limare | Retinex poisson equation: a model for color perception[END_REF] (Screened Poisson Equation) and (PDE Retinex) respectively. Figs. 234567show the obtained results using the two proposed algorithms compared to the other methods. Throughout this section, for the proposed methods, we set h = 80, while the numerical results by the variational Retinex method are reported with α = 0.0001 and β = 0.1, we set the threshold parameter t = 4 for the PDE-Retinex model, and λ = 0.0001 the tradeoff parameter for the Screened Poisson Equation, while the parameters in the implementation of Multiscale Retinex are chosen as follows: Note that for the first four grayscale images, we consider also a binarization step to show better the estimated text. In the following, we investigate a quantitative comparison in terms of OCR using the restored images in the previous tests apart the handwriting images. Optical Character Recognition tools are efficiently used to convert an image file or PDF file into editable and searchable text format. However, there are many limitations of OCR that can lead to inaccurate or missing text, the OCR engines fail usually to take into account blurry documents, colored paper, mathematical formulas or handwritten text.

The table 1 compares the character recognition rates for the grayscale images using three different types of OCR tools, the "Abbyy FineReader 12"1 , the software "Tesseract"2 and the free web-based optical character recognition software "OnlineOCR.net"3 that return a text with the best confidence rating. The rate "%" is obtained by a careful calculation of the correct characters from the OCR output. The

Conclusion

In the present paper, an existence result for the local PDE's text enhancement model proposed in [START_REF] Buades | Nonlocal image and movie denoising[END_REF] was given, then, based on the nonlocal operators properties, two nonlocal models were proposed to estimate the reflectance image. Moreover, the two nonlocal algorithms are simple, easy to implement and give quite satisfactory results compared to some available methods, which gives confidence in the efficiency of our models. Beyond the two proposed nonlocal PDEs, we could also extend these two models to other PDEs with more regularity properties.

Figure 1 :

 1 Figure 1: The mean at the neighboring pixels of the image u.

  then u ≤ v in ]0, T [×Ω Moreover, for each u 0 ∈ C(Ω), there exists a unique continuous viscosity solution u ∈ C([0, T [×Ω) of (6).

Theorem 2 .

 2 Let Ω be a bounded open subset of R 2 . Then, the problem (34) admits a unique generalized global solution u ∈ C([0, +∞[; L 2 (Ω)).

parameter N σ 1 σ 2 σ 3

 3

Figure 2 :

 2 Figure 2: From top to bottom, (a) the original image 'Proposition', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) Screened Poisson Equation [34], (e) PDE Retinex [29], (f) local laplacian model [32], (g) the local derivatives model [32], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

Figure 3 :

 3 Figure 3: From top to bottom, (a) the original image 'Exercise', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) Screened Poisson Equation [34], (e) PDE Retinex [29], (f) local laplacian model [32], (g) the local derivatives model [32], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

Figure 4 :

 4 Figure 4: From top to bottom and from left to right, (a) the original image 'Phrase', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) ScreenedPoisson Equation[START_REF] Morel | Screened poisson equation for image contrast enhancement[END_REF], (e) PDE Retinex[START_REF] Limare | Retinex poisson equation: a model for color perception[END_REF], (f) local laplacian model[START_REF] Mahani | Text enhancement by pde's based methods[END_REF], (g) the local derivatives model[START_REF] Mahani | Text enhancement by pde's based methods[END_REF], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

Figure 5 :

 5 Figure 5: From top to bottom and from left to right, (a) the original image 'Text', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) Screened Poisson Equation [34], (e) PDE Retinex [29], (f) local laplacian model [32], (g) the local derivatives model [32], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

Figure 6 :

 6 Figure 6: From top to bottom and from left to right, (a) the original image 'Exercise 7', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) Screened Poisson Equation [34], (e) PDE Retinex [29], (f) local laplacian model [32], (g) the local derivatives model [32], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

Figure 7 :

 7 Figure 7: From top to bottom and from left to right, (a) the original image 'Board', the recovered reflectance by (b) MSRetinex [40], (c) Variational Retinex [25], (d) Screened Poisson Equation [34], (e) PDE Retinex [29], (f) local laplacian model [32], (g) the local derivatives model [32], (h) the first proposed nonlocal laplacian model, (i) the second proposed nonlocal derivatives model

  table proves, as might be expected, and identically to the visual comparison in the six different examples, that the two non-local approaches can better recover the degraded text compared to the other approaches, providing then promising numerical results. In comparison, considering the greyscale images, we established that the two proposed methods can reduce the effects of the nonuniform illumination to a greater degree than the other methods. The model based on the nonlocal second derivatives seemed to enhance better color document-images than the first nonlocal model. Our two nonlocal methods, while different, performed comparably at estimating the reflectance part.

Table 1 :

 1 OCR accuracy results

	Image	Method	FineReader % OnlineOCR % Tesseract %
		the original	72.11	77.88	0
		MSRetinex	82.69	95.19	0
		Var. Retinex	95.19	95.19	82.69
		Screned Poisson Equ	91.35	82.96	0
	Proposition	PDE Retinex	0	0	0
		L.Laplacian Model	96.15	98.08	75
		L.Derivatives Model	96.15	98.08	80.77
		The proposed Model 1	97.11	98.08	75
		The proposed Model 2	98.08	96.15	76.92
		the original	0	0	0
		MSRetinex	76.67	97.5	0
		Var. Retinex	98.33	95.83	10
		Screned Poisson Equ	98.33	100	0
	Exercise	PDE Retinex	0	0	0
		L.Laplacian Model	95.83	93.33	37.5
		L.Derivatives Model	93.33	95	40.83
		The proposed Model 1	100	100	82.5
		The proposed Model 2	99.17	97.5	82.5
		the original	0	0	0
		MSRetinex	95	100	0
		Var. Retinex	100	95	60
		Screned Poisson Equ	35	0	0
	Phrase	PDE Retinex	0	60	0
		L.Laplacian Model	45	60	50
		L.Derivatives Model	50	60	55
		The proposed Model 1	100	100	100
		The proposed Model 2	100	100	100
		the original	39.2	62.4	0
		MSRetinex	92	95.2	43.2
		Var. Retinex	99.19	100	89.6
		Screned Poisson Equ	56.8	64.8	0
	Text	PDE Retinex	18.4	31.2	0
		L.Laplacian Model	68.8	99.19	85.6
		L.Derivatives Model	68.8	99.19	86.4
		The proposed Model 1	99.19	100	98.39
		The proposed Model 2	100	100	98.39

Preprint submitted to Journal-of-the-Franklin-instituteJune 7, 2018

https://en.wikipedia.org/wiki/ABBYY FineReader

https://en.wikipedia.org/wiki/Tesseract (software)

https://www.onlineocr.net/

Conflict of interest

The authors declare that they have no conflict of interest.