Nathalie Villa
email: villa@univ-tlse2.fr

Martin Paegelow

Maria T Camacho Olmedo

Laurence Cornez

Frédéric Ferraty

Louis Ferré

Pascal Sarda

Various Approaches for Predicting Land Cover in Mountain Areas

Keywords: polychotomous regression modelling, multilayer perceptron, classification, prediction, comparison

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

PREDICTING LAND COVER

From the sketch maps made by geographers or from the analysis of satellite images or aerial photographs, we can build land cover maps for a given country which can be rather precise: the studied area is then cut into several squared pixels whose sides are about 20 meters long and whose land cover is known on various dates. The type of land cover can be chosen from a pre-determined list: coniferous forests, deciduous forests, scrubs, . . . Here, we are not interested in making such maps (for satellite data analysis, see (Cardot et al., 2003)). Our purpose is to contruct a simulated land cover map at a given future date, by the use of land cover maps at older dates and of other environmental variables; on a geographical point of view, 1 prospective simulations have a great interest to help the local administrations to develop these mountain areas. The idea is then to compare different approaches in order to confront their ability to be generalized to various mountain areas.

For a given pixel, determined by its spatial coordinates, latitude (i) and longitude (j), the value of the land cover on date t, c i,j (t), is a categorical random variable depending on several variables:

• the land cover of this pixel on previous dates: c i,j (t -1), . . . , c i,j (t -T ) (time serie of length T );

• the land covers of the neighbouring pixels on previous dates: V i,j (t -1), . . . , V i,j (t -T ), where V i,j (tτ ) is a set of values of land cover on date tτ for the pixels in a neighbourhood of the pixel (i, j) (vectorial time serie);

• some environmental variables: for example, the elevation, the aspect, the proximity of roads and villages, . . . : Y 1 i,j , . . . , Y p i,j . We face here a problem of classification in which the predictors are both qualitative and quantitative and are also highly dependent (spatial time process). To solve this question, we propose to use and to compare two wellknown statistical approaches with the empirical geographic method (namely the GIS, Geographic Information System). The first of these methods is a generalized linear model in which we estimate the parameters of the model by maximizing a log-likelihood type criterion. The second one uses a supervised multilayer perceptron. By confronting these various approaches, we expect to give ideas in order to improve the GIS approach.

A comparison of these two approaches was done on two little areas: the "Garrotxes" ("Pyrénées Orientales", south west of France) and the "Alta Alpujarra Granaderia" (Sierra Nevada, Spain) where several surveys of the land cover were done at various dates. We confronted the various scenarii constructed with the real maps.

In the following, we describe the data more precisely (section 2) and present the two approaches (section 3). Then we present how we applied these methodologies on these data sets (section 4) and finally, we compare the results obtained by analyzing the advantages and the limits of the models (section 5).

DESCRIPTION OF THE DATA SETS

The areas under study stand in the moutains "Pyrénées" for the Garrotxes and Sierra Nevada "Alta Alpujarra". A big drift from the land has led to the desertion of the land under cultivation and the recovery of the fields by scrubs and forests. There is almost none human action on these areas. The aridity of the climate explains a much slower dynamic in the spanish area than in the Garrotxes: we count 3 times less pixels changing in the Alta Alpujarra than in the Garrotxes. On the contrary, the french area is considered, at least on a geographical point of view, as a dynamic area and it is then more difficult to predict the land cover.

We are given quantitative and qualitative informations through maps divided into pixels: about 241 000 pixels for the French area and 560 000 for the Spanish one (which is much bigger). For each pixel, we know:

• a categorical variable which is the land cover at different dates: 3 dates (1980, 1990 and 2000) were avalaible for the Garrotxes and 4 dates (1957, 1974, 1987 and 2001) for the Alta Alpujarra. As the land cover evolution is very slow in the Sierra Nevada (less than 25% of the pixels had changed their value between 1957 and 2001), these dates were considered as equidistant, according to geographers opinion. This categorical variable was taken from a list of several choices (8 for the Garrotxes and 9 for the Alta Alpujarra) which are of classical use in geography. These data were used to make maps of the studied area (see Figure 1);

[Figure 1 about here.]

• several environmental variables; some of them are of numeric type (the elevation, the slope, the aspect, the distance of roads and villages,. . . ) and others are of categorical type (forest and pasture management: governmental or not ? ground geological type, . . . ). The environmental variables were not the same for the Garrotxes and the Alta Alpujarra (see Figure 2 for examples of environmental variables); all these environmental variables kept the same value at all dates.

[Figure 2 about here.]

PRESENTATION OF THE TWO APPROACHES

Geographers usually estimate the land cover evolution by an empirical method which allows to introduce some expert knowledge. The so-called GIS (Geographic Information System) approach is time expensive and necessitates precise knowledge on the geographic constraints of the area under study. Roughly speaking, the method consists in two steps: at first one computes time transition probabilities for each land cover type whereas, in a second step, one uses spatial constraints (introduced by an expert) for "smoothing" the maps obtained at the first step (see (Paegelow et al., 2004) or (Paegelow and Camacho Olmedo, 2005) for further details on GIS for these data sets). In order to propose automatic alternatives to the GIS, which can take in the same model the spatio-temporal nature of the problem, two approaches have been adapted to estimate the evolution of the land cover: the first one, polychotomous regression modelling, is a generalized linear approach based on the maximum log-likelihood method. The second one, multilayer perceptron, is a popular method which has recently proved its great efficiency to solve various types of problems.

The idea is to confront a parametric linear model with a non parametric one to provide a collection of automatic statistical methods for geographers. They both have concurrent advantages that have to be taken into account when choosing one of them: the polychotomous regression modelling is faster to train than multilayer perceptrons, especially in high dimensional spaces and does not suffer from the existence of local minima. On the contrary, multilayer perceptrons can provide nonlinear solutions and are then more flexible than the linear modelling; moreover, both methods are easy to implement even for non statisticians through the pre-made softwares (for example, "Neural Network" Toolbox for neural network with Matlab).

THE MODEL

Let us now describe the statistical setting more formally. We note X i,j (t) the vector of variables that could explain the value of the land cover for a given pixel (i, j) on date t. We suppose that the time dependence is of order 1; then, X i,j (t) contains:

• for the time series: the value of the land cover for the pixel (i, j) at the previous time t -1;

• for the spatial aspect: the frequency of each type of land cover in the neighbourhood of pixel (i, j) on the previous date. Then, the shape and the size of the neighbourhood had to be chosen. For the shape, we had many choices: the simpler one was a square neighbourhood or a star-shaped neighbourhood around the pixel (i, j); the most sophisticated could use the slope to better take into account the morphological influences of the land. For the size of the neighbourhood, we had to find at which distance a pixel could influence the land use of pixel (i, j).

Moreover, for the multilayer perceptrons, in order to respect the spatial aspect of the problem, we weighted the influence of a pixel by a decreasing function of its distance to the pixel (i, j) (see Figure 3).

[Figure 3 about here.]

• environmental variables (slope, elevation, . . . ).

Let us repeat that c i,j (t) is the land cover for a given pixel on date t. We note C 1 , . . . , C K the different types of land cover. Then, for every k = 1, . . . , K, we try to estimate the probability P (c i,j (t) = C k |X i,j (t)) that the pixel (i, j) has a land cover equal to C k given the vector X i,j (t); thus, the model is of the following form :

P (c i,j (t) = C k |X i,j (t)) = f k (X i,j (t)).
(1)

Once a model was chosen through f k , these probabilities were estimated by the way of a multi-layer perceptron or a generalized linear model and we predicted the type of land cover, c i,j (t), by the rule of maximum:

argmax k=1,...,K P (c i,j (t) = C k |X i,j (t)).
In both approaches, we estimated f k thanks to a training sample. To that end, we have collected the values of the predictors and of the land cover for many pixels on various dates (see next section for more details); the observations are denoted by (X (1) , c (1) ), . . . , (X (N ) , c (N ) ).

The time and spatial aspects are taking into account together both by the polychotomous regression modelling and by the multilayer perceptron and the land cover prediction is performed in a single estimation procedure. This is not the case for the usual GIS approach which is performed in two steps: it first estimates the land cover probability by modelling a time serie and it then introduces a spacial smoothing with environmental constraints.

POLYCHOTOMOUS REGRESSION MODELLING

When we wish to predict a categorical response given a random vector, a useful model is the multiple logistic regression (or polychotomous regression) model (Hosmer and Lemeshow, 1989). A smooth version of this kind of method can be found in (Kooperberg et al., 1997). Applications of these statistical techniques to several situations such as in medicine or for phoneme recognition can be found in these two works. Their good behaviour both on theoretical and practical grounds have been emphasized. In our case, where the predictors are both categorical and scalar, we then have the derived model below.

Let us note, for k = 1, . . . , K

θ (C k |X i,j (t)) = log P (c i,j (t) = C k |X i,j (t)) P (c i,j (t) = C K |X i,j (t)) .
Then, we get the following expression

P (c i,j (t) = C k |X i,j (t)) = exp θ (C k |X i,j (t)) K k ′ =1 exp θ (C k ′ |X i,j (t))
.

(2)

Now, to estimate these conditional probabilities, we use the parametric approach to the polychotomous regression problem, that is the linear model

θ (C k |X i,j (t)) = α k + c∈V i,j (t-1) K l=1 β kl 1 1 [c=C l ] + p r=1 γ kr Y r i,j , (3)
where we recall that V i,j (t -1) are the values of the land cover in the neighbourhood of the pixel (i, j) on the previous date t -1 and (Y r i,j ) r are the values of the environment variables. Let us call δ = (α 1 , . . . , α K-1 , β 1,1 , . . . , β 1,K , β 2,1 , . . . , β 2,K , . . . , β K-1,1 , . . . , β K-1,K , γ 1,1 , . . . , γ 1,K , . . . , γ K-1,1 , . . . , γ K-1,p ), the parameters of the model to be estimated. We have to notice that since θ (C K |X i,j (t)) = 0, we have α K = 0, β K,l = 0 for all l = 1, . . . , K, and γ K,r = 0 for all r = 1, . . . , p. We now have to estimate the vector of parameters δ. For that end, we use a penalized likelihood estimator which is performed on the training sample. Let us write the penalized log-likelihood function for model (3). It is given by

l ε (δ) = l(δ) -ε N n=1 K k=1 u 2 nk , (4)
where the log-likelihood function is

l(δ) = log N n=1 P δ c (n) |X (n) . (5)
In this expression, P δ (c (n) |X (n) ) is the value of the probability given by ( 2) and (3) for the observations (X (n) , c (n) ) and the value δ of the parameter.

In expression (5), ε is a penalization parameter and, for k = 1, . . . , K,

u nk = θ δ (C k |X (n) ) - 1 K K k ′ =1 θ δ (C k ′ |X (n)
). Our penalized likelihood estimator δ ε satisfies:

δ ε = argmax δ∈ R M l ε (δ),
where M = K 2 + (K -1) * p -1 denotes the number of parameters to be estimated.

As pointed out by (Kooperberg et al., 1997) in the context of smooth polychotomous regression, it is possible that, without the penalty term, the maximization of the log-likelihood function l(δ) leads to infinite coefficients β k,l . In our model it may be the case, for example, when, for fixed k, the value of the predictor is equal to zero for all (i, j). Actually, this "pathological" case cannot really occurs in practice but for classes k with a few number of members, the value of the predictor is low and then a numerical unstability happens when maximizing the log-likelihood. Then, the form of the penalty based on the difference between the value θ δ (C k |X (n) ) for class k and the mean over all the classes has the aim of preventing this unstability by forcing θ δ (C k |X (n) ) to be not too far from the mean. On another side, for reasonable values of ǫ, we can expect that the penalty term does not affect so much the estimation of parameters while it guarantees numerical stability. Finally, numerical maximization of the penalized log-likelihood function is achieved by a Newton-Raphson algorithm.

MULTILAYER PERCEPTRON

Neural networks have a great adaptability to any statistical problems and especially to overcome the difficulties of non linear problems even if the predictors are highly correlated; thus it is not surprising to find them used in the chronological series prediction ((Bishop, 1995), (Lai and Wong, 2001) and (Parlitz and Merkwirth, 2000)). The main interest of neural networks is their ability to approximate any function with the desired precision (universal approximation): see, for instance, (Hornik, 1991).

Here we propose to estimate, in model ( 1), the function f k in the form of a multilayer perceptron with one hidden layer (see Figure 4), ψ, which is a function from R q to R that can be written, for all x in R q , as

ψ w (x) = q 2 i=1 w (2) i g x, w (1) i + w (1) i,0 ,
where q 2 in N is the number of neurons on the hidden layer, (w

(1) i ) i=1,...,q 2 (respectively (w (2) i ) i=1,...,q 2 , (w (1)
i,0 ) i=1,...,q 2 ) are in R q (resp. R) and are called weights of the first layer (resp. weights of the second layer, bias) and where g, the activation function, is a sigmoïd; for example, g(x) = 1 1+e -x .

[Figure 4 about here.]

Then, the output of the multilayer perceptron is a smooth function (here it is indefinitly continuous and derivable) of its input. This property ensures that the neural network took into account the spatial aspect of the data set, since two neighbouring pixels have "close" values for their predictor variables.

To determine the optimal value for weights w = ((w

(1) i ) i , (w (2) i ) i , (w (1)
i,0 ) i ), we minimized, as it is usual, the quadratic error on the training sample: for all k = 1, . . . , K, we chose

w k opt = argmin w∈R q 2 (q+2) N n=1 c (n) k -ψ k w (X (n) ) 2 , (6)
where c (n) and the categorical data in X (n) are written on a disjunctive form. This can be performed by classical numerical methods of the first or the second order (such as gradient descent or conjugate gradients, . . . ) but faces local minima problems. We explain in section 4 how we overcome this difficulty. Finally, (White, 1989) gives many results that ensure the convergence of the optimal empirical parameters to the optimal theoretical parameters.

PRACTICAL APPLICATION TO THE DATA SETS

In order to compare the two approaches, we applied the same methodology: we first determined the optimal parameters for each approach (training step, see below) and then, we used the first maps to predict the last one and compared the errors to real map (comparison step, see section 5).

As usual in statistical methods, there are two stages in the training step: the estimation step and the validation step.

• The estimation step consists in estimating the parameters of the models (either for the polychotomous regression or the neural network);

• The validation step allows us to choose, for both methodologies, the best neighbourhood, for polychotomous regression, the penalization parameter and, for neural network, the number of neurons on the hidden layer. Concerning the neighbourhood, we only considered square shapes so choosing a neighbourhood is equivalent, in our procedure, to determine its size.

For the Sierra Nevada, we saw that large areas are constant, thus we only used the pixels for which one neighbour, at least, has a different land cover. These pixels are called "frontier pixels"; the others were considered as constant (see Figure 5). For the generalized linear model, we used the whole frontier pixels of the 1957/1974 maps for the estimation set and the whole 1974/1987 maps for the validation set. We then constructed the estimated 2001 map from the 1987 one. For the multilayer perceptron, we reduced the training set size in order not to have huge computational times when minimizing the loss function. Then, estimation and validation data sets were chosen randomly in the frontier pixels of the 1957/1974 and 1974/1987 maps.

[Figure 5 about here.] For the Garrotxes data set, due to the fact that we only had got 3 maps and much less pixels, we had to use the 1980/1990 maps for the estimation step (only their frontier pixels for the MLP) and the whole 1990/2000 ones for the validation step. This led to a biased estimate when constructing the 2000 map from the 1990 map but, as our purpose is to compare two models and not to make significant the error rate, we do not consider this bias as important.

POLYCHOTOMOUS REGRESSION

• The estimation step produces the estimated parameter vector δ ε of the parameters δ ε of model (3) for given neighbourhood and penalization parameter ε. This step was repeated for various values concerning both neighbourhood and penalization parameter.

• Validation step:

Once given an estimated parameter vector δ ε = ( α 1 , . . . , α K-1 , β 1,1 , . . . , β 1,K , β 2,1 , . . . , β 2,K , . . . , β K-1,1 , . . . , β K-1,K , γ 1,1 , . . . , γ 1,p , . . . , γ K-1,1 , . . . , γ K-1,p ), the quantities

P (c i,j (t) = C k |X i,j (t)) = exp θ (C k |X i,j (t)) K k ′ =1 exp θ (C k ′ |X i,j (t))
, were calculated, for all k = 1, . . . , K, with

θ (C k |X i,j (t)) = α k + c∈V i,j (t) K l=1 β kl 1 1 [c=C l ] + p r=1 γ kr Y r i,j .
At each pixel (i, j) for the predicted map on date t, we affected the most probable vegetation type namely the C k which maximizes

P (c i,j (t) = C k |X i,j (t)) k=1,...,K
.

Programs were made using R programm (see (R Development Core Team, 2005)) and are avalaible on request.

MULTILAYER PERCEPTRON

We used a neural network with one hidden layer having q 2 neurons (where q 2 is a parameter to be calibrated). The inputs of the neural network were:

• For the time series, the disjunctive form of the value of the pixel;

• For the spatial aspect, the weighted frequency of each type of land cover in the neighbourhood of the pixel;

• the environmental variables.

The output was the estimation of the probabilities (1). The estimation was also made in two stages:

• The estimation step produces the estimated weights as described in (6) for a given number of neurons (q 2 ) and a given neighbourhood. For this step, the neural network was trained with an early stopping procedure which allows to stop the optimization algorithm when the validation error (calculated on a part of the data set) is starting to increase (see (Bishop, 1995)). This step was repeated for various values of both neighbourhood and q 2 .

• Validation step: once an estimation of the optimal weights was given, we chose q 2 and the size of neighbourhood, as for the previous model. Moreover, in order to escape the local minima during the training step, we trained the perceptrons many times for each value of neighbourhoud and of q 2 with various training sets; the "best" perceptron was then chosen according to the minimization of the validation error among both the values of the parameters (size of the neighbourhoud and q 2 ) and the optimization procedure results.

Programs were made using Matlab (Neural Networks Toolbox, see (Beale and Demuth, 1998)) and are avalaible on request.

COMPARISON AND DISCUSSION

The validation step led to select the following parameters (Table 1):

[Table 1 about here.]

After the two models were trained, we built the predicting map on date 2000 (Garrotxes data set) and 2001 (Alta Alpujarra data set). The performances of the two models were compared with a GIS approach.

For the Garrotxes data set, the results are summarized in Table 2 and the frequencies of errors for each land cover type were calculated on the pixels which are really of this land cover type. We focus on the 6 more frequent land cover types, since the number of agriculture pixels tends to zero. In Figure 6, we can see the three predictive maps given by our approaches and the GIS approach that can be confronted with the real map.

[Table 2 about For the Alta Alpujarra data set, the results are summarized in Table 3 (land cover type under 5 % of the area have been omitted). Predicted maps and real maps are compared in Figure 7.

[Table 3 about here First of all, the predictive maps provided by the two statistical methods are coherent, smooth and close to reality. This can also be shown through the good error rates (about 25 % -27 % for the Garrotxes data set and 9 % -12 % for the Alta Alpujarra) which are clearly a good performance considering the poverty of the data (we only had got 3 or 4 dates to train the models). Furthermore, the striking fact is that the "automatic" statistical approaches did as well (Garrotxes data set) or even much better (Alta Alpujarra) than the guided GIS approach. This is an interesting point in order to help improving the classical geographical approach to predicting land cover, and better understand the environmental changes in time and space. Moreover, the "automatic" statistical methods were much faster than the GIS as they do not use any expert knowledge which takes a long time to be modelized and needs to be remade for each area. On the contrary, the polychotomous regression modelling and the multilayer perceptron approaches did not lead on these data sets to significant differences. The first method was much faster to train and it was then quite attractive to use it. However, we think that, on a general point of view, the greater flexibility of multilayer perceptron could be usefulness to predict land cover for other data sets where a parametric model could fail.

The main advantage of the automatic statistical approaches is in the fact that they simultaneously take into account the spatio-temporal aspect of the problem and also the environmental variables. GIS works in two steps: it first predicts the number of pixels for each land cover type by a simple temporal model and then takes into account the spatial aspect and the environmental variables to allocate these pixels spatially. This could partially explain that GIS had worse performances for the Alta Alpujarra data set, as the coniferous reforestation used to be important in the 60's and has then be given up. This led the GIS to predict, in the 2001 map, much more coniferous reforestation pixels than in the real map: 18.8 % of the pixels were predicted in the coniferous reforestation type against 7.9 % for the multilayer perceptron, 9.6 % for the polychotomous regression modelling and 9.2 % for the real map. Then GIS approach had a much lower error rate on the coniferous reforestation land cover type but a bigger one for the other ones.

Finally, looking further in the missclassification rates for the various land cover types, we can see that the most dynamic land cover type were harder to train: this is the case, for instance, for the scrubs in the Garrotxes area where they tended to grow fast and became deciduous forests; this is also the case, in the Alta Alpujarra for the fallows and irrigated croplands because agricultural lands were tending to be left. These dynamics could be better predicted by adding pertinent informations for these kinds of land cover types (density of the scrubs, for example, can help knowing if they can, or not, become forests).

CONCLUSION

Finally, this work shows the great potential of the two statistical models in predictive prospection on geographical data. These models had as good performances as GIS approach and we can hope that a combination of the two points of view (statistics and GIS) can improve the land cover predictions: the empirical first step of the GIS could be improved by being replaced by one of these statistical approaches. This issue, that is of big interest for geographers, is still under study as the GIS approach was performed through pre-made programs and has then to be totally re-though to that aim.

Another aspect that has to be worked on is the form of the data: for example, we underlined that an information on the density of the scrubs is needed to better understand their evolution. This could help geographers to better understand what is of interest for predicting the land cover evolution for their future studies.

Version 3 .

 3 The Matworks Inc.(2)Bishop, C. (1995). Neural Networks for Pattern Recognition. New York: Oxford University Press.(3) Cardot, H.,Faivre, R., Goulard, M. (1993). Functional approaches for predicting land use with the temporal evolution of coarse resolution remote sensing data. Journal of Applied Statistics, 30: 1185-1199. (4) Hornik, K. (1991). Approximation capabilities of multilayer feedforward networks. Neural Networks, 4(2): 251-257. (5) Hosmer, D., Lemeshow, S. (1989). Applied logistic regression. New York: Wiley. (6) Kooperberg C., Bose, S., Stone, J. (1997). Polychotomous regression. Journal of the American Statistical Association, 92: 117-127. (7) Lai, T. and Wong, S. (2001). Stochastic neural networks with applications to nonlinear time series. Journal of the American Statistical Association, 96(455): 968-981. (8) R Development Core Team (2005). R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. (9) Paegelow, M., Camacho Olmedo M.T. (2005). Possibilities and limits of prospective GIS land cover modeling -a compared case study: Garrotxes (France) and Alta Alpujarra Granadina (Spain). International Journal of Geographical Information Science, 19(6): 697-722. (10) Paegelow, M., Villa, N., Cornez, L., Ferraty, F., Ferré, L. and Sarda, P. (2004). Modélisations prospectives de l'occupation du sol. Le cas d'une montagne méditerranéenne. Cybergéo, 295. (11) Parlitz, U., Merkwirth, C. (2000). Nonlinear prediction of spatiotemporal time series. In ESANN'2000 proceedings. Bruges, Belgium: 317-322. (12) White, H. (1989). Learning in artificial neural network: a statistical perspective. Neural Computation, 1: 425-464.

Figure 1 :OUTPUTFigure 4 :

 14 Figure 1: Land cover for the Garrotxes (1980 -left) and for the Alta Alpujarra(1957 -right)

Table 1 :

 1 Parameters selected by the validation step

		Garrotxes Alta Alpujarra
	Poly. regression
	Size of neighbourhood	9	1
	ǫ	10	0.1
	ML perceptron
	Size of neighbourhood	7	4
	q 2	8	30
	perceptron size	19-8-7	35-30-9

ACKNOWLEDGEMENTS

The authors are grateful to the Ministerio de Cencia y Tecnologia who supports this research (Plan nacional de investigacion cientifica, Desarollo e innovacion tecnologica, BIA 2003 01499).

The authors are also grateful to the anonymous referees for their detailed and constructive comments and suggestions which have substantially improved the manuscript.

BIBLIOGRAPHY

(1) Beale, M., Demuth, H. (1998). Neural network toolbox user's guide.

List of Figures

Land cover for the Garrotxes (1980 -left) and for the Alta Alpujarra (1957 -right)


Senior Member, IEEE Ronan Fablet

Franc ¸ois Rousseau

Joint interpolation of multi-sensor sea surface temperature fields using non-local and statistical priors

This work addresses the joint analysis of multi-source and multi-resolution remote sensing data for the interpolation of high-resolution sea surface geophysical fields. As case-study application, we consider the interpolation of sea surface temperature fields. We propose a novel statistical model, which combines two key features: an exemplar-based prior and statistical priors. The exemplar-based prior, referred to as a non-local prior, exploits similarities between local patches (small field regions) to interpolate missing data areas from previously observed exemplars. This nonlocal prior also sets an explicit conditioning between the multi-sensor data. Two complementary statistical priors, namely a prior on the spatial covariance and a prior on the marginal distribution of the high-resolution details, are considered as sea surface geophysical fields are expected to depict specific spectral and marginal features in relation to the underlying turbulent ocean dynamics. We report experiments on both synthetic data and real SST data. These experiments demonstrate the contributions of the proposed combination of non-local and statistical priors to interpolate visually-consistent and geophysically-sound SST fields from multi-source satellite data. We further discuss the key features and parameterizations of this model as well as its relevance with respect to classical interpolation techniques.

I. INTRODUCTION

The ocean surface is monitored with a dense network of satellites. Various satellites record images of multiple ocean parameters at different resolutions. As a peculiar example, sea surface temperature (SST), which is the temperature of the thin ocean's upper layer, is a peculiar example. Microwave (MW) radiometry provides low resolution observation (0.5 • -0.25 • , 25/50km) while infrared (IR) sensor delivers high resolution SST measurements (up to 0.02 • , 1km). SST monitoring is of high interest as SST observations are key observations for a wide range of studies, including weather forecasting, ocean circulation, global warming impacts. R. Fablet is with the Department of Signal and Communications at Institut Telecom/Telecom Bretagne; UMR LabSTICC ; Technopôle Brest-Iroise, Brest, France, email: ronan.fablet@telecom-bretagne.eu. F. Rousseau is with with the Department of Image and Information Processing at Institut Telecom/Telecom Bretagne; UMR LaTIM ; Technopôle Brest-Iroise, Brest, France, email: francois.rousseau@telecom-bretagne.eu [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] DRAFT

Beyond their intrinsic resolution, the different satellite sensors also differ in their sensitivity to the atmospheric conditions, and especially the cloud coverage. Whereas microwave radiometers involve very low missing data rates, IR sensor may result in high missing data rates (up to 90% over several consecutive days in some regions).

Operational level-4 products provide interpolated fields, which typically rely on covariance-driven kriging-based approaches (e.g., [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF], [START_REF] Stark | An operational, high resolution, real time, global sea surface temperature analysis system[END_REF]). The joint exploitation of the multi-source multi-scale data remains however a challenge to produce geophysically-consistent interpolated fields.

In this context, we address here the interpolation of missing data in high-resolution geophysical fields under the assumption that a simultaneous low-resolution observation is available. Such multi-scale interpolation is expected to depict geophysically-consistent features, with a view to conforming to the underlying turbulent dynamics, especially: i) consistent high-resolution textured patterns, ii) non-Gaussian marginals and iii) specific spectral signatures.

We propose a novel model to address this multi-scale interpolation jointly accounting for the above-mentioned constraints. The key idea is to exploit a non-local patch-based framework [START_REF] Freeman | Markov Random Fields for Super-resolution and Texture Synthesis[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF], [START_REF] Tartavel | Variational Texture Synthesis with Sparsity and Spectrum Constraints[END_REF]. Such exemplar-based representation provides an implicit texture model, which directly exploits available observations to interpolate missing data areas.

II. PROBLEM STATEMENT AND RELATED WORKS

The state-of-the-art techniques for missing data interpolation in earth remote sensing mainly involve two main categories of approaches. The most popular techniques, by far, involve covariance-based geostatistical filtering.

This category includes kriging methods as well as optimal interpolation schemes. They are widely used both in land remote sensing [START_REF] Qi | The assimilation of Jason-2 significant wave height data in the north Indian ocean using the ensemble optimal interpolation[END_REF], [START_REF] Wang | Improving LAI mapping by integrating MODIS and CYCLOPES LAI products using optimal interpolation[END_REF] and ocean remote sensing [START_REF] Qi | The assimilation of Jason-2 significant wave height data in the north Indian ocean using the ensemble optimal interpolation[END_REF], [START_REF] Tandeo | Linear Gaussian State-Space Model with Irregular Sampling: Application to Sea Surface Temperature[END_REF], including applications in most operational products [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF], [START_REF] Stark | An operational, high resolution, real time, global sea surface temperature analysis system[END_REF]. These approaches rely on the definition and/or estimation of the spatial covariance of the field to be interpolated. The interpolation is stated as the computation for each grid point of the conditional mean of the field given the observed grid points and the considered covariance model. In general, this geostatistical interpolation simultaneously performs the interpolation and some smoothing of the original data. By nature, through the parameterization of the covariance model, it involves some trade-off between the fidelity to the observed data and the reconstruction of a spatially-consistent field. Hence, interpolated fields typically involve a noticeable smoothing of the truly observed data. The same difficulty arises from the second category of approaches which involves matrix factorization [START_REF] Alvera-Azcarate | Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions: application to the Adriatic Sea surface temperature[END_REF], especially EOF-based (Empirical Orthogonal Functions, also referred to as PCA in the image analysis community) approaches, which relies on the projection onto a precomputed basis functions. By contrast, we develop here a different point of view. We rather aim at reconstructing or emulating a spatially-and geophysically-consistent field with a strong fidelity to the observed data. This is referred to as image inpainting in the image processing literature [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], [START_REF] Le Meur | Hierarchical Super-Resolution-Based Inpainting[END_REF]. This is also viewed as a mean to truly handle multi-source data acquired at different resolutions. Model-based interpolation, such as polynomial or spline-based interpolation, may be regarded as possible solutions [START_REF] Yi | A fast iterative interpolation method in f-k domain for 3-D irregularly sampled GPR data[END_REF], [START_REF] Zhang | A simple method for soil moisture determination from LST-VI feature space using nonlinear interpolation based on thermal infrared remotely sensed data[END_REF], but they cannot deal with textural patterns.

In the image processing field, image inpainting and super-resolution have been particularly active research topics over the last decade [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], [START_REF] Freeman | Markov Random Fields for Super-resolution and Texture Synthesis[END_REF], [START_REF] Glasner | Super-resolution from a single image[END_REF], [START_REF] Romano | Single Image Interpolation Via Adaptive Nonlocal Sparsity-Based Modeling[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF]. A variety of models and algorithms have been introduced, including January 28, 2016 DRAFT for instance variational, Bayesian, non-local setting. As geophysical fields may be viewed as textured images, we focus on texture-related models. Among others, we may distinguish two categories of approaches: exemplar-based or patch-based models and explicit statistical models. Exemplar-based schemes state the inpainting and super-resolution issue as the sampling from a collection of local patches (also referred to as a patch dictionnary). From the original work of Efros et al. [START_REF] Efros | Image Quilting for Texture Synthesis and Transfer[END_REF], the formulation of these exemplar-based strategies as non-local models [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Peyre | Non-local Regularization of Inverse Problems[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF] has become the state-of-art for a wide range of inverse problems in image processing. Applications to remote sensing data have also been investigated [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF], including missing data interpolation [START_REF] Lorenzi | Inpainting strategies for reconstruction of missing data in VHR images[END_REF], [START_REF] Yu | An effective interpolation method for MODIS land surface temperature on the Qinghai-Tibet Plateau[END_REF]. By contrast, explicit statistical models rely on specific statistical features such as spatial covariances (or equivalently, spectral power densities) [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF], [START_REF] Byung | Super-resolution texture synthesis using stochastic PAR/NL model[END_REF], marginal distributions [START_REF] Heeger | Pyramid Based Texture Analysis/Synthesis[END_REF], scale-invariance [START_REF] Chainais | Virtual super resolution of textured images using multifractal stochastic processes[END_REF]. Recently, the coupling of patch-based representation and of spectral characteristics has been investigated for texture synthesis.

Regarding ocean surface geophysical fields, such statistical features are particularly relevant as ocean turbulence is expected to resort in specific spectral and non-Gaussian characteristics. Besides, for any given case study region, existing image time series provide the mean for building representative dictionaries of exemplars. Therefore, the combination of non-local models and statistical priors, which has not been investigated in previous work to our knowledge, appears as a natural paradigm to address sea surface field reconstruction from multi-source and multiscale satellite data. This work aims at addressing this issue and our main contributions are two-fold:

• the combination of non-local and statistical priors within a mathematically-sound formulation for the texturebased interpolation of geophysical fields;

• the development and evaluation of non-local models for the joint interpolation of multi-source/multi-resolution SST data.

We detail in the subsequent these contributions.

III. PROPOSED APPROACH

Formally, the considered missing data interpolation is stated as an inverse problem, i.e. the reconstruction of unknown fields X (1) and X (2) from partial observations Y (1) and Y (2) , where the subscripts refers to satellitederived gridded fields from two different sensors. These sensors may differ both in the associated spatial resolution and in the missing data rate. We consider all fields are gridded and interpolated on some reference grid Ω, referred to hereafter as the high-resolution grid. We denote respectively by Ω (1) and Ω (2) the missing data region for each sensor, and, r (1) and r (2) the actual resolution associated with each sensor. Subscript (1) refers to the sensor with the lowest resolution. By convention, we assume below for the sake of simplicity that r (1) ≥ 1 and r (2) = 1. This amounts to considering below that all scales are defined up to the actual resolution grid.

In this section, we describe the proposed framework as follows. We first introduce the considered multi-resolution representation. Second, we introduce non-local and patch-based priors. Third, the considered statistical priors are detailed. Fourth, we present the proposed interpolation model along with the associated numerical scheme.

January 28, 2016 DRAFT

A. Multi-resolution representation

With a view to handling information sources at different resolutions, we consider a wavelet decomposition of the analysed fields [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. It comes to decompose field X over an orthonormal basis of functions {Ψ s,i }, which are scaled and translated versions of a mother wavelet Φ(.):

∀p ∈ Ω, X(p) = s,i X s,i Ψ s,i (p) (1)
where {X s,i } are the wavelet coefficients for field X. Indices s and i refer respectively to a scale index and to a spatial position at scale s. We consider here a dyadic decomposition, such that scale indices are of the form 2 j . For a given field, the computation of the associated wavelet coefficients are issued from the projections onto the basis functions, which are implemented as decimation and filtering steps iterated from the fine scale to the coarsest scale.

It may be noticed that the wavelet decomposition guarantees the exact reconstruction of field X from its wavelet coefficients {X s,i }. We let the reader refer to [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] for further details on the wavelet decomposition. We introduce the associated wavelet-based projection operator P s * W1 for scale s * , which refers to the reconstruction of field X from all scales greater than s *

∀p ∈ Ω, P s * W (X)(p) = s>s * ,i X s,i Ψ s,i (p) (2)
It may be noted that, by construction, field X -P s * W (X) refers to a detail field, denoted by δP s * W (X) which only involves information at scales finer than s * .

For the targeted issue, we denote by s LR the reference low-resolution scale in the dyadic wavelet decomposition, given as the dyadic scale log 2 r (1) the closest to the low resolution r (1)2 . Given that the two sensors observe the same geophysical condition, we assume that the two fields X (1) and X (2) share the same low-resolution projection at scale s LR such that we impose

P s LR W X (1) = P s LR W X (2) (3)

B. Multimodal patch-based interpolation and non-local projection

As stressed in the introduction, patch-based and non-local priors are particularly appealing to develop data-driven image processing techniques, by contrast to model-driven techniques [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF], [START_REF] Freeman | Markov Random Fields for Super-resolution and Texture Synthesis[END_REF], [START_REF] Lorenzi | Inpainting strategies for reconstruction of missing data in VHR images[END_REF]. This appears especially relevant when the definition of explicit analytical a priori models is a complex issue, as for sea surface turbulent dynamics whose characteristics may vary in space and time. We here focus on non-local regularization or projection [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Peyre | Non-local Regularization of Inverse Problems[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF]. Given a field X, it comes to compute a non-local mean P N LM ({A k })(X), with respect to a reference set of patches based on patch similarities. The non-local mean for a patch around a pixel p ∈ Ω is given by

P N LM ({A k })(X)(p) = k∈K w k (X, {A l } , p) A k ( 4
)
where k is a patch index within the reference patch dictionary, A k the k th P × P patch in this dictionary. Typical values for the patch width P are between 5 and 13. w k (I, A k , p) is a weighing factor which weighs the contribution of the k th reference patch A k in the non-local mean at point p. It relies on the similarity between the P × P patch around pixel p in field X, denoted by X(N p ), and reference patch A k . Patch similarities are generally evaluated according to a classical Euclidean distance [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF], such weight w k (X, {A l } , p) resorts to a normalised weight as follows:

w k (X, {A l } , p) = g ( X(N p ) -A k ) l g ( X(N p ) -A l ) (5)
where g(.) is a kernel function. Non-local regularisation [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Peyre | Non-local Regularization of Inverse Problems[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF] typically considers a Gaussian kernel

g(x) = exp(-γ.x 2 )
with γ a smoothing parameter to be defined. In texture synthesis, nearest-neighbor kernels may be preferred to avoid some over smoothing of textured details. Different strategies may also be considered for the combination of the non-local means computed for different image positions. Non-local denoising and segmentation [START_REF] Buades | A non-local algorithm for image denoising[END_REF], [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF] typically involves the computation of the non-local mean for every pixel to derive the non-local image projection. By contrast, in texture synthesis and inpainting [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], [START_REF] Freeman | Markov Random Fields for Super-resolution and Texture Synthesis[END_REF], [START_REF] Le Meur | Hierarchical Super-Resolution-Based Inpainting[END_REF], one usually achieves the non-local projection for partially or non-overlapping patches with a predefined, optimized or randomized selection of the pixels at which the non-local means are actually evaluated (Eq.4). The later, which is the more generic, is selected in this work. We apply the non-local means to randomly selected positions and, for each selected position p, we update the non-local image projection with the (P -B) × (P -B) interior of the non-local mean P × P patch (4)

with B the width at the boundary. It may noted that, in our implementation, we use for each non-local projection a different random sequence of patch positions with a view to ensuring shift invariance [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF].

Here, we apply the non-local setting to the joint interpolation of the two fields X (1) and X (2) associated with two different sensors. Hence, the evaluation of patch similarities naturally relies on the information from each field and the reference patch dictionary {A (1)

k , A (2)
k } is formed by concomitant patch exemplars form each sensor. Non-local projection (4) becomes a multi-modal non-local projection. For the sake of simplicity and readability, we drop the explicit reference to patch dictionary {A (1)

k , A (2)
k } in our notations. For field X (1) , we denote by P N LM,X (2) (X (1) ) its non-local projection given reference patch dictionary A (1)

l , A (2)
l and field X (1) : ∀p ∈ Ω

P N LM,X (2) X (1) (N p ) = k∈K w k X (1) , X (2) , p A (1) k (6)
where the normalized non-local weights involve joint similarities for the two fields w k (X (1) , X (2) , p) ∝ g α (1) X (1)

(N p ) -A (1) k + α (2) X (2) (N p ) -A (2) k (7)
where α (1) and α (2) are scalar values which weigh the relative contribution of each sensor. Regarding the targeted issue, we explore different strategies for the multi-modal non-local projection:

• Projection P N LM,X (2) X (1) : the non-local projection of field X (1) given field X (2) ; January 28, 2016 DRAFT

• Projection P N LM, ∇X (2) δX (1) : the non-local projection of detail field δX (1) given field ∇X (2) . This strategy is motivated by recent observations that the variability of the high-resolution detail δX (1) is driven by the magnitude of the low-resolution gradient ∇X (2) [13].

In the reported experiments, we compare these different strategies and evaluate their pros and cons.

C. Statistical constraints

Geophysical fields are expected to depict specific statistical features [START_REF] Xu | Global variability of the wavenumber spectrum of oceanic mesoscale turbulence[END_REF], [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], [START_REF] Ba | Statistical descriptors of ocean regimes from the geometric regularity of Sea Surface Temperature obser-vations[END_REF], [START_REF] Bernard | Conformal invariance in two-dimensional turbulence[END_REF], which may be used as additional priors to further constrain their reconstruction from partial and noisy observations. In this work, we focus on two categories of statistical features: covariance (or equivalently spectral) features and non-Gaussianity.

Sea surface geophysical processes, especially for the mesoscale range (i.e., for horizontal scales from tens of kilometres to hundreds of kilometres), involve turbulent processes characterised by specific spectral signatures [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], [START_REF] Mccaffrey | Estimates of ocean macro-turbulence: structure function and spectral slope from Argo profiling floats[END_REF]. Their radial power spectral density generally involves power laws [START_REF] Xu | Global variability of the wavenumber spectrum of oceanic mesoscale turbulence[END_REF], [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF]. A typical example is the homogeneous and isotropic turbulence for incompressible flows [START_REF] Kolmogorov | The local structure of turbulence in in-compressible viscous fluid for very large reynolds numbers[END_REF] which resorts to Kolmogorov's 5/3 power law for the power spectral density of the energy of the velocity field. We then set a power law prior on the radial spectral density of geophysical fields X (1) and X (2) for the fine-scale range (i.e., scales lower than r (1) ). Let us denote by R * the reference radial power spectral density and S X,X the power spectral density of field X. The radial power spectral density of field X, denoted by R X is a radial integration of power spectral density S X,X :

R X (r) = ω =r S X,X (w)dw (8)
It may be noted that power spectral density S X,X is the Fourier transform of the spatial covariance of field X. The considered spectral prior comes to a define a reference radial power spectral density R * . In the reported numerical experiments, we consider both a non-paramatric densities and power-law parameterization as:

R * (r) = β • r ξ (9)
where γ and α are the parameters of this spectral prior. In the reported experiments, we calibrate these parameters from the spatial covariance of field X (1) for scales greater than r (1) (typically, greater than 100km for the considered experiments on SST data). Given this prior, the fine-scale spectral constraint for field X (2) amounts to:

R X ( 2) (r) = R * (r), ∀r ∈ [r (2) , r (1) ] (10)
This spectral constraint can be solved for in the Fourier domain through the computation of a transformed field

T R * r (1) (X) such that its Fourier transform T R * r (1) (X) is given by: T R * (X) = R * R X X (11)
Numerically, the Fourier transform is computed here using a windowed FFT. It might be pointed out that the Fourier transform amounts to a linear filtering and a projection onto the manifold of the two-dimensional fields accounting for spectral constraint [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF].

January 28, 2016 DRAFT

In addition to this spectral constraint, ocean turbulence is widely acknowledged to depict non-Gaussian marginals in multiscale decomposition [START_REF] Bernard | Conformal invariance in two-dimensional turbulence[END_REF], [START_REF] Fablet | Random walk models for geometry-driven image super-resolution[END_REF], [START_REF] Nieves | Common turbulent signature in sea surface temperature and chlorophyll maps[END_REF] as considered here. Turbulent flows tend to advect tracers to form fronts or high-gradient structures at different scales. This in turn results in non-Gaussian heavy-tailed marginals of the detail fields in the considered multi-scale representation. Generalized Gaussian distributions GG σ,α appear as relevant, yet simple, parametric models to account for the observed non-Gaussian signatures:

GG σ,α (x) ∝ exp - x 2σ ν ( 12
)
where σ is a scaling parameter and ν a power exponent. For ν = 2, it resorts to a standard zero-mean Gaussian distribution. For a given field X (2) , we fit a generalized Gaussian distribution to the detail field δX (2) -P s LR W X (2) according to a maximum likelihood (ML) criterion. Let us denote by GG σ * ,α * the reference distribution. The marginal constraint then amounts to:

X (2) -P s LR W X (2) ∝ GG σ * ,α * (13)
To solve for this constraint for field X (2) , we proceed as follows. Given current fields X (1) and X (2) , we fit a generalized Gaussian distribution GG σ, α to the current detail field X (2) -X (1) according to a ML criterion. The transformed field T GG * (X (2) ) resorts to:

T GG * (X) = P s LR W X (2)
+ Ψ σ * ,α * Ψ -1 σ, α X (2) (p) -P s LR W X (2) (p) (14)
where Ψ σ,α is the cumulative distribution function of generalized Gaussian distribution GG σ, α, Ψ -1 σ, α the associated inverse cumulative distribution function.

D. Joint interpolation and blending framework

The joint interpolation and blending of two fields X (1) and X (2) from partial observations Y (1) and Y (2) is stated as the following joint non-local projection:

           X (1) = arg min X (1) X (1) -P N LM,X (2) X (1) 2 Ω (1) X (2) = arg min X (2) X (2) -P N LM,X (1) X (2) 2 Ω (2) (15)
January 28, 2016 DRAFT where • 2 Ω stands for the quadratic norm restricted to spatial domain Ω. The above minimization is to be solved for under the four following constraints:

                                                              
Interpolation constraints:

X (1) Ω \ Ω (1) = Y (1) Ω \ Ω (1) X (2) Ω \ Ω (2) = Y (2) Ω \ Ω (2)
Low-resolution constraint:

P s LR W X (1) = P s LR W X (2)
Spectral constraints:

R X (1) = R *
Marginal constraint:

X (2) -P s LR W X (2) ∝ GG σ * ,α * (16)
Problem formulation [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF] refers to the application of the non-local projections to field X (1) given field X (2) . As described in Section III-B, for field X (2) , we consider a different strategy with non-local projection P N LM, ∇X (2) δX (1) applied to the detail field. This formulation refers to the two key contributions of this work, namely:

• the combination of non-local priors to statistical priors through both non-local projections [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF] and spectral and marginal constraints in (16);

• the application of non-local models for the joint interpolation of multi-source/multi-resolution SST data highlighted by both interpolation constraints in [START_REF] Glasner | Super-resolution from a single image[END_REF] and non-local projections (15) at low-and high-resolution.

We summarize in Table ?? in the Supplementary material the different variables and parameters considered in this study.

Symbol Description Default value γ Scale parameter of the Gaussian kernel, Eq.( 5) 1

α (1) , α (2) Relative weights of modality (i) in the non-local weight, Eq.( 7

SUMMARY OF THE PARAMETERS CONSIDERED IN OUR MODEL

To solve for this constrained minimization, we consider an iterative procedure. From given initialization X (1),0 and X (2),0 , it comes to iterate the non-local projections of the current solutions X (1),n and X (2),n and the projection of the resulting fields P N LM,X (2) X (1),k and P N LM,X (1) X (2),k onto the four additional constraints [START_REF] Glasner | Super-resolution from a single image[END_REF].

Numerically, this approach involves three key components:

January 28, 2016 DRAFT

• The projection onto interpolation and low-resolution constraints (first-two constraints in ( 16)): we solve for this projection using a gradient-based descent to retrieve the field accounting for the interpolation constraint while minimizing the quadratic norm to the expected low-resolution field (see Algorithm [START_REF] Alvera-Azcarate | Reconstruction of incomplete oceanographic data sets using empirical orthogonal functions: application to the Adriatic Sea surface temperature[END_REF] in the Supplementary Material for details);

• The projection onto the constraints on the spectral and marginal constraints. Based on transforms ( 14) and [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF], we use a gradient-based descent to jointly solve for these two constraints. Each iteration of this gradientbased descent involves a projection onto interpolation and low-resolution constraints (see Algorithm [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF] in Supplementary Material for details);

• The overall algorithm iterates three main steps: i) non-local projections of fields X (1) and X (2) followed by a projection onto interpolation and low-resolution constraints, ii) projections of resulting field updates onto the marginal and spectral constraints, iii) the update of the reference low-resolution field as the mean of the low-resolution projection of the updated fields X (1) and X (2) as well as the update of reference generalised Gaussian distribution GG σ * ,α * for the high-resolution details of field X (1) (cf. ( 13)). The procedure is described in the Supplementary Material by Algorithm (1).

The initialization of fields X (1) and X (2) from partial observations Y (1) and Y (2) may exploit classical interpolation techniques. In the reported experiments, we consider a initial variational interpolation using the AMLE (Absolutely Minimizing Lipschitz Extension) [START_REF] Caselles | An axiomatic approach to image interpolation[END_REF]. An additional random component is added in missing data area Ω (2) for highresolution field X (2) . We use a white Gaussian noise, such that initial field X (2) involves the expected variance for high-resolution detail X (2) -P s LR W X (2) .

IV. RESULTS

We report numerical experiments to evaluate the proposed framework for two types of data: i) high-resolution SST data from ODYSSEA analysis [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF] for which we simulate missing data and a low-resolution condition, ii) real AMSR-E and SEVIRI daily SST data.

A. ODYSSEA SST data

Considered data and expriments: ODYSSEA SST fields are gridded and interpolated composite fields produced by Ifremer/CERSAT [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF] onto a grid with a 0.02 • spatial resolution (∼ 2km). ODYSSEA analysis combines both microwave sensor data (including TMI and AMSR-E data) and infrared sensor data (including AVHRR, VISSR and SEVIRI data). It involves different processing steps, especially sensor-specific bias correction, sensor intercalibration and optimal interpolation of missing data areas. The selected field (January, 8th 2012) and case study area, off South Africa, are chosen for the the quality of the composite field and the richness of the observed SST structures (e.g., filaments, fronts).

From the considered ODYSSEA SST field, we first evaluate the relevance of the proposed model for the interpolation of a high-resolution field given a known low-resolution condition, i.e. observation Y (1) does not involve missing data and X (1) = Y (1) . We proceed as follows. For a given dyadic downsampling factor κ, here 16

January 28, 2016 DRAFT for the real real (black,-) and interpolated (red, -) fields along with the target generalized Gaussian distribution (magenta, -) and a Gaussian distribution with the same variance (cyan,-). We superimpose image level-lines to illustrate the underlying spatial geometry.

(i.e., r (1) = 2km and r (1) = 32km and s LR = 4), we derive the low-resolution condition X (1) as P s LR W X (2) the low-resolution projection of the considered ODYSSEA field with a Daubechies-8 wavelet transform at scale s LR = log 2 (κ). Regarding high-resolution observation Y (2) , we generate randomly-distributed missing data. With a view to evaluating the performance of the model and the relative importance of the non-local and statistical priors, we may consider a random missing data pattern with large missing data area (as in Fig. 1 with a mean radius of 150 pixels ) or small missing data areas (here, 10-pixel-wide missing data areas).

The considered parameter setting for the non-local priors is as follows. We select 5000 11 × 11 patches off the missing data area. The non-local means is computed for the first five neighbors (Eq.7). The scale parameter γ of the Gaussian kernel is set to 1 (Eq.5) and relative weights for α (1) and α (2) (Eq.7) are set to 0.5 and 0.5 for non-local projection P N LM,X (2) X (1) , and to 0.1 and 0.9 for non-local projection P N LM, ∇X (2) δX (1) . The reference power spectral density R * is the radial spectral density of the high-resolution field. Similarly, the reference marginal distribution GG σ * ,α * is the generalized Gaussian distribution fitted to the detail of the true high-resolution field according to a Maximum Likelihood criterion (σ * = 0.08 and α * = 1.44).

Comparison of stochastic interpolation models:

We report in Fig. 1 the interpolation results for the large missing January 28, 2016 DRAFT 2) a a function of the low-resolution gradient ∇X (1) (center), standard deviation of the angular difference between gradient fields ∇X (1) and ∇X (2) as a function of the low-resolution gradient ∇X (1) (right).

For each subplot, we report the features associated with the real high-resolution field (black,-), the low-resolution field (blue,-), and four different interpolation models, namely using only the covariance and marginal constraints (C+M, magenta,-), using only the non-local projection (NL, light green, -), using the proposed model (NL+C+M, red,-), using the proposed model with the non-local projection applied to detail field 2) rather than field X (2) (NL+C+M+dX, dark green,-).

δP s LR W X (
area example. These interpolation results are issued from the non-local projection of field X (2) given low-resolution field X (1) . As illustrated, the interpolated field depict to the expected spectral and marginal characteristics. With a view to evaluating the relative importance of the non-local an statistical priors, we consider four different texturebased interpolation models:

• a model relying only on the spectral and marginal priors (model C + M )3 ;

• a model using only the non-local prior (model N L) 4 ;

• the proposed model combining the non-local and statistical priors (model N L + C + M );

• the proposed model combining the non-local and statistical priors with the non-local projection applied to detail field δP s LR W X (2) rather than field X (2) (model N L + C + dX). As depicted in Fig. 2, we perform a quantitative comparison of these four models according to six statistical features 5 :

Model C+M Model NL+C+M Model NL+C+M+dX

the radial spectral density, the multifractal spectrum, the marginal distribution of high-resolution details P s LR W X (2) , the distribution of the magnitude of gradient ∇X (2) , the standard deviation of high-resolution detail P s LR W X (2) as a function of the norm of low-resolution gradient ∇X (1) and the standard deviation of the angular differences between low-resolution and high-resolution gradients as a function of the norm of low-resolution gradient ∇X (1) .

These different features are selected to analyse the extent to which we recover the inter-scale characteristics of the reference field. These results clearly stress the relevance of the combination of a non-local prior with covariance and marginal priors. The sole use of a non-local prior leads to high-resolution field which depicts too much fine-scale energy. The combination with the covariance-based prior allows us to further regularize the non-local projections and to match the expected spectral properties. It also leads to a better fit in terms of gradient statistics and inter-scale angular statistics, as the non-local prior exploits an actual conditioning of the high-resolution detail from the lowresolution field. This is further illustrated by the visual inspection of the gradient field of the different interpolation results (Fig. 3). Little difference is observed between the proposed model whether we apply the non-local projection to high-resolution field X (2) or to high-resolution detail field P s LR W X (2) .

Comparison to a kriging-based interpolation: We report a comparison Fig. 4 to a kriging-based interpolation 5 It might be noted that classical mean square error (MSE) statistics do not apply for texture-based interpolation and super-resolution for large subsampling factors [START_REF] Chainais | Virtual super resolution of textured images using multifractal stochastic processes[END_REF], [START_REF] Freeman | Markov Random Fields for Super-resolution and Texture Synthesis[END_REF] as considered here. For textured fields, small spatial translations of the high-resolution details generate large MSE values, so that the interpolant, which minimizes the MSE would resort to a non-textured field as illustrated in Fig. 4 with a kriging-based interpolation January 28, 2016 DRAFT

for the example considered in Fig. 1. We apply the kriging-based interpolation to the detailed field P s LR W X (2) .

The considered covariance model is a stationary isotropic Matern model fitted to the true detailed field P s LR W X (2) [33]. We resort to a Matern covariance with a standard deviation of 0.01, a power exponent of 3 and a spatial range of 6 pixels. Under a Gaussian hypothesis, kriging-based interpolation comes to derive a linear prediction of the missing data from the truly observed ones. This linear prediction refers to the mean of the conditional Gaussian distribution of the the missing data given the truly observed ones. This mean value then minimizes the mean square reconstruction error. For points, which are far from any truly observed data, the linear prediction converges to the mean of the Gaussian process, that is to say that, in terms of mean square error statistics, the best prediction is the mean of the process is no close-enough points is truly observed. This is observed in Fig. 4. Given the spatial range of the fitted covariance model, the linear prediction of the high-resoluton detail for points which are further than 10 pixels from the boundary of the missing data area is close to 0, so that the kriging-based interpolation resorts to a relative mean square error close to 1, which means it cannot bring additional information compared to the low-resolution field. As such, the interpolated field does not depict spectral and statistical patterns similar to the true high-resolution field (Fig. 5). This is intrinsically different from the considered interpolation strategy. We do not aim at minimizing a mean square reconstruction error statistics and retrieve a relative mean square error of 1.3 larger than the kriging-based interpolation due to small spatial misalignements between the true and simulated high-resolution details. Meanwhile, as shown in (Fig. 4 & Fig. 5), we emulate a spatially-consistent field, where the high-resolution details are statistically-consistent at both observed and missing data points. Similar experiments have been performed for random small missing areas with a mean radius of 10 pixels, for which little difference is reported between the kriging-based and proposed interpolation models both in terms of mean square error statistics (0.61 vs. 0.67) and of the statistical features reported in Fig. 5.

B. Joint analysis of AMSR-E and SEVIRI SST data

We consider an application of the proposed framework to the joint analysis of AMSR-E and SEVIRI SST fields.

AMSR-E SST field involves a typical resolution grid of 0.2 • , whereas we consider daily composite SEVIRI field with a 0.1 • grid resolution. From Fig. 6, one can notice that their radial power spectral density differ from about 1 • , with a significantly lower spectral energy of the AMSR-E field. This difference relates to the interpolation and filtering applied to AMSR-E fields. The situation reported in Fig. 6 depicts a low missing data rate for the AMSR-E field, but a much higher missing data rate, especially off South Africa for the SEVIRI field. Besides, due to the associated geostationary orbit, no SEVIRI data is available, whatever the atmospheric conditions, for the extreme south-east area of the case-study region. In a similar way, contrary to SEVIRI data, AMSR-E data do not provide any information close to the sea shore. Regarding the parameterization of the proposed model, we proceed as follows. For the reference radial power spectral density, we consider a power law fitted for scales [0.8 We report in Fig. 6 the interpolation results issued from the proposed model with a non-local projection applied to detail field δP s LR W X (2) . Though the AMSR-E field involve few missing data, we benefit from the greater resolution of the SEVIRI field close to the shore to interpolate the AMSR-E field in these areas. This is particularly January 28, 2016 DRAFT noticeable off South Africa for the Aghulas current, characterized by warm waters going southward along the eastern coast of South Africa (see also, the first row of Fig. 7). It may also be outlined that the constraint of a common low-resolution component between the interpolated AMSR-E and SEVIRI fields result in more coherent scatterplots of the low-vs high-resolution data (Fig. 6, bottom right). From the visual inspection of the zoom onto an area off South Africa, involving an important missing data rate for the SEVIRI field (Fig. 7), the gradient fields clearly stress the greater relevance of the combination of non-local and statistical priors in the south-east of the considered area. Besides, it may be noted that the application of the non-local projection to field X (2) produces unexpected tiling effects. This is due to a very low number of cold patches in the reference dataset. By contrast, applying the non-local projection to field δP s LR W X (2) leads to much more satisfactory interpolation results, as we are provided with more representative patch datasets for the detail field and a conditioning by the low-resolution gradient. The gradient fields also clearly stress the difference between the information carried out by AMSR-E and SEVIRI fields in terms of presence and location of frontal structures, with a much coarser resolution of the AMSR-E fields. We report a second example of joint AMSR-E/SEVIRI interpolations for a date at which AMSR-E field also involves a large missing data area off South Africa (Fig. 8). Following the conclusions drawn from the reported experiments, we only depict the zoom on the region off South Africa and the interpolation results for the proposed model with a non-local projection applied to detail field δP s LR W X (2) . This example shows that we can interpolate For each column, we depict the observed data, the interpolation results using the proposed model with a non-local projection applied to field X (2) (Model NL+C+M), the interpolation results using the proposed model with a non-local projection applied to detail field δP s LR W X (2) (Model NL+C+M+dX), the interpolation results using only the statistical priors (Model C+M).

visually-consistent fields both at AMSR-E and SEVIRI resolutions and jointly benefit from the two types of data to jointly improve the interpolation of the low-resolution and high-resolution fields.

V. DISCUSSION

In this work, we have presented a novel model for the joint interpolation and blending of multi-source and multiresolution remote sensing images. As case-study, we have focused on sea surface temperature. The key contributions are two-fold: the development of non-local and exemplar-based approaches for this multimodal interpolation issues 2) .

and the combination of these non-local priors to second-order statistical priors. The later is particularly motivated by the fact that sea surface fields are expected to depict specific fine-scale patterns associated with the advection of the geophysical tracers by turbulent flows. To our knowledge, this is the first work which addresses such combination between patch-based models and statistical priors for multi-resolution interpolation issues.

Reported results for both synthetic and real data demonstrate the relevance of this model compared to the sole use of non-local or statistical priors. Regarding the later, the prior on the spatial covariance of the high-resolution field is of greater importance than the prior on the marginal of the high-resolution details. We have also evaluated that the non-local projection applied to the detail field, rather than to the high-resolution field, offers a greater robustness with respect to possibly weakly represented high-resolution patches (in our case, patches corresponding to cold waters, Fig. 7). The representativeness of the reference patch dataset is obviously a key issue for the effectiveness of the non-local prior. In this respect, we considered reference dataset with 5000 patches. Larger datasets may be considered at the expense of a greater computational cost. One may also explore spatialized version of the patch dataset, where patch similarities account not only for visual similarities but also for space and/or time similarities.

Model parameters for the statistical priors can be fitted from the observed data and/or exploit priors on the considered turbulent dynamics for the spectral prior. The parameters of the non-local prior, especially the size of the patches, were here set empirically. Local patches between 9 × 9 and 13 × 13 were proven relevant. Larger patches typically result in border effects. The automatic learning and tuning of the patch similarity measure could be further explored as developed for other applications of non-local priors in computational imaging [START_REF] Rousseau | A non-local approach for image super-resolution using intermodality priors[END_REF].

We reported a comparison of the proposed interpolation schemes to a kriging interpolation, classically used for January 28, 2016 DRAFT ocean remote sensing data. As kriging-based model states the interpolation as a linear prediction of the missing data points from truly observed ones, it does not address the reconstruction of textural patterns, when the width of the missing data areas is above the spatial range of the considered covariance model (Fig. 5). In operational applications, kriging and optimal interpolation involve a trade-off between the fidelity to the observed data and the interpolation of the missing data areas, with a view to delivering spatially-homogeneous fields. The proposed interpolation models develop a different strategy. We aimed at keeping all the relevant observations, while filling the missing data areas with geophysically-sound field patterns. Interestingly, we showed that the MSE statistics are comparable to a kriging interpolation for small missing data areas. Future work will further explore and evaluate this intrinsic difference with the classical geostatistical tools. In particular, we may consider interpolation ensembles from multiple random initializations. This would provide a mean to evaluate a mean field with persistent (deterministic) structures and some associated variabilities.

From a methodological point of view, our future work will also explore several extensions. Regarding the spectral prior, we only consider here radial (isotropic) covariance structure. The extension to parametric and non-parametric two-dimensional covariance structure could be explored. The same holds for other parametric models of the marginal distribution of the high-resolution detail. The extension to spatio-temporal fields will also be of key interest, including exemplar-based data assimilation issues [START_REF] Tandeo | Combining analog method and ensemble data assimilation: application to the Lorenz-63 chaotic system[END_REF]. Such extension might require manifold learning and kernel learning strategies to define a patch similarity measure adapted to the considered space-time processes [START_REF] Zhao | Analog forecasting with dynamics-adapted kernels[END_REF].

Regarding the application to ocean remote sensing data, we focused here on SST datasets issued from two sensors. The extension of the proposed model to the joint interpolation of three or more data sources should be investigated. For instance, for SST data, operational ODYSSEA product typically aims at merging the daily data issued from at least 6 different sensors [START_REF] Autret | Product User Manual for ODYSSEA Level 3 and 4 global and regional products[END_REF]. Such developments should involve particular methodological issues and could benefit from a classical coarse-to-fine implementation strategy to decompose this multi-sensor/multi-scale fusion as a sequence iteration of two-scale fusions. Applications to other sea surface geophysical fields (e.g., ocean colour, sea surface salinity) appear direct. Beyond the joint processing of multi-sensor data for same geophysical fields (here, SST), one could also explore the relevance of the proposed model for the joint analysis of different geophysical fields for which a mutual conditioning might be expected, but difficult to make explicit through a numerical model. The joint interpolation of SST and ocean colour data may provide a typical example. SST data might provide a mean to constrain the interpolation of ocean colour data, which often involve very large missing data rate in temperate areas.

  Parameters of the spectral prior, Eq.(9)t -σ, νParameters of the prior on the marginal of the details, Eq.(12) λ I , λ M , λ S Gradient descent coefficients for the interpolation, marginal and spectral constraints, Algorithms 2 & 3 0.1 TABLEI

Fig. 1 .

 1 Fig. 1. Interpolation results for a large simulated missing data area: first row, low-resolution SST field (subsampling factor of 16 from a 8-Daubechies wavelet of the ODYSSEA SST field) with superimposed the contour of the missing data area in blue (left), high-resolutionODYSSEA SST field (resolution grid of 2km) (center), interpolated field using the proposed model (right); second row, radial spectral densities of the low-resolution (blue, -), real high-resolution (black,-) and interpolated high-resolution (red, -) fields (left), multifractral spectrum of the low-resolution (blue, -), real (black,-) and interpolated high-resolution fields (center), distribution of high-resolution detail δP s LR W X(2)

Fig. 2 .

 2 Fig. 2. Statistical analysis of interpolation results for the synthetic case considered in Fig.1: first row, radial spectral density (left), multifractal spectrum (center), distribution of high-resolution detail δP s LR W X (2) (right); second row, distribution of the norm of the gradient ∇X

Fig. 3 .X ( 2 )

 32 Fig.3. Gradient fields of the interpolation compared in Fig.2: from left to right and top to bottom, low-resolution field with superimposed the contour of the missing data area, reference high-resolution field and four different interpolation models, namely using only the covariance and marginal constraints (C+M), using only the non-local projection (NL+C+M), using the proposed model, using the proposed model with the non-local projection applied to detail field δP s LR W X(2) rather than field X (2) (NL+C+M+dX).

Fig. 4 .

 4 Fig. 4. Visual comparison to a kriging-based interpolation for the example considered in Fig.1: first row, interpolated fields; second row, associated gradient magnitude fields. We report results using the proposed model (model NL+C+M, left), the proposed model with the non-local projection applied to detail field δP s LR W X (2) (model NL+C+M+dX, center) and a kriging-based model (right).

Fig. 5 .

 5 Fig. 5. Quantitative comparison of the interpolation results reported in Fig.4: we consider the same statistical features as in Fig.2.

Fig. 6 .

 6 Fig. 6. Reconstruction of a daily SST field from a joint analysis of SEVIRI and AMSRE data: first row, AMSR-E field with few missing data (left), SEVIRI field field (center), radial power spectral density of the different fields (right); second row, interpolated AMSR-E field (left), interpolated SEVIRI field (center), scatterplot of AMSR-E and SEVIRI data for the original (black, *) and interpolated (red, •) fields (right).

Fig. 7 .

 7 Fig. 7. Zoom of joint AMSR-E/SEVIRI interpolation results for the situation considered in Fig.6: first row, observed and interpolated AMSR-E fields; second row, observed and interpolated SEVIRI fields; third row, gradient of the AMSR-E fields, fourth row, gradient of the SEVIRI fields.

Fig. 8 .

 8 Fig.8. Joint AMSR-E/SEVIRI interpolation results for a second situation involving large missing data areas for both sensors: first row, observed fields, second row, interpolated fields. For each row, we depict from left to right, the AMSR-E field, the SEVIRI field, the gradient of the AMSR-E field and the gradient of the SEVIRI field. We report a zoom on the same area as in Fig.7off South Africa, which involves the largest missing data rate. The reported interpolation results are issued from the proposed model with a non-local projection applied to detail field δP s LR W X(2) .

Franc

  ¸ois Rousseau is a Professor at Télécom Bretagne / Institut Mines Télécom since 2015. He obtained his PhD in 2003 (University of Rennes) and his Habilitation in 2014 (University of Strasbourg). He worked at the University of California San Francisco in 2004?2005, and then at ICube lab (Université de Strasbourg) as a CNRS research scientist between 2006-2015. His research focuses on image analysis, and more specifically on image reconstruction, segmentation, temporal data and computational anatomy.

  • , 3 • ] and extrapolated to the finer scales. Here, in Algorithm (??), the spectral constraint is applied globally and not only in the missing data area with a view to filtering noisy patterns observed in the daily SEVIRI composite. As patch dataset, we randomly sample joint patches in the observed AMSR-E and SEVIRI fields for the considered data and

	January 28, 2016	DRAFT

By convention, throughout the manuscript, notation P refers to a projection operator in the image domain. Here, notation P W refers to a wavelet-based projection operator. In the subsequent, we introduce other projection operators corresponding to the non-local, covariance-based and margin-based priors.

. is the ceiling function, such that a is the integer part of number a plus one.[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] DRAFT

This comes to apply Algorithm (1) in the Supplementary Material without the non-local projection step.

This comes to apply Algorithm (1) in the Supplementary Material without projections onto the spectral and marginal constraints.[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] DRAFT

[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] DRAFT

ACKNOWLEDGEMENTS

This work was supported by ANR (Agence Nationale de la Recherche, grant ANR-13-MONU-0014). We are grateful to E. Autret, B. Chapron, J.F. Piolle from Ifremer/CERSAT for providing the SST data.

SUPPLEMENTARY MATERIAL

We provide in this annex the detailed iterative algorithm used to solve for constrained minimization [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]). More precisely, we describe in the subsequent:

• Algorithm (1): the algorithm of function projection Interp LR which solves for the projection onto the interpolation and low-resolution constraints (first two constraints in ( 16));'

• Algorithm (2): the algorithm of function projection PSD PDF used to solve for the projection onto the spectral and marginal constraints (last two constraints in ( 16));'

• Algorithm (3): the general iterative algorithm of the proposed approach for the joint interpolation and blending of two multi-sensor fields Algorithm 1 Description of the function which solves for the projection of field X onto interpolation constraint X(Ω) = Y (Ω) and low-resolution constraint P s W (X) = X LR (the first two constraints in [START_REF] Glasner | Super-resolution from a single image[END_REF]. It uses as input parameter high-resolution field X, low-resolution field X LR , missing data region Ω * , scale s, coefficient λ and precision .

function PROJECTION INTERP LR(X,X LR ,Ω,s,λ, )

end function [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] DRAFT Algorithm 2 Description of the function which solves for the projection of field X (1) and X (2) onto marginal and spectral constraints (the last two constraints in in [START_REF] Glasner | Super-resolution from a single image[END_REF]. It uses as inputs fields X (1) and X (1) , missing data regions Ω (1) and Ω (2) , marginal distribution GG σ * ,ν * , radial power spectral density R * , scale s, coefficients λ I , λ M and λ S and precision .

function PROJECTION PDF PSD(X (1) ,X (2) , Ω (1) , Ω (2) ,

X LR ← P s W (X (2) ) Common low-resolution component at scale s, cf. ( 3)

Transform associated with marginal constraint (13), cf. ( 14)

Gradient-based update of X (2) 2) , s, λ I , %% Spectral constraint ( 13)

Gradient-based update of X (2) %%Projection onto interpolation and low-resolution constraints (first two constraints in ( 16))

end while return X (1),n+1 and X (2),n+1

end function January 28, 2016 DRAFT Algorithm 3 Description of the algorithm for the proposed joint interpolation and blending of two fields X (1) and X (2) from constrained minimization [START_REF] Galerne | Random phase textures: Theory and synthesis[END_REF]. It uses as inputs fields initial fields X (1),0 and X (2),0 , missing data regions Ω (1) and Ω (2) , marginal distribution GG σ * ,ν * , radial power spectral density R * , scale s, coefficients λ I , λ M and λ S , precision and iteration number N .

procedure JOINT INTERPOLATION(X (1) ,X (2) , Ω (1) ,Ω (2) , GG σ * ,ν * ,R * ,s,λ I ,λ M ,λ S , ,N)

for n = 0 to N -1 do %% Non-local projection for X (1) followed by the projection onto interpolation and low-resolution constraints, cf.( 16)

,n X (1),n Non-local projection for field X (1) , cf. ( 6)

,n+1 , X LR , Ω (1) , s, λ I , %% Non-local projection for X (2) followed by the projection onto interpolation and low-resolution constraints, cf.( 16)

Non-local projection for field X (2) , cf. ( 6) 2) , s, λ I , %% Projection onto marginal and spectral constraints, cf. ( 16)

,n+1 ← PROJECTION PDF PSD(X (1) , X (2) , Ω (1) , Ω (2) , GG σ * ,ν * , R * , s, λ I , λ M , λ S , ) %% Update of the common low-resolution field and resulting high-resolution fields

Mean low-resolution field X (1),n+1 ← X (1),n+1 -P s LR W X (1),n+1 + X (1,2),n Update of field X (1) given the new low-resolution field X (2),n+1 ← X (2),n+1 -P s LR W X (2),n+1 + X (1,2),n Update of field X (2) given the new low-resolution field GG σ * ,ν * ← ML estimate of the generalized Gaussian distribution of Y (1) -X (1,2),n+1 Ω (1) end for end procedure January 28, 2016 DRAFT