Charlotte Dion

Simone Hermann

Adeline Samson

Mixedsde: an R package to fit mixed stochastic differential equations

Keywords: R, simulation, data analysis, nonparametric estimation, parametric estimation, Bayesian method

Stochastic differential equations (SDEs) are useful to model continuous stochastic processes. When (independent) repeated temporal data are available, variability between the trajectories can be modeled by introducing random effects in the drift of the SDEs. These models are useful to analyse neuronal data, crack length data, pharmacokinetics, financial data, to cite some applications among other. The R package focuses on the estimation of SDEs with linear random effects in the drift. The goal is to estimate the common density of the random effects from repeated discrete observations of the SDE. The package mixedsde proposes three estimation methods: a Bayesian parametric, a frequentist parametric and a frequentist nonparametric method. The three procedures are described as well as the main functions of the package. Illustrations are presented on simulated and real data.

Introduction

Continuous stochastic processes are usually observed discretely in time (with equidistant time points or not) leading to times series, although their intrinsic nature is of continuous time. While discrete time stochastic models such as auto-regressive models (ARMA, GARCH, ...) have been widely developed for time series with equidistant times, more and more attention have been focused on Stochastic Differential Equations (SDEs). Examples of applications where SDEs have been used include dynamics of thermal systems [START_REF] Bacher | Identifying suitable models for the heat dynamics of buildings[END_REF], solar and wind power forecasting [START_REF] Iversen | Probabilistic forecasts of solar irradiance using stochastic differential equations[END_REF], neuronal dynamics [START_REF] Ditlevsen | Estimation in partially observed stochastic Morris-Lecar neuronal model with particle filter and stochastic approximation methods[END_REF], pharmacokinetic/pharmacodynamic (PK/PD) [START_REF] Hansen | Predicting plasma glucose from interstitial glucose observations using bayesian methods[END_REF], crack growth [START_REF] Hermann | Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model[END_REF].

Depending on the applications, independent repeated temporal measures might be available. For examples, drug concentration of several subjects is usually observed in PK; dynamics of several neurons is measured along time; time to crack lengths can be measured repeatedly in crack growth study. Each trajectory represents the behaviour of a unit/subject. The functional form is similar for all the trajectories. Fitting the overall data simultaneously obviously improves the quality of estimation, but one has to take into account these variabilities between experiments. This is the typical framework of mixed-effects models where some parameters are considered as random variables (random effects) and proper to each trajectory. Hence the random effects represent the particularity of each subject. Some parameters can also be considered as common to all the trajectories (fixed effects).

1

In this work the model of interest is thus a mixed-effects stochastic differential equation (MSDE), mixed-effects for both fixed and random effects. The mixedsde package has been developed to estimate the density of the random effects from the discrete observations of M independent trajectories of a MSDE. It is available from https://github.com/charlottedion/mixedsde.

More precisely, we focus on MSDE with linear drift. We consider M diffusion processes (X j (t), t ≥ 0), j = 1, . . . , M with dynamics ruled by SDE, for t ∈ [0, T ] dX j (t) = (α j -β j X j (t))dt + σa(X j (t))dW j (t) X j (0) = x j (1) where (W j ) 1...j...M are M independent Wiener processes, (α j , β j ) are two (random) parameters, σa(X j (•)) is the diffusion coefficient with a a known function and σ an unknown constant. The initial condition x j is assumed fixed (and known) in the paper with possibly different values for each trajectory.

In the package, we restrict the models to the two famous SDEs with linear drift, namely the Ornstein-Uhlenbeck model (OU) with a(x) = 1 and the Cox-Ingersoll-Ross model (CIR) with a(x) = √ x. For the CIR model, we assume that x j > 0, σ > 0, α j > σ 2 /2 and β j > 0 to ensure that the process never crosses zero.

The random parameters are denoted φ j and belong to R d with either d = 1 or d = 2:

• (d = 1) φ j = α j random and for all j = 1, . . . , M , β j = β fixed

• (d = 1) φ j = β j random and for all j = 1, . . . , M , α j = α fixed

• (d = 2) φ j = (α j , β j ) random
The φ j 's are assumed independent and identically distributed (i.i.d.) and independent of the W j 's. The mixedsde package aims at estimating the random effects φ j and their distribution whose density is denoted f , from N discrete observations of the M trajectories (X j (t)) j from equation ( 1) at discrete times t 0 = 0 < t 1 < . . . < t N = T (not necessarily equidistant). To the best of our knowledge, this is the first package in R language dedicated to the estimation of MSDE. Estimation procedures for MSDE have been proposed in the non-parametric and the parametric frameworks, with a frequentist and a Bayesian point of view. The parametric approaches assume Gaussian random effects φ j . Among other references, for parametric maximum likelihood estimation, we can cite [START_REF] Ditlevsen | Mixed Effects in Stochastic Differential Equation[END_REF] Three estimation procedures are implemented in the mixedsde package: a kernel nonparametric estimator [START_REF] Dion | Bidimensional random effect estimation in mixed stochastic differential model[END_REF], a parametric maximum likelihood estimator [START_REF] Delattre | Maximum likelihood estimation for stochastic differential equations with random effects[END_REF] and a parametric Bayesian estimator [START_REF] Hermann | Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model[END_REF]. The parametric frequentist and Bayesian approaches assume the random effects Gaussian. The Bayesian approach seems the most appropriate method for a small time of observation T and a small number of trajectories M . The nonparametric approach can be used when no prior idea on the density is available and when T and M are both large enough. Finally, the parametric frequentist estimation can be used when data are symmetric with a large number of discrete observations. This paper reviews in Section 2 the three estimation methods. An overview of the mixedsde package is given in Section 3 through a description of the main functions and of other related companion functions. The practical use of this package is illustrated in Section 4 on simulated data and in Section 5 on one real dataset in neuronal modelling.

Density estimation in mixed stochastic differential models

We briefly recall the methodology of the three estimators implemented in the mixedsde package. We start with the nonparametric approach, then the frequentist parametric Gaussian method and finally the Bayesian parametric Gaussian method.

Nonparametric estimation of the random effects density

The first step of the nonparametric approach is to estimate the random effects. The idea is to maximize the likelihood of the process X ϕ j solution of the stochastic differential equation with fixed ϕ. Assuming continuous observations of (X j (t), 0 ≤ t ≤ T ), the likelihood function is obtained with the Girsanov formula:

T (ϕ) = exp T 0 α -βX ϕ j (s) σ 2 a 2 (X ϕ j (s)) dX j (s) - 1 2 T 0 (α -βX ϕ j (s)) 2 σ 2 a 2 (X ϕ j (s))
ds .

Maximizing the likelihood yields to the following estimator of φ j

A j := V -1 j U j (2)
where U j and V j are the two sufficient statistics of the model. They are explicit depending on the form of the random effects:

• α j random and β known

U j := T 0 1 σ 2 a 2 (X j (s)) dX j (s) + β T 0 X j (s) σ 2 a 2 (X j (s)) ds, V j := T 0 1 σ 2 a 2 (X j (s)) ds
• β j random and α known

U j := - T 0 X j (s) σ 2 a 2 (X j (s)) dX j (s) + α T 0 X j (s) σ 2 a 2 (X j (s)) ds, V j := T 0 X j (s) 2 σ 2 a 2 (X j (s)) ds
• (α j , β j ) random, denote b(x) = (1, -x) t with u t the transposition of vector u. Here U j is a column vector with size 2 × 1 and V j = (V j,k, ) k, ∈{1,2} a 2 × 2 symmetric matrix:

U j := T 0 b σ 2 (X j (s))dX j (s), V j := T 0 b b t σ 2 (X j (s))ds. (3)
Truncated versions of this estimator have been introduced for theoretical reasons. In the bidimensional case φ j = (α j , β j ), [START_REF] Dion | Bidimensional random effect estimation in mixed stochastic differential model[END_REF] propose the following estimator

A j := A j 1 Bj , B j := {V j ≥ κ √ T I 2 } = {min(λ 1,j , λ 2,j ) ≥ κ √ T } (4)
with I 2 the 2×2 identity matrix and λ i,j , i = 1, 2 the two eigenvalues of the symmetric non negative matrix V j , and κ a numerical constant calibrated in practice. In the one-dimensional case φ j = β j with α = 0, Genon-Catalot and Larédo (2016) propose

A j := A j 1 Vj ≥κ √ T (5)
with κ a numerical constant calibrated in practice.

Based on these estimators of the φ j 's, we can proceed to the second step, the estimation of their density f . Several nonparametric estimators of f have been proposed (see Comte et al., 2013, for example). In the package mixedsde, we focus on the kernel estimator of f . Let us introduce the kernel function K : R d → R, with d = 1, 2 depending on the dimension of φ j . We assume K to be a C 2 function satisfying

K(u)du = 1, K 2 = K 2 (u)du < +∞, (∇K(u)) 2 du < +∞ (with ∇K the gradient of K). A bandwidth h ∈ (R + ) d , for d = 1, 2, is used to define the function K h (x) = 1 h K x h , x ∈ R d .
Note that in the bidimensional case, h = (h 1 , h 2 ) and the two marginal bandwidths are different.

The nonparametric estimator of the density f of φ j is

f h (x) = 1 M M j=1 K h (x -A j ). (6)
and the estimator

f h (x) = 1 M M j=1 K h (x -A j ) is computed when the truncated estimator Âj is different than A j .
In the mixedsde package, Gaussian kernel estimators are implemented with the R-functions density (available in package stats) when d = 1 and kde2d (available in package MASS) when d = 2 with an automatic selection of the bandwidth h. Note that when there is only one random effect, the bandwidth is selected by cross-validation with the argument bw="ucv", or as the default value given by the rule-of-thumb if the chosen bandwidth is to small. Note that the estimator is unstable for small variance of the random effects.

It is important to notice that the two random effects are not assumed independent. When there is only one random effect, the fixed parameter has to be entered by the user.

The computation of A j = V -1 j U j does not require the knowledge of σ 2 as it appears both in U j and V j . It requires however the evaluation of the two continuous integrals in U j and V j while observing the trajectories (X j ) at discrete times (t 0 , t 1 , . . . , t N ). For ∆ k = t k+1 -t k , k = 0, . . . , N -1, the non-stochastic integrals T 0 g(X j (s))ds for any function g are approximated by

T 0 g(X j (s))ds ≈ N -1 k=0 g(X j (t k ))∆ k .
For the stochastic integrals, we use the following simple discretization

T 0 g(X j (s))dX j (s) ≈ N -1 k=0 g(X j (t k ))(X j (t k+1 ) -(X j (t k )))∆ k .

Frequentist parametric estimation approach

In this section and the following one, we assume that the random parameters φ j are Gaussian:

• when d = 1, φ j ∼ N (µ, ω 2 ) with µ ∈ R • when d = 2, φ j ∼ N (µ, Ω) with µ ∈ R 2 and a diagonal covariance matrix Ω = diag(ω 2 1 , ω 2
2 ). For the bidimensional case d = 2 we estimate by maximum likelihood the parameters θ := (µ, Ω). We define the likelihood function assuming first that the trajectories are continuously observed, similarly to the nonparametric approach (Section 2.1). Thanks to the Girsanov formula, the likelihood function of the j th trajectory X j is

L(X j , θ) = 1 det(I 2 + ΩV j ) exp - 1 2 (µ -V -1 j U j ) R -1 j (µ -V -1 j U j ) exp 1 2 U j V -1 j U j with R -1 j = (I 2 + V j Ω) -1
V j and I 2 is the 2 × 2 identity matrix. For the case d = 1, the parameters to estimate are θ := (µ, ω, ψ) where ψ denotes the fixed effect α or β. We adopt the subscript r for the value of random, equal to 1 or 2, and c for the position of the common fixed effect (thus 2 or 1). The likelihood function of the j th trajectory X j is

L(X j , θ) = 1 1 + ω 2 V j,r,r exp - 1 2 V j,r,r (1 + ω 2 V j,r,r ) -1 (µ -V -1 j,r,r (U j,r -ψV j,c,r )) 2 × exp ψU j,c - ψ 2 2 V j,c,c exp 1 2 (U j,r -ψV j,r,c ) 2 V -1 j,r,r
with the notations U , V from (3). Details on this formula are available in the Appendix 6.

The likelihood function is defined as L(θ) = M j=1 L(X j , θ). The maximum likelihood estimator θ := ( µ, Ω, ψ) when d = 1 and θ := ( µ, Ω) when d = 2 is defined by

θ = arg max θ L(θ) = arg max θ M j=1 L(X j , θ)
This estimator is not explicit. In the mixedsde package, the function optim is used to maximize numerically the likelihood. Initialization of optim is obtained by computing the mean and the variance of the estimators A j of the random parameters (see equation ( 2)). Sufficient statistics U j and V j are discretized as explained in Section 2.1.

Note that this parametric approach requires the knowledge of σ 2 to compute the sufficient statistics U j and V j because V j appears alone in R j . We plug the following estimator of σ 2

σ 2 = 1 M M j=1 1 N N -1 k=0 ((X j (t k+1 ) -X j (t k )) 2 ∆ k a 2 (X j (t k )) . (7)
Selection of (non-nested) models can be performed with the BIC criteria, defined by -2 log L( θ)+ 2 log(M ) for model with one random effect and -2 log L( θ) + 4 log(M ) with two random effects and the AIC criteria defined by -2 log L( θ) + 2 for one random effect and -2 log L( θ) + 4 for two random effects. These asymptotic criteria indicate the trade-off between maximizing fit and minimizing model complexity. Note that their theoretical properties are guaranteed only when σ 2 is known.

Bayesian parametric approach

For the Bayesian approach we assume similarly to the frequentist parametric estimation method a Gaussian distribution for φ j , with a diagonal covariance matrix Ω = diag(ω 2 1 , ω 2 2 ). In this method, we estimate in the same time the diffusion coefficient σ. The parameters of interest are thus θ = (µ, Ω, σ) and we want to estimate their posterior distribution p(θ|(X j (t k )) j=1,...,M,k=1,...,N ). Let denote X 1:M = (X j (t k )) j=1,...,M,k=1,...,N in the following.

We now introduce prior distributions implemented in mixedsde package for the parameters θ:

µ ∼ N (m, V ), V = diag(v) ω 2 i ∼ IG(α ω,i , β ω,i ), i = 1, 2 σ 2 ∼ IG(α σ , β σ ),
where IG is the inverse Gamma distribution which is conjugate to the normal likelihood and m, V, α ω,i , β ω,i , α σ , β σ are hyperparameters fixed by the user. The case of only one random effect is nested by setting ω 2 1 or ω 2 2 equal to zero. The aim is to calculate the posterior distribution p(θ|X 1:M ) which is not explicit for the whole vector of parameters. Therefore, we simulate it through a Gibbs sampler (see e.g. [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]. Here, we have a true transition density of both processes that is used for the likelihood, see [START_REF] Iacus | Simulation and Inference for Stochastic Differential Equations[END_REF]. For a general hierarchical diffusion approach based on the Euler approximation, see [START_REF] Hermann | Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model[END_REF].

Analogically to the frequentist approach, there is a first step: sample from the full conditional posterior of the random effects p(φ j |(X j (t k )) k=1,...,N , θ), j = 1, . . . , M . This is done by a Metropolis Hastings (MH) algorithm.

The second step is the estimation of the hierarchical parameters µ and Ω. Full conditional posteriors p(µ|φ 1 , . . . , φ M , Ω) (resp. p(Ω|φ 1 , . . . , φ M , µ)) are Gaussian (resp. inverse Gamma) and can, for example, be found in [START_REF] Hermann | Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model[END_REF].

The last step of the Gibbs sampler is sampling from the full conditional posterior of σ 2 . For the CIR model, this is also conducted by a MH step. For the OU model, the inverse Gamma distribution is conjugate to the normal likelihood. The full conditional posterior distribution is given by

σ 2 |X 1:M , φ 1 , ..., φ M ∼ IG   α σ + M N 2 , β σ + 1 2 M j=1 N k=1 β j 1 -e -2βj ∆ k X j (t k ) - α j β j -X j (t k-1 ) - α j β j e -βj ∆ k 2   .
In the case of one random effect, there is one additional Gibbs sampler step for the fixed effect, that is also conducted through a MH algorithm.

In the package, the starting values for the Gibbs sampler are set equal to the mean of the prior distributions. In all the MH algorithms, one each has to choose a proposal density. In the package mixedsde, we use a normal density for all location parameters with mean equal to the last chain iteration and a proposal variance that has to be chosen. For the CIR model, the proposal distribution for σ 2 is chosen by

√ σ 2 ∼ N ( σ 2 prev , variance)
where σ 2 prev is the previous value of σ 2 . The remaining question is how to choose the suitable proposal variance. This variance controls the chain dependence and the acceptance rate. If the variance is small, the acceptance rate is large and the chains gets very dependent. If the proposal variance is large, only few candidates are accepted with the advantage of weakly dependent chains. This problem is solved in the package with an adaptive Metropolis-within Gibbs algorithm [START_REF] Rosenthal | Optimal proposal distributions and adaptive MCMC[END_REF] using the proposal distribution N (0, e 2l ) with l the logarithm of the standard deviation of the increment. This parameter is chosen so that the acceptance rate is approximately 0.44 which is proposed to be optimal in the Metropolis-within Gibbs sampler. It is proposed to add/subtract an adoption amount δ(n) = min(0.1, n -1/2 ) to/from t after every 50th iteration and adapt the proposal variance if the acceptance rate is smaller than 0.3 or larger than 0.6.

Predictions

In many cases, one is not only interested in parameter estimation but also in the prediction for future observations. The first step is the prediction of a future random effect φ pred . The simulation of a new random effect is direct for the frequentist parametric approach sampling from N ( µ, Ω). For the nonparametric approach, first note that f h is an estimator given on a discrete grid {x 1 , . . . , x n }, i.e. a vector of corresponding {p 1 , . . . , p n } after normalisation. Simulating from the estimator f h can therefore be performed simulating a discrete variable from vector {x 1 , . . . , x n } with (normalized) probabilities {p 1 , . . . , p n }. For the Bayesian approach, a new φ pred is sampled from the predictive distribution p(φ pred |X 1:M ) = p(φ pred |µ, Ω)p(µ, Ω|X 1:M ) d(µ, Ω) where the posterior of µ and Ω is approximated by the results of the Gibbs sampler. This distribution is not explicit, we propose to sample over a grid through inversion method, equal to the nonparametric case.

Given a new random effect φ pred , we are able to simulate predictive trajectories. This is performed using the transition density p(X(t k )|X(t k-1 ), φ pred , σ 2 ) for the frequentist approach. The starting points of the process x j are the observed ones. For the Bayesian approach, we implement two prediction settings. Firstly, analogously to the frequentist approach a new trajectory is simulated using the transition density p(X(t k )|X(t k-1 ), φ pred , σ 2 ) where φ pred is sampled from the MCMC posterior distribution p(φ|X 1:M ). Secondly, we can calculate the predictive distribution

p(X(t i )|X 1:M ) = p(X(t i )|φ pred , σ 2 )p(φ pred , σ 2 |X 1:M ) d(φ pred , σ 2 )
in each time point. We can then calculate only the quantiles for a prediction interval or to draw directly samples from the predictive distribution. For this predictive distribution, we take the starting point x j = x 0 to be the same for all series. If the starting points would vary, this is an additional random effect whose density has to be estimated. This is not implemented in the estimation procedure and will, therefore, left out for the prediction.

It is then interesting to compare the new trajectories to the real ones and if the number of new trajectories is large enough, to compute an empirical confidence interval.

Overview of the mixedsde functions

This Section presents an overview of the functions implemented in the package. Illustrations of the code are given in Section 4.

Data

Data is a matrix X of size M ×N for M trajectories with N time points, not necessarily equidistant but similar for the M trajectories. The vector of times is a vector times of length N . Real datasets are available on the package, and detailed on Section 5.

To lead a simulation study, the function mixedsde.sim allows to generate a list with a M × N matrix X of M trajectories on the interval [0, T ] with N equidistant points (default value 100) and a vector times with the equidistant times. This function leans on function sde.sim available via package sde from the web address http://cran.r-project.org/package=sde, to simulate SDE. One has to choose: model either OU or CIR; random that fixes the position and the number of random effects: random = 1 for α j random, random = 2 for β j random or random = c(1,2) for α j and β j random; σ the diffusion coefficient; invariant, default value 0 means that X 0 is 0 (default) or fixed by the user, value 1 means that X 0 is generated from the invariant distribution (see details in the package documentation); density.phi to choose the distribution of the random effect (see package documentations).

Main function

Main function is mixedsde.fit producing estimation of the random effects and their common density. Inputs of mixedsde.fit are • X a M × N matrix containing the trajectories by rows

• times the vector of observations times

• model the chosen model either OU or CIR

• random that fixes the position and the number of random effects: random = 1 for α j random, random = 2 for β j random or random = c(1,2) for α j and β j random

• estim.method the estimation method: nonparam (see Section 2.1), paramML (see Section 2.2) or paramBayes (see Section 2.3)

• fixed the value of the fixed effect β (resp. α) when random = 1 (resp. random = 2), default 0. (Only for the frequentist approaches).

• estim.fix 1 if the fixed effect is estimated, default 0. (Only for the frequentist parametric approach when random=1 or 2). • gridf the x-axis grid on which the random effect distribution is computed: we recommend a fine grid with at least 200 points, default value is a sequence of length 500 starting in 0.8 × min j φ j and ending in 1.2 × max j φ j . (Only for the frequentist approaches)

• prior the list of prior parameters m, v, alpha.omega, beta.omega, alpha.sigma, beta.sigma for paramBayes method: default values are calculated based on the estimations (A j ) j for the first min(3, M • 0.1 ) series and main estimation is only made with the remaining M • 0.9 . (Only for the Bayesian approach)

• nMCMC length of the Markov chain for paramBayes method. (Only for the Bayesian approach)

In the following we describe the related methods, proposed in the package, they are summarised in Table 1.

Outputs

Output of mixedsde.fit is a S4 class called Freq.fit for the frequentist approaches and Bayes.fit for the Bayesian approach. Results of the estimation procedure are available as a list applying function out to the Freq.fit (resp. Bayes.fit) object. Elements of Freq.fit are:

• sigma2 estimator σ 2 of the diffusion coefficient

• estimphi estimator (A j ) j of the random effects • estimphi.trunc truncated estimation ( A j ) j of the random effects • estim.fixed estimator of the fixed effect if random = 1 or 2, estim.method = paramML; estim.fix = 1, default 0 • gridf the x-axis grid on which the random effect distribution is computed

• estimf estimator of the density of the random effects (for both paramML and nonparam methods)

• cutoff binary M -vector of binary values indicated the truncated trajectories, default FALSE when no truncation

• estimf.trunc truncated estimation of the density of the random effects

• mu estimation of Gaussian mean of the random effects (only for paramML method)

• omega estimation of Gaussian variance matrix of the random effects (only for paramML method)

• aic and bic AIC and BIC criteria (only for paramML method)

• index index of trajectories used for the estimation, excluded are trajectories with V j = 0 or V j = +∞ (one random effect) or det V = +∞ (two random effects), trajectories containing negative values for CIR model, Elements of Bayes.fit are:

• sigma2 trace of the Markov chain simulated from the posterior of σ 2

• mu trace of the Markov chain simulated from the posterior of µ

• omega trace of the Markov chain simulated from the posterior of ω 2

• alpha trace of the Markov chain simulated from the posterior of α j , nMCMC×M matrix if α is random effect, nMCMC×1 otherwise

• beta trace of the Markov chain simulated from the posterior of β j , nMCMC×M matrix if β is random effect, nMCMC×1 otherwise

• burnIn a proposal for the burn-in phase

• thinning a proposal for the thin rate

• ind.4.prior the indices used for the prior parameter calculation, M + 1 if prior parameters were specified.

Outputs burnIn and thinning are only proposals for a burn-in phase and a thin rate. The proposed burnIn is calculated by dividing the Markov chains into 10 blocks and calculate the 95% credibility intervals and the respective mean. Starting in the first one, the block is taken as burn-in as long as the mean of the current block is not in the credibility interval of the following block or vice versa. The thinning rate is proposed by the first lag which leads to a chain autocorrelation of less than 80%. It is not easy to automate these choices, so it is highly recommended by the authors to verify the chains manually.

Command plot() applied to a Freq.fit object produces a frequencies histogram of (A j (T )) j (one or two according to the number of random effects) with the estimated density (red curve) and the truncated estimator if available (dotted grey red curve) and a quantile-quantile graph with the quantiles of the A j 's versus the quantiles of a normal sample of the same length, with the same empirical mean and standard deviation. This illustrates the normality of the sample. Applying this function to the nonparametric results indicates if the Gaussian assumption of the parametric approach is appropriate. When plot() is applied to a Bayes.fit object, one can choose four different options, named style. The default value is chains, it plots the Markov chains for the different parameter values. acf leads to the corresponding autocorrelation functions, density to the approximated densities for each parameter and cred.int leads to the credibility intervals of the random parameters with the input parameter level with default 0.05. For all options, with the input parameter reduced = TRUE, the burn-in period is excluded and a thinning rate is taken, default is FALSE. There is also a possibility to include the prior means in the plots by lines with plot.priorMean = TRUE, default is FALSE.

In the Bayesian estimation the influence of prior parameters is interesting, thus for the Bayes.fit object, there is a second plot method, named plot2compare where three estimation objects can be compared. For reasons of clarity, only the densities are compared, with the default reduced = TRUE. Here, there is also a possibility to include true.values, a list of the true parameters for the comparison in a simulation example.

Command summary() applied to a Freq.fit object computes the kurtosis and the skewness of the distribution, σ 2 , the empirical mean and standard deviation computed from the estimators (A j ) j , µ, Ω (and the fixed effect α or β), AIC, BIC criteria for the frequentist MLE method. When applied to a Bayes.fit object, it computes means and credibility interval (default level 95%) for each parameter (µ, Ω, σ, α, β). Here, there is also a possibility to choose the burn-in and the thinning rate manually by the input parameters burnIn and thinning.

Command print() applied to a Freq.fit object returns the use or not of the cutoff and the vector of excluded trajectories. When applied to a Bayes.fit object, it returns the acceptance rates of the MCMC procedure.

Validation methods

Validation of a mixed model, obtained with function valid, is an individual validation. Indeed, the validation of trajectory j is obtained comparing it to M new trajectories simulated with parameters (α, β) fixed to the estimator A j (or A j ) in the frequentist approaches and to the posterior means in the Bayesian approach. Inputs of the function are • Freq.fit or Bayes.fit object • plot.valid 1 to generate a figure (default value is 1)

• numj a specific individual trajectory to validate (default all trajectories)

• Mrep number of simulated trajectories (default value 100)

Each observation X numj (t k ) is compared with the Mrep simulated values (X 1 numj (t k ), . . . , X Mrep numj (t k )), for k = 1, . . . , N .
Outputs are the list of the (X

(t k )), for k = 1, . . . , N .
This is an empirical method. The recent work [START_REF] Kuelbs | Limit Theorems for Quantile and Depth Regions for Stochastic Processes High dimensional probability VII-Progress in Probabiliy[END_REF] on depth and quantile regions for stochastic processes (see for example [START_REF] Zuo | General notions of statistical depth function[END_REF] for depth functions definitions) should provide the theoretical context for a more extensive study. This could be done in further works. (and optional plot parameters). Function pred applied to a Freq.fit object returns a list with predicted random effects phipred, predicted trajectories Xpred and indexes of the corresponding true trajectories indexpred (see Section 2.4 for details of simulation). If plot.pred = TRUE (default) three plots are produced. Left predicted random effects versus estimated random effects. Middle: true trajectories. Right predicted trajectories and their empirical 95% prediction intervals (default value level=0.05). The prediction can also be done from the truncated estimator f h based on the A j (5), if the argument pred.trunc = 1.

Prediction methods

Function pred applied to a Bayes.fit object returns a S4 class object Bayes.pred. The first element of this class is Xpred, which depends on the input parameters. Including the input trajectories = TRUE, matrix Xpred contains the M drawn trajectories by rows (see first method described for the Bayesian approach in Section 2.4). Default is trajectories = FALSE which leads to the calculation of the predictive distribution explained in Section 2.4. With the input only.interval = TRUE (default), only the quantiles for the 1level prediction interval are calculated, stored in qu.l and qu.u. Input only.interval = FALSE provides additionally Xpred containing sample.length (default 500) samples from the predictive distribution in each time point of the observations (except the first). In both cases, with plot.pred = TRUE, two figures are produced. On the left side, the data trajectories are compared with the prediction intervals and on the right side, the coverage rate is depicted which is stored in entry coverage.rate, namely the amount of series covered by the prediction intervals for each time point. The last class entry estim stores the results from the Bayes.fit object in a list. Other input parameters are burnIn and thinning which allow for the choice of other burn-in phase and thinning rate than proposed in the Bayes.fit object.

For the Bayes.pred class object, two plot methods are available. plot() repeats the figures that are created with the plot.pred = TRUE command in the pred method. plot2compare() compares up to three Bayes.pred objects, where in a first figure the prediction intervals are presented in colors black, red and green and the observed data series in grey and in a second figure the corresponding coverage rates are compared. With the input parameter names a vector of characters to be written in a legend can be indicated.

Package mixedsde through simulated examples

In this part two simulated examples are given to illustrate the strengths of each proposed method. Two datasets are simulated according to:

1. CIR model with one non-Gaussian random effect β j ∼ Γ(1.8, 0.8), α j = 1, T = 50, M = 200, N = 1000: model1 <-"CIR"; random1 <-2; fixed1 <-1; sigma1 <-0.1 ; M1 <-200; T1 <-50; N1 <-1000; X01 <-1; density.phi1 <-"gamma"; param1 <-c(1.8, 0.8); simu1 <-mixedsde.sim(M=M1, T=T1, N=N1, model=model1, random=random1, fixed=fixed1, density.phi=density.phi1, param=param1, sigma=sigma1, X0=X01) X1<-simu1$X; phi1 <-simu1$phi; times1 <-simu1$times 2. OU model with one Gaussian random effect α j ∼ N (3, 0.5 2 ), β j = 5, T = 1, M = 50, N = 500: model2 <-"OU"; random2 <-1; sigma2 <-0.1; fixed2 <-5; M2 <-50; T2 <-1; N2 <-500; X02 <-0; density.phi2 <-"normal"; param2 <-c(3, 0.5); simu2 <-mixedsde.sim(M=M2, T=T2, N=N2, model=model2, random=random2, fixed=fixed2, density.phi=density.phi2, param=param2, sigma=sigma2, X0=X02) X2 <-simu2$X; phi2 <-simu2$phi; times2 <-simu2$times Example 1 has non Gaussian random effect, the nonparametric method is the most appropriate approach. Example 2 has T small and Gaussian random effect, nonparametric method is therefore not the most appropriate approach. Parametric methods should performed well with a preference to the Bayesian as the number of trajectories M 2 = 50 is not large. Nonparametric estimation fits well the histogram of (A j ) (left plot) and we see that the random effects are non-Gaussian (right plot). Because we are working on simulated data, we can compare the estimations with the true random effects and the true f :

# comparison of the true f and its estimation gridf1 <-outputsNP$gridf f1 <-dgamma(gridf1, param1 [1], param1 [2]) # true density function fhat <-outputsNP$estimf # nonparametric estimated density function plot(gridf1, f1, type='l', lwd=2, xlab='', ylab='')

# plot of the true density lines(gridf1, fhat, col='red') # plot of the nonparametric estimated density # comparison of the true random effects and their estimations phihat1 <-outputsNP$estimphi # estimated random effects plot(phi1, phihat1, type = "p", pch = 18, xlab='', ylab='') # true vs estimated abline(0, 1) It yields to Figure 2. On the left plot, the estimated density (dotted curve) is very close to the true density f (plain line). The right plot shows that A j is a good estimation of φ j . This confirms that the nonparametric approach performs well for this settings. Validation of the MSDE is produced by function valid. The two graphs on the right of The small number of observations makes the frequentist estimation harder, nevertheless here, the histogram seems pretty well fitted by the parametrically estimated density. Because we are working on simulated data, we can compare the estimations with the true random effects and the true f : # comparison of the true f and its estimation gridf2 <-outputsP$gridf f2 <-dnorm(gridf, param [1], param [2]) # true density fhat_param <-outputsP$estimf # parametric estimated density plot(gridf2, f2, type = 'l', lwd = 2, xlab = '', ylab = '') # plot of the true density lines(gridf2, fhat_param, col='red', lty = 2, lwd = 2) # plot of the estimated density # comparison of the true random effects and their estimations phihat2 <-outputsP$estimphi # estimated random effects plot(phi2, phihat2, type="p", pch=18, xlab='', ylab='') # true vs estimated abline(0, 1) It yields to Figure 4. It shows that estimation of the density is satisfactory (left) and estimation of the random effects is very good (right).

Validation of the MSDE is produced by function valid. For example the individual validation of the first trajectory is plotted Figure 5, the first two graphs on the left, using This illustrates the good estimation of the random effects: a beam of trajectories with the true one in the middle and the lining up of the quantiles.

Finally, we can predict some trajectories using pred. Predictions are shown on Figure 6, as a result of predPOU <-pred(estim_param) Beam of predicted trajectories (right) is close to the true one (middle). The lining up of the predicted random effects versus the estimated random effects (left) shows the goodness of the prediction from the estimated density, thus of the estimation of the density.

Bayesian estimation

Bayesian method is applied to Example 2. Priors are constructed from the true values, but default values can be used.

prior2 <-list( m = c(param2 [1], fixed2), v = c(param2 [1], fixed2), alpha.omega = 11, beta.omega = param2 [2]^2*10, alpha.sigma=10, beta.sigma = sigma^2*9) estim.method <-'paramBayes' estim_bayes <-mixedsde.fit(times2, X2, model = 'OU', random = 2, estim.method, prior = prior2, nMCMC = 10000) print(estim_bayes) [1] "acceptance rates for random effect:" Min. 1st Qu. Median Mean 3rd Qu. Max. 0.5607 0.5675 0.5694 0.5705 0.5730 0.5828 [1] "acceptance rate for fixed effect:" "0.3591"

The fixed effect β has a small acceptance rate, explaining the dependent chain (Figure 7 top right). This is due to a very sharp likelihood because of the large amount of observations (N • M ) in comparison to the random effect (N ).

Predictions in the Bayesian framework and the corresponding Figure 8 is obtained by pred.result <-pred(estim_bayes)

Figure 8 shows the beam of simulated data trajectories together with the 95% prediction interval. Coverage rates are shown on the right plot and we see that the intervals hold the level.

5 Package mixedsde through a real data example A real dataset is available (neuronal.data.rda) through lists of a matrix X and a vector times. We detail below the analysis of this dataset, following the next steps: run the two random effects model with both the parametric and nonparametric procedure; choose the number of random effects depending on the variability of the estimators (A j,1 , A j,2 ), on the shape of f h and the variance Ω.

Neuronal data (neuronal.data)

Neurons are the basement of nervous system and each neuron is connected with around 1000 other neurons. They are communicating through emission of electrical signal. We focus on the dynamic of the neuron membrane potential between two spikes emission measured in volts as the difference of ions concentration between the exterior and the interior of the cell. Data are obtained from one single neuron of a pig. Data are composed of M = 240 membrane potential trajectories with N = 2000 equidistant observation times. Time step is δ = 0.00015 [s] and observation time is T = 0.3 [s]. Data are uploaded using data(neuronal.data). They are presented on Figure 9.

These data have been previously analysed with a Ornstein-Uhlenbeck model with one additive random effect (α j ): [START_REF] Picchini | Parameters of the diffusion leaky integrate-and-fire neuronal model for a slowly fluctuating signal[END_REF] and [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF] use parametric methods assuming the normality of the random effect, and Dion (2014) with a nonparametric method. Here α j represents the local average input that the neuron receives after the j th spike. The initial voltage (the value following a spike) is assumed to be equal to the resting potential and set to zero: x j = 0. Parameter β j (non negative) is the time constant of the neuron. It was fixed in [START_REF] Picchini | Parameters of the diffusion leaky integrate-and-fire neuronal model for a slowly fluctuating signal[END_REF] and [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF].

In this new analysis, we assume that both α j and β j may change from one trajectory to another because of other neurons or environment influence, for example. Indeed, the form of each trajectory lead us to think that this is the good model: each one has its mean-revert value and its own speed to reach this value. There is no reason to assume that the speed is always the same, but looking at the trajectories the stationary state seem reached nearly at the same time, thus the second random effect should have a small variance.

Fitting with MSDEs

Our goal is also to compare the two models OU and CIR, both with two random effects, and two approaches: the nonparametric density estimation and the parametric density estimation. Let us remark that for the CIR model the algorithm removes two trajectories: 168 and 224, because they contain negatives values. For two random effects the command is estim <-mixedsde.fit(times, X=X, model= model, random = random, estim.method = estim.method)

and they can be found in the help data file (command ?neuronal.data). We first apply the two frequentist approaches on models with two random effects. Kurtosis and skewness of the distribution of the estimation A j of the random effects given in Table 2 are not closed to a symmetric distribution. The bidimensional density of (α j , β j ) is estimated for both models with the parametric and nonparametric methods running function mixedsde.fit. Figure 10 gives the 4 estimated marginals. The blue (black) is for the OU model and the green (grey) for the CIR model. The dotted lines are the estimations from the parametric method, the plain lines for the nonparametric estimation. Parametric and nonparametric estimators are close, except for the second random effect with the OU model. Indeed, parametric estimation produces a very small variance for the second random effect, suggesting it could be fixed. Would this assumption be valid, it explains the difference with the nonparametric estimator which is not stable if the variance is to small. Estimation of σ is σ = 0.0136 for the OU model and σ = 0.163 for the CIR model.

To compare with previous literature results, we focus on the OU model. To select the number and the position of the random effect, we run the code with one random effect, additive or multiplicative: random=1 or random = 2, for both models estimating the common fixed parameter. Estimators of the means µ 1 , µ 2 and standard deviations ω 1 , ω 2 are given in Table 3 and compared to values obtained in [START_REF] Picchini | Parameters of the diffusion leaky integrate-and-fire neuronal model for a slowly fluctuating signal[END_REF] and [START_REF] Picchini | Stochastic differential mixed-effects models[END_REF]. Criteria AIC and BIC are also given in Table 3. The preferred model is the one minimizing both criteria. Thus, the OU model with one additive random effect φ j = α j and β = 37 seems to be the best model to describe these data. The summary method gives for the kurtosis: 4.55 and for the skewness -0.95. Also σ = 0.0136. Estimated densities obtained for this model with β = 37 are given in Figure 11. The dotted line is the result of the parametric estimation and the plain line of the nonparametric estimation, plotted on the histogram of the A j (T )'s. The nonparametric estimation detects a left tail that is not detected by the parametric one. Otherwise both estimators are very close.

The OU model with random= 1 is then validated with valid function. Figure 12 illustrates the result for a random trajectory (number 232): 100 simulated trajectories (black) and true trajectory (X 232 , red) (left plot) and quantiles of the true measurement among the 100 simulated points at each time points versus uniform quantiles. The qq-plot is satisfactory (compared to the graph obtained on simulated data Figure 5). Finally some prediction plots are performed (not shown) with the pred method and they confirm that model OU with random = c(1,2) with the parameters obtain from the parametric estimation, and the OU model with random = 1 and β = 37 produce very close trajectories and could be both validated.

We then apply the Bayesian procedure. As already mentioned, for the Bayesian procedure, large data sets are a problem because of the very long running time. Therefore, we thin the data

Discussion

In this paper we illustrate the functionality of the package mixedsde for inference of stochastic differential equations with random and/or fixed effects. This package, and mainly the function misedsde.fit, can be used to choose the best model to fit some data. It allows to compare two models: OU or CIR with one or two random effects. The three estimation methods can be used to help the decision maker. Nevertheless each method can be more appropriate to a specific situation, as explained before: the Bayesian method is recommended for a small number of observations, the frequentist nonparametric is a good tool with two random effects and no prior available. In particular the frequentist parametric proposes for a large sample, an estimation of the fixed effect and of the parameters of the Gaussian distribution for the fixed effect when there is only one. A neuronal dataset is studied with the three methods. To enrich this package and this study, one could implement the parameter estimation method developed in [START_REF] Delattre | Mixture of stochastic differential equations with random effects: application to data clustering[END_REF] for random effects distributed according to a Gaussian mixture distribution. Furthermore, other real data should be investigated with the present package.

  ; Picchini et al. (2010) (Hermite expansion of the likelihood); Delattre et al. (2013) (explicit integration of the Girsanov likelihood) or Delattre et al. (2016) (mixture of Gaussian distributions for the random effects); for parametric Bayesian estimation, we can cite Oravecz et al. (2009) (restricted to Ornstein-Uhlenbeck) and Hermann et al. (2016) (general methodology); for non-parametric estimation, we can cite Comte et al. (2013); Dion (2014); Dion and Genon-Catalot (2015) (kernel estimator and deconvolution estimators).

  Prediction (see Section 2.4) is implemented in function pred. Main inputs of the function are • Freq.fit or Bayes.fit object • invariant TRUE if the new trajectories are simulated according to the invariant distribution • level level of the empiric prediction intervals (default 0.05) • plot.pred TRUE to generate a figure (default TRUE)

Figure 1 -Figure 2 -

 12 Figure 1 -Simulated example 1 (CIR with one Gamma random effect), nonparametric estimation.Left: histogram of estimated random effects (A j ) and nonparametric estimation of f . Right: qqplot of (A j ) versus a Normal sample (true distribution is Gamma)

Figure 1

 1

Figure 5

 5 are obtained by choice1 <-floor(runif(1,1,M)) validationCIR <-valid(estim_nonparam, numj=choice1) plot(estim_param) valid(estim_param)

Figure 3 -Figure 4 -Figure 5 -Figure 6 -

 3456 Figure 3 -Simulated example 2 (OU with one Gaussian random effect) frequentist parametric estimation. Left: histogram of the (A j ) and Gaussian parametric estimation of f . Right parametric qqplot of (A j ) versus a Normal sample

Figure 7 -

 7 Figure 7 -Simulated example 2 (OU with one Gaussian random effect) Bayesian estimation. Markov chains of µ 1 , β, ω 2 1 and σ 2

Figure 7

 7 Figure 7 is produced by plot(estim_bayes) Traces of the Markov chains of µ 1 , β, ω 2 1 and σ are plotted, showing that all chains converge and have the correct location. Command print() yields acceptance rates of the MH algorithm:

Figure 8 -

 8 Figure 8 -Simulated example 2 (OU with one Gaussian random effect) Bayesian estimation. Left: predicted trajectories in black and 95% prediction interval in grey (green). Right: coverage rates: amount of observed values covered by the prediction intervals

Figure 11 -

 11 Figure 11 -Neuronal data, OU model, α random, β fixed, frequentist approaches. Estimation of the density f : parametric estimation in blue (black) plain line, non-parametric estimation blue (black) dotted line and histogram of the A j 's

Table 1 -

 1 Summary of the different methods for the two S4-classes Freq.fit and Bayes.fit resulting of the function mixedsde

	class	Freq.fit	Bayes.fit
	method	out	out
	method	plot	plot
	method	-	plot2compare
	method	print	print
	method summary	summary
	method	pred	pred
	method	valid	valid

  1 numj (t k ), . . . , X Mrep numj (t k )). If plot.valid=1, two plots are produced. Left: plot of the Mrep new trajectories (black) and the true trajectory number numj (in grey/red). Right: quantile-quantile plot of the quantiles of a uniform distribution and the N quantiles obtained comparing X numj (t k ) with the Mrep simulated values (X 1 numj (t k ), . . . , X

	Mrep
	numj

Acknowledgements

The authors would like to thank Vincent Brault and Laurent Bergé for technical help on the package. This work has been partially supported by the LabExPERSYVAL-Lab(ANR-11-LABX-0025-01).

The second author, Simone Hermann, was financially supported by Project B5 "Statistical methods for damage processes under cyclic load" of the Collaborative Research Center "Statistical modeling of nonlinear dynamic processes" (SFB 823) of the German Research Foundation (DFG).

Based on the best model selected by the frequentist approach, the OU model with one random effect φ j = α j is fitted. No prior knowledge is available, we therefore leave this information out and let the algorithm take the first 10%, i.e. 24, series for the calculation of the prior parameter, as described in Section 3.2. Figure 13

Appendix

When there is one random effect, what is the likelihood function and the MLE of the fixed effect?

Assume that we are in the case of random = 1, thus the process is

Let us compute the log-likelihood function when φ j = ϕ fixed. We omit the subscript j in the following. We use the same notation as for random=c (1,2), where V is a symmetric matrix size 2 × 2. We have :

We assume that the random effect is Gaussian with density f ξ , and denote ξ = (µ, ω) and θ := (µ, ω, α). Thus,

We find:

Finally after simplification we get:

Thus for random=1 we get

(9) Then, when random = 2 the roles of α and ϕ are exchanged. To implement a general formula, we note: r for random: 1 or 2, and c for the number of the common effect. We denote ψ the fixed effect and we the get the general formula:

[µ -V -1 r,r (U r -ψV c,r )] 2 + (U r -ψV r,c ) 2 2V r,r .

(10)


Benoit Arbelot

Romain Vergne

Thomas Hurtut

Joëlle Thollot

Pierre Bénard

Holger Winnemoeller

Automatic Texture Guided Color Transfer and Colorization

Keywords: Descriptors (according to ACM CCS): I, 4, 3 [Computer Graphics]: Enhancement-Filtering I, 4, 7 [Computer Graphics]: Feature Measurement-Texture

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In this paper, we propose a method to automatically apply local color transfer and colorization between images. Manually colorizing a grayscale image, or tuning colors to obtain a desired ambiance is challenging, tedious and requires advanced skills. Exemplar-based methods offer an intuitive alternative by automatically changing colors of an input image according to a reference image (the exemplar) containing the desired colors. The main challenge of these methods is to accurately match content between the input and reference image.

The first color transfer algorithms were based on global approaches reshaping the input image color histogram to match the histogram of the reference image. While these approaches can be simple and successful with carefully chosen image pairs, they often mismatch regions in the input and reference images, and are not suited for the colorization problem when the input image does not have a color histogram to begin with.

Alternatively, local approaches (soft-)segment an image into several subregions that can be processed independently. Colors are then added or transferred between similar regions. Those regions can be either manually provided, or automatically computed based on image descriptors.

Our approach is automatic and relies on regions defined as areas of similar textural content. This choice was driven by the fact that textures can be found everywhere in nature, and thus in a lot of photographs. Moreover, perceptual studies showed that the early stages of human vision are composed of several filters to analyze textures and color variations in our visual field [YJ * 93,Bal06]. This suggests that textures are important when observing images and should be a pertinent basis for local color transformations. Furthermore, textures can be efficiently described by a summary of first and second order statistics, and present an attractive middle ground between low-level descriptors (luminance, chromaticity) that cannot efficiently describe textured regions, and high-level descriptors (object and region semantic) that are complex, error-prone and slow to compute.

To apply color transfer between textured regions, our descriptors are computed on a large scale to be able to characterize large textures, but they must also preserve image structures. Existing methods for texture and structure decomposition are not well suited for our application: edge-aware image descriptors (such as bilateral filtering) have trouble analyzing highly contrasted textures and may introduce discontinuities in the color transfer. The alternative consists in detecting variations of the descriptors themselves (such as region covariance), but in that case, image edges are smoothed, leading to halos in the transfer.

Our solution to estimate texture properties is based on a texture analysis, followed by an edge-aware processing to compute edgeaware texture based descriptors. Our main contribution is to compute accurate textural information while preserving image structure. We use it in a generic framework for local color transfer and colorization between images based on textural properties.

Related Work

In this section, we review previous work on color transfer and colorization, before discussing several approaches to extract and analyze textures for image manipulation.

Color Transfer. An extensive review of color transfer methods can be found in [FPC * 15]. Color transfer consists in changing the colors of an input image to match those of a reference image. It was first introduced in [RAGS01] as a simple histogram reshaping, where the mean and variance of each channel are transferred separately, using the decorrelated Lαβ color space. This rather straightforward method can be surprisingly effective with well chosen input images. A rotation component was added in the matching process by Xiao and Ma [XM06], allowing the transfer to be done in a correlated color space (such as RGB). Instead of processing each channel independently, Pitié et al. [START_REF] Pitié | Automated colour grading using colour distribution transfer[END_REF] proposed to tightly match the 3-dimensional histograms using iterative 1-dimensional matchings. While the matching offered by this approach is very good, it is almost "too good" for the color transfer application as it tends to produce artifacts by forcing the input to have exactly the same number of pixels of each color as the reference. Finally, a more recent approach based on multiscale histogram reshaping was proposed in [START_REF] Pouli | Progressive color transfer for images of arbitrary dynamic range[END_REF] where the user can control how tightly the histograms should be matched. Overall, these global methods are simple, but histogram matchings do not ensure colors to be transferred between similar regions. When such automatic methods fail, manual segmentations can be provided to locally transfer between selected regions [DX09, AP10, LSZ12].

In order to automatically apply a local color transfer, Tai et al.

[TJT05] used mixtures of Gaussians to segment the input images and transfer colors between regions of similar luminance. A method to color grade videos based on color transfer between sequences was proposed in [BSPP13]. Their color transformation segments the images using the luminance and transfer chrominance between shadows, mid-tones and highlight regions. In a similar vein, Hristova et al. [HLMCB15] partition the images into Gaussian distributed clusters considering their main features between light and colors. Color-based segmentation was also used in [FSDH15] to extract color palettes and transfer between them using optimal transportation. While more accurate than global transfers, these approaches are still only based on first order information to segment the image and do not take higher order information to match regions between images. Consequently, regions with different textural properties but similar luminance cannot be distinguished.

Other approaches similar to Image Analogies [HJO * 01] have been applied to color transfer [SPDF13, OVB * 15]. However they differ from our approach as they use an additional input to compute the transformation.

Colorization. Colorization deals with the problem of adding colors to a grayscale image. One of the first approaches to tackle this issue relies on user input scribbles being extended via optimization across regions of similar luminance [START_REF] Lischinski | Colorization using optimization[END_REF]. This optimization is used with automatically generated scribbles in a lot of examplebased colorization methods [ICOL05, GCR * 12, KCP15]. Because they rely on a luminance-based optimization in their final step, these methods tend to have trouble with highly contrasted textures where the optimization does not propagate colors properly. More recently, Jin et al.

[JCT14] proposed a randomized algorithm to better match color distributions between user segmented regions.

Closer to our approach, some other methods rely on higher-order information to transfer the chrominance between pixels containing similar statistics [WAM02, CHS08, BT12, BTP14]. However, they often produce halos due to the window used in the statistics computation. These methods also rely on an energy minimization which typically makes them slow and hard to use on large images.

Texture Analysis. Many different descriptors have been used to manipulate images according to their textural content. Previous automatic colorization methods used SURF, Gabor features, or the histogram of oriented gradients as base tools for texture analysis [CHS08, GCR * 12, KCP15]. These descriptors are known to be discriminative, but also computationally and memory intensive due to their high number of features. Similarly, the shape-based texture descriptors introduced in [XDG10, XHJF12], although offering multiple invariants, are too complex for an image manipulation application where we expect to compute results in a reasonable time for relatively large images. The recent approaches proposed in [XYXJ12,CLKL14] precisely separate texture from structure using a relative total variation, but their descriptors are not accurate enough to discriminate textures among themselves. Finally, Karacan et al. [KEE13] proposed to use region covariance as a texture descriptor for image smoothing. Our method also relies on a variant We finally use these distance maps for both color transfer (C1) and colorization (C2), where attributed colors depends on pixel similarities.

of this descriptor, as it is compact and efficient in describing textural properties. One main drawback is that most of these descriptors tend to be unreliable around image edges and texture transitions, especially when estimated on large neighborhoods. For that reason, we also briefly describe edge-aware filtering methods that could be used to solve this issue.

Edge-aware filters are crucial to preserve image structures when smoothing, denoising, enhancing details, or extracting textural information from images. A well known approach regarding that goal is the bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF], which efficiently smoothes images while mostly preserving luminance edges. However, it tends to locally introduce halos and gradient reversal artifacts which can modify textural properties. The guided filter [START_REF] He | Guided image filtering[END_REF] offers a different approach by using a linear transform of a guidance image to filter an image but may also produce halos around edges. The anisotropic diffusion [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] or the unnormalized bilateral filter [APH * 14] are more appropriate for our descriptors, since they avoid both halos and gradient reversal when large scale diffusions are needed.

Overview

Our approach for automatically editing image colors based on textural content is summarized in Figure 2. First, descriptors are computed for the input and reference images in three steps (A): covariance matrices of several local image features are computed over a coarse scale to roughly characterize the textural content of each region (A.1). A multi-scale gradient descent then locally displaces descriptors in order to recover texture edges lost during the coarse scale analysis (A.2). Finally, an edge-aware filter is applied to obtain descriptors that accurately discriminate homogeneous textural regions while preserving detailed texture transitions (A.3).

Our descriptors allow the computation of similarities between pixels. As such, they also enable soft segmentations of the input and reference images, where smooth and sharp structures are preserved. This is illustrated in Figure 2 (B), where the vegetation is automatically isolated in both the input and reference images. Finally, similarity maps locally control the transfer of colors between images (C1) or colorize regions according to similar textural content (C2). The remainder of the paper is organized as follows: Descriptors are described in Section 4 and local color manipulation algorithms are detailed in Section 5. Results and comparisons are then presented in Section 6 before concluding in Section 7.

Edge-aware Texture Descriptors

Local Texture Descriptors

We want to analyze the textural information surrounding each pixel in both the input and the reference images. To that end, we chose to use region covariance [TPM06, KEE13] as it is an efficient and compact way of describing image regions. Region covariance captures the underlying texture by computing a small set of second order statistics on specific image features such as the luminance or the gradient. Let us consider a pixel p, described by a d-dimensional feature vector z(p). The region covariance is defined as the following d × d covariance matrix:

Cr(p) = 1 W ∑ q∈N p r (z(q) -µ µ µr)(z(q) -µ µ µr) T wr(p, q),
where N p r is a square neighborhood centered on p of size (2r + 1) × (2r + 1) and µ µ µr is a vector containing the mean of each feature inside this region. Unlike [START_REF] Tuzel | Region covariance: A fast descriptor for detection and classification[END_REF], we add a Gaussian weighting function with standard deviation r/3 that ensures de- scriptors to be smoothly defined from pixel to pixel: wr(p, q) = exp( -9 q-p 2 2r 2

) . Note that this weight function should also be used to compute the mean features µ µ µr. W is the normalization factor: W = ∑ q∈N P r wr(p, q). We typically use r ∈ [20, 30] and rely on a 6dimensional feature vector based on luminance derivatives to capture coarse scale textural content on natural images:

z(p) = L(p) ∂L(p) ∂x ∂L(p) ∂y ∂ 2 L(p) ∂x 2 ∂ 2 L(p) ∂y 2 ∂ 2 L(p) ∂x∂y T ,
where L(p) denotes the luminance of pixel p. In practice, each feature is first centered and normalized (i.e. we substract its mean and divide by its standard deviation) to equally contribute to the analysis. Note that other features, such as color derivatives, could also be used for different applications. For our color manipulations, we found that luminance carried most of the relevant texture information, especially in natural images.

As explained in [HCS * 09, KEE13], region covariances only describe second-order statistics, which can be a limitation when describing textural content as it cannot separate two distributions which only vary with their mean. Moreover, computing distances between covariance matrices is expensive because they do not lie in a Euclidean space. We thus follow the solution proposed by Karacan et al. [KEE13] who use the Cholesky decomposition to transform covariance matrices into vectors that can be easily compared and enriched with first-order statistics. Our descriptor is then rep- resented by:

S S Sr = L 1 r • • • L d r µ µ µr , (1)
where L i r is the i th column of the lower triangular matrix Lr obtained with the Cholesky decomposition Cr = LrL T r at scale r and µ µ µr are the first-order mean features in the corresponding region.

Visualizations of our descriptors are shown in Figure 3 where we can see that their values are similar when computed on the same types of regions. On the other hand, these values are dissimilar between different regions, making our descriptor able to discriminate different textural regions. Figure 4 shows how descriptors are affected by the scale r. Small scales (b) preserve edges but tend to produce noisy descriptors. Conversely, larger scales successfully describe uniform regions but fail to accurately preserve sharp texture transitions that often occur inside images. This is shown in (c), where the sharp transition between trees and sky is blurred when computing the descriptor with a large neighborhood. This phenomenon is due to the fact that on these particular pixels both tree and sky features are mixed to compute the descriptor, which then tend to represent this transition as a third texture. However, this is problematic for our color manipulation applications, where such descriptors will produce halos around edges. Note that we cannot integrate luminance edges in the weight function wr (as in the bilateral filter for instance). Indeed, this would prevent highly contrasted textures to be accurately captured since such textures would be fragmented into multiple pieces. For our purpose, we need both constraints to be satisfied: homogeneous descriptors inside regions and sharp texture edges preserved.

Multiscale Gradient Descent

To prevent texture transitions from being blurred, we propose to use a multiscale gradient descent algorithm to give these regions c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association. valid descriptors. Intuitively this multiscale gradient descent locally propagates relevant descriptor values (occurring inside homogeneous textural regions) to replace irrelevant ones (occurring around region borders). In order to do so, we use the variance of the descriptors to guide a gradient descent as this variance is low on homogeneous regions and high around texture edges. This gradient descent will then replace descriptors with high variance by those contained in uniform regions. Formally, the variance of a pixel p is computed as follows:

Vr(p) = 1 W ∑ q∈N p r
(S S Sr(q) -ν ν νr)(S S Sr(q) -ν ν νr) T wr(p, q) ,

where S S Sr(p) is the descriptor at pixel p and ν ν νr is the weighted average of the descriptors over the neighborhood N p r . The gradient descent replaces the descriptors on either side of the variance (e.g. texture edges) by descriptors with lower variance, consequently sharpening descriptor edges. Figure 6 (top) shows the pseudo-code of the gradient descent, where the returned map contains the coordinates of the descriptor that should be used for each pixel. The result is shown in the top row of Figure 5, where initial descriptors (a) are replaced by descriptors from homogeneous regions by following the gradient of the variance (b). The result obviously depends on the scale at which descriptors are computed. On large scales, complex texture transitions are smoothed out and consequently, some descriptors might be incorrectly attributed to different regions. This is illustrated in the bottom row of Figure 5, where the red pixel located in the sky (a) is mistakenly associated with the descriptor of a tree (b) after the gradient descent pass. Our solution to preserve complex texture changes with large scale descriptors is to use a multiscale gradient descent, where the scale of both descriptor and variance are gradually increased to guide the gradient descent of the initial (large scale) descriptor.

Figure 6 (bottom) shows the pseudo-code of the proposed multiscale gradient descent process. The idea is to iteratively apply gradient descents, from fine to coarse scales, in order to propagate desciptors from homogeneous regions while preserving complex texture edges. At small scales, the descent accurately preserves edges, but quickly falls into local minima. Increasing scales slowly select pixels further and further away from the detailed edges, ensuring that the descriptors are consistent. In practice, the number of iterations used for a given scale is set to the size of the neighborhood (small and large scales may respectively lead to small and large propagations). Note that, even if small scale descriptors are needed to compute the variance, the resulting new coordinates only modify the coarse scale descriptor. The result is shown in Figure 5 (c). The obtained descriptor (top) better preserves complex texture transitions. The red pixel (bottom) now successfully takes descriptor values of a homogeneous region inside the sky.

Unnormalized Bilateral Filtering

Gradient descent ensures the precise capture of textural properties around each pixel, even near texture edges. Yet, descriptors might still contain some variations that do not appear in the original image. These might happen around U-shaped texture transitions (as in the left part of Figure 5 (c)) or when a region cannot be properly defined by its textural content (such as a fine edge on a uniform background). This has to be prevented since any variations in the descriptors might lead to color changes during transfer or colorization. In a last step, we thus smooth the descriptor using an edge-aware filter to perfectly fit to the image structure. To that end, we adapt the unnormalized bilateral filter [APH * 14], such that it iteratively smoothes the descriptor according to luminance variations. This filter is simple, efficient, and introduces very little halos if any. However, any other edge-aware filter could have been used [TM98, HST13, PM90]. Formally, we use the unnormalized bilateral filter as follows:

S S S ub f (p) = S S S(p)+ ∑ q∈N p Gσ s (q -p)Gσ l (L(q) -L(p))(S S S(q) -S S S(p)) 2πσ 2 s , (3)
where Gσ(x) = exp(

-x 2 2σ
2 ) is a standard Gaussian kernel. σs and σ l respectively control the influence of spatial distances and lumi- nance variations. In practice, we iteratively apply Equation (3) with rather small values of σs and σ l (typically 2 and 0.05) in order to accurately diffuse descriptors on large neighborhoods. Figure 7 shows the effect of the filter on a problematic region, where the descriptors do not precisely follow edges around the palm tree (a). The unnormalized bilateral filter accurately brings back the leaf edges, as shown in (b). The last image (c) shows the effect of the filter when applied on the original descriptor (i.e. without gradient descent). In that case, halos are propagated inside regions and create unreliable descriptors.

Effect of Each Step on Color Transfer Result

An example showcasing the effect of each of the previous steps on a color transfer result is given in Figure 8. Initial descriptors (a) tend to produce strong color halos around texture transitions. The multiscale gradient descent (b) produces much sharper transitions, but may also introduce discontinuities around U-shaped regions. When applied alone, the unnormalized bilateral filtering (c) smudges halos instead of suppressing them. Combining multiscale gradient descent and unnormalized bilateral filtering (d) creates a clean result where even strong initial halos are efficiently removed.

Local Color Manipulation

Now that we have obtained reliable descriptors, we propose to use them for color manipulations by defining transfer functions that only rely on similar pixels between the input and reference images.

Pixel Similarity

We define a similarity measure based on the L2 Euclidean distance between two descriptors: where S S S(p) and S S S(q) are the descriptors at locations p and q and σ d is the standard deviation that controls how close descriptors should be to contribute to the similarity measure. Note that other metrics could have been used as detailed in [HCS * 09, KEE13], but we did not find any significant differences for our purpose. An example of similarity measure is shown in Figure 9, where pixels (b), (c) and (d) are compared with all the other pixels of the input image (a). We can observe that trees, sky and grass regions are accurately selected and distinguished in the results.

Dσ d (p, q) = exp -S S S(p) -S S S(q) 2 2σ 2 d , (4)

Color Transfer

The main idea for transferring colors between images is to rely on local histogram matchings between input and reference images, where both sets of color points are defined by their texture similarities. The matching process is based on a translation and scaling of the distribution in a decorrelated color space, as originally proposed by Reinhard et al. [START_REF] Ashikhmin | Color transfer between images[END_REF]. Input and reference images are therefore first transformed into the uncorrelated and perceptually uniform CIE-Lab color space before being processed. The following transfer function is then applied on each channel c ∈ {L, a, b} separately:

Tσ d (p) = std re f (p) std in (p) c in (p) -µ in (p) + µ re f (p), (5)
where superscripts "in" and "re f " denote the input and reference images. "µ" and "std" are the weighted mean and standard deviations respectively, computed as follows, according to the similarities of the pixel p of the input image:

µ img (p) = 1 W ∑ q c img (q)Dσ d (p in , q img ) std img (p) = 1 W ∑ q (c img (q) -µ img (p)) 2 Dσ d (p in , q img ),
where img ∈ {in, re f }, q iterates over img and W is the normalization factor: W = ∑ q Dσ d (p in , q img ). A color transfer example is shown in Figure 10 (top) where we can observe the effect of the σ d parameter. When σ d is small, colors are transferred only between highly similar regions, such as the sea or the clouds of the input and reference images here. Wider and wider regions are considered when increasing σ d , leading to results closer to the global matching of [START_REF] Ashikhmin | Color transfer between images[END_REF].

Colorization

Histogram matching techniques cannot be used directly for colorizing images that do not contain chrominance channels. In this case, we assign the mean chrominance of the reference image to each input pixel, weighted by our similarity measure:

Cσ d (p) = ∑ q c re f (q)Dσ d (p in , q re f ) ∑ q Dσ d (p in , q re f ) . ( 6
)
Note that this transfer function is applied on chrominance channels only, although the luminance could also be modified depending on the purpose. A colorization example is shown in Figure 10 (bottom). Large values of σ d tend to average colors on large regions and consequently create pale and monochrome results. Therefore σ d should be kept small enough for colorization purpose, in order to only average colors over regions of highly similar descriptors.

Implementation & Performances

We fully implemented our color manipulation functions on the GPU using Cuda. All the results presented in this paper were obtained with a NVIDIA Quadro 6000 graphics card. In practice, we first precompute the descriptors S S S for both the input and reference images before applying a transfer or a colorization. However, Equations (5) and 6 require to iterate over all the pixels of the input image, and compute the similarities with the whole reference for each of them in order to obtain the weighted mean and standard deviations. A naïve implementation of these equation leads to extensive computation times.

To achieve reasonable speed, we propose to quantify similarities using a user-defined distance τ that controls how close two descriptors should be to be considered as equal. Considering a particular τ = 0, 20min τ = 0.01, 21s τ = 0.1, 6s τ = 0.5, 1s

Figure 11: Optimization impact. Color transfer results for increasing τ values for 512 × 512 images. The lower τ, the higher the speed-up and the probability of quantization artifacts. In this example, τ = 0.1 allows for a fast transfer which can be used for efficient results exploration with minimal visual artifacts.

input pixel p, all the other pixels p i such as Dσ d (p, p i ) < τ are processed using the same similarity function. That way, increasing τ decreases the total number of iterations needed to obtain the result.

The effect of this optimization can be seen in Figure 11, where important speed-up is achieved without visual impacts. High values of τ tend to produce quantization artifacts, but may be used to interactively explore the result space.

To summarize, the user can tune the following parameters to achieve the desired results:

• rmax controls the size of the window on which descriptors are computed and thus defines the scale at which textures are estimated. Typically, we found that rmax = 21 works well for natural images of resolution 512 × 512. • σs and σ l respectively control the influence of spatial distances and luminance variations when smoothing the descriptor with the unnormalized bilateral filter. All the results in the paper were done with σs = 2 and σ l = 0.05. The number of iterations used for this filter depends on the complexity of texture edges. We typically used 500 iterations for our results. • σ d controls how strongly the weight between two pixels is influenced by their distances in the descriptors space. In practice, we respectively used σ d = 1 and σ d = 0.2 for most color transfer and colorization results. • τ controls the quantization step. In our results, we used τ = 0.01 as it provides a good speed-up while keeping a good visual quality in almost every case.

The timings of our algorithm for the image in Figure 11 where "Desc." stands for the descriptors computation and "Transfer" stands for the color transfer or colorization step. Typically the colorization takes longer because of the lower σ d used which create less similarity between pixels (see Equation (4)), leading to more computation during the transfer step. Note that our code was designed to be strongly flexible, while maintaining decent timings as much as possible. We believe further Input Reference Our Result optimization could still provide significant speed-up. For example a bilateral grid could be used to compute the filtering step much faster during the descriptors computation. The transfer step could also be further optimized by considering subsampled versions of the reference and input images which would greatly reduce the transfer time for large images. However this unoptimized code still allows for computation times similar to other color transfer or colorization methods relying on image descriptors.

Results

Results and comparisons presented in the paper and in the supplemental materials were all made with the default parameters given in the previous section.

Color Transfer Results

Figure 13 (top) shows the results of our color transfer against other state-of-the-art methods. The results of [START_REF] Ashikhmin | Color transfer between images[END_REF] were computed with our own implementation of their method. The results of [START_REF] Pitié | Automated colour grading using colour distribution transfer[END_REF][START_REF] Pouli | Progressive color transfer for images of arbitrary dynamic range[END_REF] were computed using the available code on the authors webpage, we used a full match (100%) for [START_REF] Pouli | Progressive color transfer for images of arbitrary dynamic range[END_REF]. The results of [FSDH15] were provided by the authors. The results of [HLMCB15] were taken from the authors webpage and drove our choice of images.

These results show that global approaches [RAGS01, PKD07] tend to produce saturated colors due to the stretching of the input color histogram. Furthermore, global histogram matchings match regions of similar colors and luminance, failing in transferring colors between similar textured regions if they have highly different luminance or colors. This is showcased in the bottom row where the orange color of the reference buildings is transferred to the input sky. The progressive approach of [PR11] also fails to accurately preserve the colors of the reference in their results. Local approaches based on color information [FSDH15,HLMCB15] lead to better results, but also fail in matching regions of similar textural content because they define similar regions by their luminance and color distributions.

Our approach successfully matches those regions, as shown in the third row, where the flower field of the reference is matched to the grass of the input (making it yellow); or in the fourth row where the buildings of the reference are matched to those of the input (making them orange). Figure 12 shows two more examples where the matching between different regions is clearly effective thanks to our descriptors. Those results show that the method of [WAM02] based on luminance matching fails when the input images are too complex: different regions with similar luminance get the same colors, such as the building and clouds in the first example. The method of [CHS08] uses SURF descriptors and Gabor filters which are strongly discriminative, leading to efficient colorization when the input and reference images have identical or very similar content. However, they have to crop image borders and colors often smudge in their results. The method from [GCR * 12] produces better results by combining superpixel segmentation and similarly robust descriptors (i.e. image intensity, standard deviation features, Gabor filters and SURF features). While achieving better results than previous methods, they still fail to distinguish between intricate regions such as the clouds and the sky in the first and fourth row, or the river and the land in the fourth row. Finally, the method from [BT12] uses descriptors based on standard-deviation, discrete Fourier transform and cumulative histograms of image patches. It is very prone to halos due to the window used in the descriptors computation and was improved in [PAB * 15] where a new luminance-chrominance model was used to better propagate colors. While this model is very good to avoid artifacts, the final colors are not always faithful to the reference image colors, as seen in the first and third rows.

Colorization Results

As seen in the last column, our approach accurately matches cor- responding textures and produces colorful results: sky, cloud, vegetation, mountain and building colors of the references are successfully transferred into the input images. Figure 12 shows two more results demonstrating a clear separation between regions of the input image and correct color associations from the reference image.

To evaluate the coherency of our descriptors, we also tried to colorize a desaturated image using the original color image as reference. These results are shown in Figure 14, where we can observe that color differences between the reference and the output images depends on σ d : the lower σ d , the higher the fidelity. This is due to the fact that input and reference images have exactly the same descriptors in that case. In the limit case, when σ d → 0, only one pixel will be taken into account when comparing descriptors (cf. Equation ( 6)) and the result will be equal to the reference. The pixel-wise difference between the result and the reference image was computed as the sum of the absolute RGB differences. Note that, when input and reference images differ, σ d should be given a higher value to avoid color artifacts.

Combining Colorization and Transfer

Since our framework is the same for colorization and color transfer, we can easily apply a combination of both to a grayscale input by adding chrominance via colorization, while modifying the luminance by transferring only the luminance from the reference image. The results of this approach can be seen in Figure 15. They show that this combination can produce a result closer to the style of the reference image, while still using only the input luminance. Comparing this to the result of the color transfer (which also transfers luminance), we see that color transfer remains more colorful because the chrominance information of the input image is also used, however it requires a color image as input which is more restrictive.

Discussion and Future Works

In this paper, we presented a generic framework for both color transfer and colorization. Our edge-aware descriptor accurately captures similar textural content in images while being robust to texture transitions. It allows local color transfer and colorization between similar regions of an input and reference images. Our method suffers from two main limitations, as described below.

(1) Considering colorization, the input and reference images should be similar enough to produce coherent results. If a particular region in the input image does not have any correspondence in the reference one, the similarity function (based on a Gaussian distance) tends to give the same weight to all pixels, resulting in a monochrome colorization. Note that this is equivalent to increasing σ d for this particular region, as seen in Figure 10 (bottom-right). This problem also occurs for color transfer but is much less visible since the mean and variance are only used to modify the histogram. To prevent this, one possibility would be to automatically detect mismatched regions and ask the user to disambiguate the transfer by providing more specific reference images.

(2) The proposed descriptors efficiently capture texture regions and their transitions, but are not able to detect higher-level semantic information such as faces, man made objects or background and foreground. Our descriptors might be altered by such objects, thus affecting the quality of the transfers. Again, this is most visible in colorization results, as shown in Figure 16. The yellow color obtained in the top left part of the image is due to the electric wires that are associated to the warning sign of the reference. The wheels of the motorbike contain fine structures associated to the girl's hat, resulting in a bluish color. One way to mitigate these issues would be to rely on more complex, but slower, descriptors combining both semantic and texture information.

Input

Reference Our Result

Figure 16: Failure case. Semantic information such as man-made objects or faces may locally modify the descriptors and produce incoherent colorizations. For example here, the motorbike wheels are colored in blue because they are detected as a texture resembling the one of the input girl hat.

Despite these limitations, we believe that our descriptor constitutes a good basis that could contribute to other applications such as tone mapping, edge-aware image decomposition, and color content modification of videos. Our method could also easily be extended to handle scribbles in the descriptors processing step. This would allow the user to fine tune the descriptors and distance maps and customize them as he pleases.

Figure 2 :

 2 Figure 2: Pipeline overview. Edge-aware descriptors are first computed to accurately describe the textural content of the input and reference images (A). They are then used to compute per-pixel distances and allows similar regions to be associated, as shown for the vegetation in (B).We finally use these distance maps for both color transfer (C1) and colorization (C2), where attributed colors depends on pixel similarities.

Figure 3 :

 3 Figure 3: Texture descriptors. Patches taken from several regions of the image in Figure 2 (top) and their respective descriptors computed for the central pixel of the window (bottom). Yellow and blue values correspond to positive and negative values respectively. Patches from similar regions have similar descriptors.

Figure 4 :

 4 Figure 4: Descriptors scales. Small scales lead to noisy descriptors (b). Large scales lead to more homogeneous descriptors and smooth sharp texture transitions. For visualization clarity, only the first element of S S Sr is shown (i.e. the first value of L 1 r ) but the rest of the set presents the same behavior.

Figure 5 :

 5 Figure 5: Gradient descent illustration. (a) A zoom in the sky/trees transition of the image shown in Figure 4. (b) A gradient descent guided by the variance of the coarse scale descriptor tends to sharpen edges (top), but may mistakenly assign descriptors to the wrong side of the edges: The red sky pixel (bottom) is considered as part of the trees here. (c) A gradient descent gradually performed at multiple scales (from fine to coarse) better preserves complex texture transitions. The red pixel is now successfully assigned to the sky.

Figure 6 :Figure 7 :

 67 Figure 6: Multiscale gradient descent algorithm.

  c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association.

Figure 8 :

 8 Figure 8: Processing effect on the transfer result. In this example, initial descriptors are blurry and create strong color halos above the trees in the transfer result (a). The Multiscale Gradient Descent (MGD) prevents the apparition of halos but some incorrect edges remain in U-shaped transitions between the sky and the trees (b). The Unnormalized Bilateral Filtering (UBF) accurately preserves the structure but smudges halos instead of suppressing them when used alone (c). The combination of both MGD and UBF leads to a cleaner result as shown in (d).

Figure 9 :

 9 Figure 9: Similarity maps. (a) Input image luminance. The green, yellow and red pixels are compared with all pixels using Equation (4) to obtain the corresponding similarity maps (b), (c) and (d). The similarity measure allows the three regions to be accurately discriminated. Similarities were computed with σ d = 1 in these examples.

B= 2 Figure 10 :

 210 Figure 10: Impact of σ d on transfer functions. Top: color transfer example. When increasing σ d , more and more pixels are considered as similar, resulting in a transfer close to a basic global histogram matching. Bottom: colorization example. As colors are obtained from the weighted average of similar pixels in the reference image, increasing σ d tends to produce a monochrome result.

Figure 12 :

 12 Figure 12: Transfer and colorization results. Different colors are clearly associated with different regions based on their textural content.

Figure 13 (

 13 Figure 13 (bottom) compares the results of our colorization against other state-of-the-art methods. The results of [WAM02, CHS08, GCR * 12] were taken from [GCR * 12]. The results of [BT12, PAB * 15] were computed using the code provided by the authors, using the default parameters suggested in their code.

Figure 13 :Figure 14 :

 1314 Figure 13: Comparison with previous methods. Top and bottom respectively compare color transfer and colorization results with previous state-of-the-art methods. See the text for more details.

Figure 15 :

 15 Figure 15: Combining colorization and luminance transfer. Our framework allows for an easy combination of colorization and luminance transfer. This combination provides a good style transfer between the input and reference images. While less colorful than a color transfer result, this result only requires a grayscale input. In those results, σ d = 0.5.

  using those parameters are described in the following table:

	Algorithm	Image Size	Desc. Transfer Total
	Colorization	512 × 512	19s	47s	66s
	(σ d = 0.2)	1024 × 1024	78s	490s	568s
	Color Transfer	512 × 512	19s	21s	40s
	(σ d = 1)	1024 × 1024	78s	132s	210s

  c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association.

	Input	Reference	[RAGS01]	[PKD07]	[PR11]	[FSDH15]	[HLMCB15]	Our Result
	Input	Reference	[WAM02]	[CHS08]	[BT12]	[GCR

*

12] [PAB

* 15]

Our Result

c 2016 The Author(s) Eurographics Proceedings c 2016 The Eurographics Association.

Acknowledgements

We thank the anonymous reviewers for their relevant comments. We also thank Oriel Frigo for helping generating the comparison results in Figure 13, and Olivier Le Meur for his helpful comments. Finally, for providing access to their code, we thank Fabien Pierre, Raj Kumar Gupta and all authors whose code was accessible online. The input image used in Figures 2,3, 4, 5, 9, 11, 15 and the house reference image in Figure 10 are courtesy of freebigpictures.com. Color transfer comparisons were made using the images from [HLMCB15]. Colorization comparison images were taken from [GCR * 12]. This work is partially supported by the French National Research Agency, MapStyle project [ANR-12-CORD-0025].