Variational Exemplar-based Image Colorization

Aurélie Bugeau, Vinh-Thong Ta, and Nicolas Papadakis Abstract-In this paper, we address the problem of recovering a color image from a grayscale one. The input color data comes from a source image considered as a reference image. Reconstructing the missing color of a grayscale pixel is here viewed as the problem of automatically selecting the best color among a set of colors candidates while simultaneously ensuring the local spatial coherency of the reconstructed color information.

To solve this problem, we propose a variational approach where a specific energy is designed to model the color selection and the spatial constraint problems simultaneously. The contributions of this paper are twofold: first, we introduce a variational formulation modeling the color selection problem under spatial constraints and propose a minimization scheme which allows computing a local minima of the defined non-convex energy. Second, we combine different patch-based features and distances in order to construct a consistent set of possible color candidates. This set is used as input data and our energy minimization allows to automatically select the best color to transfer for each pixel of the grayscale image. Finally, experiments illustrate the potentiality of our simple methodology and show that our results are very competitive with respect to the state-of-the-art methods.

I. INTRODUCTION

I MAGE colorization consists in recovering a color image from a grayscale one with initial color input data. This application attracts a lot of attention in the image-editing community in order to restore or to colorize old grayscale movies or pictures. A first category of colorization algorithms requires a huge amount of user intervention to manually add initial colors to the grayscale image. The colorization process is then performed by propagating the input color data to the whole image. The color propagation is done by using methods such as quadratic optimization, diffusion or variational techniques (see for instance [START_REF] Levin | Colorization using optimization[END_REF]- [START_REF] Buades | Conditional image diffusion[END_REF] and references therein). The main drawback of these methods is that the initialization of the colorization process is done completely manually.

Another colorization approach consists in transferring color from one (or many) initial color image considered as example. The color prediction is performed from the luminance (grayscale) channel of both images and the final color is computed with methods such as image statistics or machine learning approaches. For instance, in [START_REF] Welsh | Transferring color to greyscale images[END_REF] (and later [START_REF] Morimoto | Automatic colorization of grayscale images using multiple images on the web[END_REF] that uses multiple images as references) the color transfer is performed by minimizing a distance based on simple statistics on luminance images. In [START_REF] Irony | Colorization by example[END_REF], [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] and [START_REF] Liu | Automatic grayscale image colorization using histogram regression[END_REF], more sophisticated methods A. Bugeau is with University of Bordeaux, LaBRI, UMR 5800, F-33400 Talence, France. email: aurelie.bugeau@labri.fr. phone: +33(0) 5 4000 6900 fax: +33(0) 5 4000 6669

V.-T. Ta is with IPB, LaBRI, UMR 5800, F-33400 Talence, France email: vinh-thong.ta@labri.fr. phone: +33(0) 5 4000 6900 fax: +33(0) 5 4000 6669 N. Papadakis is with CNRS, IMB, UMR 5251, F-33400 Talence, France email: nicolas.papadakis@math.u-bordeaux1.fr phone: +33(0) 5 4000 2116 fax: +33(0) 5 4000 2123 are proposed based on numerous and complex steps such as image segmentation, classification, locally weighted histogram regression, machine learning techniques, color refinement and image descriptors to capture textures or complex structures. One drawback of such approaches for image colorization is the spatial coherency during the color transfer which can lead to possible inconsistent colorization in the final result. To overcome this limitation, [START_REF] Welsh | Transferring color to greyscale images[END_REF] adds user interaction in the colorization process and [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] uses a parameter to enforce the local spatial relationship between pixels.

In this paper, we focus on this second category of methodology based on reference images. We propose a simple method with very few steps which simultaneously performs the color transfer while enforcing the spatial coherency of the reconstructed colors. The experiments presented in this paper (Section IV) show that our simple methodology can be competitive with the aforementioned more complex approaches.

Exemplar-based Methods: Most of exemplar-based colorization techniques make the assumption that pixels with similar intensities or similar neighborhood should have similar colors. Such approaches are also called exemplar-based (or patch-based) methods. They have recently received a lot of attention in the image processing community. Indeed, first introduced for texture synthesis [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF], this concept has been derived for many applications such as image denoising or image inpainting where the corrupted or the missing data are restored from the values of the other pixels contained in the image. In the case of image denoising, the so-called NLmeans filter [START_REF] Buades | A non-local algorithm for image denoising[END_REF], which corresponds to a weighted averaging filter [START_REF]Image denoising methods : A review[END_REF], uses weights proportional to the similarity measure between the different patches of the image to restore. For image inpainting, the same idea is used. For a given missing pixel, the best corresponding pixel (or region) to copy/paste during the inpainting process is found by using exemplarbased similarities (see for instance in [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF], [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] and references therein).

The main underlying idea behind all these exemplar-based approaches is that the best a priori for an image is the image itself.

For all the aforementioned applications, the common problem can be stated as follows. Given a point defined in a domain, how to select the best match(es) from a set of unordered examples (defined in a d-dimensional vector space) in order to obtain the best reconstruction of the missing data ? In the case of non-local image denoising, it consists in computing the weights between the different candidates in order to realize a weighted average. For image inpainting or image colorization, this problem corresponds to choosing a single candidate to realize the copy/paste operation.

One drawback of such approaches is that discontinuities or spatial incoherence in the reconstructed image may be visible since most of the algorithms are of greedy type and each pixel is processed independently. In the case of image inpainting, several methods have been proposed to overcome this issue by optimizing an energy defined over the whole missing area [START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF]- [START_REF] Wexler | Space-time video completion[END_REF]. To the best of our knowledge, such strategy has not yet been studied in the specific context of image colorization.

Another particularity of exemplar-based methods is that they almost always use the SSD measure (sum of square differences) in order to compute a pixel-wise distance between the patches to compare. Nevertheless, natural images are composed of different textures and image structures which are not well captured by the SSD distance. To take into account the different type of information contained in natural images, it seems relevant to use and combine different patch-based features and their associated metrics.

Previous Works and Main Contributions: In a previous work [START_REF] Bugeau | Patch-based image colorization[END_REF], we have proposed an image colorization framework which combines different features and distances in order to select a set of different color candidates for each pixel of the grayscale image. For each pixel of the grayscale target image T , a set of possible colors is found in the color source image S. This set is constructed by comparing the luminance patches of both images according to different features and distances. One of the candidates is then selected by choosing the median value on the principal axis given by a dimensional reduction (PCA) on the different color candidates A color regularization step is finally performed to smooth this arbitrary point-wise choice and ensure the spatial coherency of the reconstruction.

This paper is an extension of the method proposed in [START_REF] Bugeau | Patch-based image colorization[END_REF] but relies on a framework that is fundamentally different. We address the exemplar-based image colorization problem under a variational point of view.

In this paper, contrary to [START_REF] Bugeau | Patch-based image colorization[END_REF], the candidate selection and color regularization problems are simultaneously solved through a variational energy minimization. Our new approach then avoids the arbitrary pre-selection of a single candidate. To that end, we design an energy whose minimizers automatically select the best color for each pixel from a set of candidates, while ensuring the spatial coherency of the reconstruction, as it enforces neighbor pixels to have similar colors in the final result.

Paper Organization: Section II reviews the general image colorization framework introduced in [START_REF] Bugeau | Patch-based image colorization[END_REF]. Section III presents a variational model for the automatic candidate selection and the related optimization algorithm. Experimental results and energy analysis are shown in Section IV and, finally, conclusions are given in Section V.

II. GENERAL FRAMEWORK FOR EXEMPLAR-BASED IMAGE COLORIZATION

Most of exemplar-based methods for image colorization proposed in the literature can be summarized into a general pipeline. We now detail this pipeline in order to introduce the proposed variational exemplar-based image colorization method. In the following S and T will respectively denote the color source image and the grayscale target image to colorize.

A. General Image Colorization Pipeline

Exemplar-based image colorization methods can be decomposed into three main steps.

1)

Step-1, Pre-Process S and T : This step consists in converting both S and T to a luminance-chrominance color space such as lαβ [START_REF] Welsh | Transferring color to greyscale images[END_REF], Lab [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] or Y CbCr [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF].

As mentioned in the introduction, exemplar-based image colorization methods rely on a copy/paste of the colors from the source image to the target image. More precisely, the chrominance channels will be transferred from one to the other. In order to select the color to transfer, the luminance of S and T should be comparable. Ideally, the transfer process should take into account the global difference in luminance of the two images, i.e., these two images should have similar global statistics such as the mean and the variance. For this purpose, a luminance mapping can be performed [START_REF] Hertzmann | Image analogies[END_REF].

2)

Step-2, Predict Color from S and Transfer to T : In order to have acceptable computation time in this second step, only a reduced set S n = {q 1 , . . . , q n } of n possible candidate pixels are extracted from the source image S with n ≪ #S, where #S corresponds to the cardinality of the set S.

The second step aims at choosing for each pixel p ∈ T the best pixel q ∈ S n which minimizes a distance ρ on the luminance neighborhoods

q = arg min q∈Sn ρ(p, q) .
Finally, the chrominance of q is transferred to p, so that T becomes a final true luminance-chrominance image.

3) Step-3, Post-Process T : This last step generally consists in a simple color space re-conversion, i.e., from the luminancechrominance to the RGB color space. A color refinement step may also be realized in order to enhance the final visual quality of the result [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF].

B. Proposed Image Colorization Pipeline

We now describe how each step of the pipeline is achieved in our image colorization approach.

1)

Step-1: As done in [START_REF] Bugeau | Patch-based image colorization[END_REF], we propose to convert S and T to the convex Y U V color space to obtain the luminancechrominance decomposition S = (Y S , U S , V S ) and T = (Y T , U T , V T ). This empirical choice comes from our experiments. Indeed, by using non standard and non convex color spaces, such as the lαβ [START_REF] Welsh | Transferring color to greyscale images[END_REF], our experiments have shown that such spaces can lead to numerical issues. More precisely, ass the projection on non convex sets is not unique, variational methods can finally provide visually inconsistent colorization results.

Finally, similarly to [START_REF] Welsh | Transferring color to greyscale images[END_REF], a luminance mapping [START_REF] Hertzmann | Image analogies[END_REF] is performed. The luminance Y S of the source image is modified in order to map the one of the target image Y T so that the two luminance images become comparable.

2)

Step-2: As described previously, only n ≪ #S samples are randomly selected from S on a jittered grid leading to the reduced set of candidate pixels S n = {q 1 , . . . , q n }. From now, only these n pixels of S will be used to colorize the target image T .

Our color prediction model consists in selecting the best chrominances to transfer to each pixel p ∈ T . Contrary to previous approaches in image colorization, we do not only predict one unique chrominance per pixel. Instead, several possible chrominances are kept according to different image features and distances. Natural images contain different types of complex structures, redundancies and textures. In order to perform colorization on natural images, every algorithm should integrate such information. By using several different features and distances, our process will adapt itself to the different properties of the considered image.

While many metrics to compare patches exist, we here only concentrate on L = 3 features f l and their associated distances ρ l with l = {1, 2, 3}. Each of these three features is computed for different size σ of patches. For each feature among all the possible sizes, the set of possible patch sizes is denoted as Σ l (with σ ∈ Σ l ). The choice of the metrics and patch descriptors has been made experimentally, but any other feature could be integrated straightforward.

Each chosen feature captures a different characteristic of the grayscale image. Given a luminance image Y , Y (p) = Y (x, y) is the luminance value of a pixel p ∈ Y . We denote a patch of size σ 2 centered at pixel p = (x, y) as

P σ (p) = {Y (x + k, y + l)} ∀k, l = -σ/2, . . . , σ/2
For any pixel p ∈ Y and for different size σ of patches, we compute: i) the variance

f 1 (p, σ) = ν(P σ (p))
which differentiates constant and textured areas; ii) the amplitude spectrum of the 2D discrete Fourier transform

f 2 (p, σ, ξ) = DFT(P σ (p), ξ) 2
where ξ = (ξ 1 , ξ 2 ) states as the frequency. This feature gives information on the high frequencies of the most important structures; iii) the cumulative histograms containing H bins with {h 1 , . . . , h H } and ∀i = 1, . . . , H

f 3 (p, σ, h i ) = i j=1 q∈Pσ(p) δ(Y (q), h j ) #P σ (p)
where δ(a, h j ) = 1 if the color a belongs to the jth bin h j , i.e., a ∈ h j and δ(a, h j ) = 0 otherwise. This cumulative histogram allows to take into account the luminance distribution inside the patch. From the features space, we define the associated distances ρ l with l = {1, 2, 3} such that for two pixels p ∈ T and q ∈ S n , we have:

ρ 1 (p, q, σ) = |f 1 (p, σ) -f 1 (q, σ)| , ρ 2 (p, q, σ) = ξ |f 2 (p, σ, ξ) -f 2 (q, σ, ξ)| and ρ 3 (p, q, σ) = H h=1 |f 3 (p, σ, h) -f 3 (q, σ, h)|.
As we are comparing 1-dimensional grayscale images, we can notice that the ℓ 1 distance ρ 3 on the cumulative histograms is equivalent to the ℓ 1 Wasserstein distance between the distributions.

From the introduced features and metrics, we can define the set of chrominances for each pixel p ∈ T from which the final chrominance will be chosen. For each point p ∈ T , each feature l = 1, . . . , L and each size σ ∈ Σ l , we select the candidate q l,σ such that:

q l,σ = arg min q∈Sn ρ l (p, q, σ).
Therefore, by finding the best source point in S n for all the L features and distances, we obtain the chrominance set C(p) for each p ∈ T . As we only intend to colorize the chrominance channels U an V of the Y U V image, the chrominance set is composed of d = 2-dimensional elements:

C(p) = {(U (q l,σ ), V (q l,σ )) , ∀l = 1, . . . , L and ∀σ ∈ Σ l } .
For the sake of clarity, we rewrite this set as:

C(p) = {c i (p) = (U (q i ), V (q i )) , ∀i = 1, . . . , N }
where N = L l=1 #Σ l . We now have a reduced set of possible chrominance for each pixel p of the grayscale image.

The best candidate for each point could be taken as the one having the smallest distance ρ l . Unfortunately, as the distances are of different nature, they are not directly comparable. Each feature will allow one to measure different image characteristics and we have no prior on which feature should be used in a particular area. In [START_REF] Bugeau | Patch-based image colorization[END_REF], the problem of selecting the best color to transfer among the set of possible colors is solved, for each pixel independently, using a dimensional reduction (PCA) of the chrominance candidates C(p). For each pixel, the median value on the main axis of the PCA is considered as the single color candidate. In this paper, in order to avoid such an arbitrary selection, we propose to use the spatial coherency to automatically chose the best candidate within a variational framework that will be presented in Section III.

3) Step-3: The post-processing step mainly aims at reconverting the luminance-chrominance image to the RGB color space. The generated color images are generally drab. In order to improve the shininess of the colored images, we here propose to first apply a mapping [START_REF] Hertzmann | Image analogies[END_REF] of the chrominances from U T to U S and V T to V S before converting T to RGB space: the chrominance histograms are linearly shifted and scaled to fit the histograms of S. Obviously this solution is not ideal as there is no reasons for these histograms to be completely similar at the end of the process but it gave us satisfactory experimental results.

III. VARIATIONAL EXEMPLAR SELECTION FRAMEWORK

In this section, we describe a variational method for solving the general following problem. Given a point p ∈ T , the problem is to select the best match from a set of unordered examples C(p), defined in a d-dimensional vector space, in order to obtain a regular reconstruction of the missing information. In other words, we want to estimate u : p → u(p) ∈ C(p) while imposing spatial consistency between u(p) and u(q), with q ∈ N(p) where N(p) is the local neighborhood of p.

In the following, we define the variational model and describe how to integrate the spatial constraints within an energy functional. We next present the associated numerical scheme used to minimize the proposed functional.

A. General Variational Model

To address the problem of finding the best candidate for a pixel p ∈ T while adding some spatial consistency, the strategy of [START_REF] Bugeau | Patch-based image colorization[END_REF] consists in first selecting the best candidate (in an arbitrary sense) among the set C(p). Next, during the post-processing step, the solution is regularized to enforce the spatial constraints. Such approach highly depends on the chosen candidate and can induce blurry results. In this work, we propose to reformulate the problem into a single process, i.e., without a first candidate selection step, the selection being induced by the regularization.

1) Data Modeling: In order to introduce the proposed data model, let us now detail what occurs for d = 1dimensional candidates. It is well known that for an ordered set A = {a 1 , . . . , a N } with a i ∈ IR and a i ≤ a i+1 , we have:

median(A) = a [N/2] = min x∈IR N i=1 ||x -a i || 1 = min x∈IR N i=1 |x -a i |,
where [.] computes the integer part and ||.|| 1 stands for the ℓ 1 norm. Therefore selecting, for each point p, the median of the set C(p) of size N is equivalent to minimizing the following energy:

E 1 (u) = p∈T N i=1 |u(p) -c i (p)| (1)
where c i (p) is the i-th candidate of set C(p) and u is the 1dimensional vector to reconstruct. Notice that the number of candidates N could be different for each point p.

A related model has been considered in [START_REF] Graber | Online 3d reconstruction using convex optimization[END_REF] for depth map fusion involving ordered depth candidates of dimension d = 1. Unfortunately, when working with d > 1-dimensional data, the ℓ 1 norm separates the distance between dimensions and the median value does not correspond anymore to one of the candidates. This behavior is a main issue, as we would like to select a single element of the d = 2-dimensional candidates C(p) for each point p ∈ T .

Instead of relying on the median operator through the ℓ 1 norm, we rather consider the ℓ 2 norm to have a metric that compares the d-dimensions conjointly. In order to have a smooth derivable data term, we consider the ℓ 2 2 norm. Our problem of selecting one single candidate per point can be re-casted as follows with i = 1, . . . , N :

E 2 (u) = T min i u -c i 2 2 = T min i d k=1 u k -c k i 2 2 . (2)
In the right side of the equation, u k and c k i stand for the k th dimension of the data. In the image colorization problem,

c 0 i (p) = U (q i ) and c 1 i (p) = V (q i ).
2) Data Term Relaxation: In order to minimize the data term (2), we now consider a relaxation by increasing the state dimension with a new variable w = {w i (p), i = 1, . . . , N } which represents the probability of choosing the i-th candidate for each point p ∈ T , such that

T min i u -c i 2 2 = min w∈W T i w i u -c i 2 2 . (3)
where w has binary values and is defined in the set

W = w(p)={w i (p)} N i=1 , w i : T → {0; 1}, s.t. N i=1 w i (p)=1 .
As W is not a convex set, we consider a relaxation and now define w = (w 1 , . . . , w N ) in the convex set, such as

W = w(p)={w i (p)} N i=1 , w i : T → [0; 1], s.t. N i=1 w i (p)=1 .
Nevertheless, since we can not ensure that the minima of our data model will be reached for binaries values, we add an additional non convex term that enforces the binary nature of w:

E 3 (u, w) = T N i=1 w i u -c i 2 2 + α 2 T N i=1 w i (1 -w i ) (4)
with α > 0. This additional term will enforce the binary conditions as it has two minima in 0 and 1 for w with values in the interval [0; 1].

3) Spatial Constraints and TV Regularization: The energy (4) has trivial minima in its current form. In this paper, we focus on the image colorization problem which involves consistent and spatially coherent image results. To model this behavior, a Total Variation (TV) based regularization of the unknown variable u is introduced in the energy to smooth the solution while keeping sharp discontinuities when required by the data. The final variational model for automatic candidate selection is given by

E(u, w) = T Gu 1 + λ 2 T N i=1 w i u -c i 2 2 + α 2 T N i=1 w i (1 -w i ) ( 5
)
where G is the discrete gradient operator defined w.r.t the chosen neighborhood N (p) with p ∈ T . The parameters λ and α weight the influence of each energy term. 4) Links with the Approach of Arias et al.: Our problem as strong connections with the method of [START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF] dedicated to patch selection for image inpainting. Using our notations, the energy they propose reads:

p∈T N i=1 w i u(p) -c i (p) 2 2,P -α T N i=1 w i log(w i ). ( 6
)
The first term involves a spatial consistency since . 2,P states as the ℓ 2 norm defined on a patch P centered on the pixel p ∈

T . The second term allows to maximize the weight entropy and enforces the weights to be binary. One can note that, in the model ( 6), the candidates c i (p) for a given pixel p ∈ T are chosen within the target image T itself on the contrary of our approach where another image is taken as reference.

The interesting aspect of the model ( 6) is that explicit expressions can be obtained to alternatively update u and w (recovering the NL-means weights for w [START_REF] Buades | A non-local algorithm for image denoising[END_REF]). However, in our current problem involving a TV-based regularization of u, we no longer have an explicit expression for u when computing the associated Euler-Lagrange equations.

We will therefore have to rely on gradient descent based minimization schemes. In the same way, as we would like to recover binary weights, the explicit update of w seems unadapted to our problem. Hence, we will also consider a gradient descent strategy on w.

In this context, the weight entropy model leads to numerical issues since the derivation of this term gives 1 + log(w i ) and tends to infinity when w i tends to 0. An ǫ > 0 parameter should therefore be added to prevent numerical instabilities. Notice that we tested this entropy term. It can give satisfactory results but it involves a difficult tuning of the time-steps parameters in an iterative update process. As a consequence, we rather consider our own non-convex term w i (1w i ).

B. Energy Minimization

In order to minimize (5), we first recall the dual formulation of the TV:

T |Gu| 1 = max z∈B T uG * (y),
where G * is the transpose of the discrete operator G and the dual variable z at point p ∈ T belongs to

IR 2d if u(p) ∈ IR d .
The convex set B is then defined, with k = 1, . . . , d, as:

B = z = (z k,1 , z k,2 ), z k,l : T → IR, s.t. d k=1 2 l=1 z 2 k,l ≤ 1 .
We end up with the following non convex problem:

min u∈A,w∈W max z∈B T uG * z + λ 2 T N i=1 w i u -c i 2 + α 2 T N i=1 w i (1 -w i ) (7)
where the unknown variable u is assumed to take its values in a predefined convex set A characterized with minimal and maximal values for each dimension d:

A = {u = {u k }, k = 1 . . . d and u k : T → [m k ; M k ] ⊂ IR}.
For instance, in the case of color images in Y U V space, m = {0, -0.436, -0.615} and M = {1, 0.436, 0.615}. The problem ( 7) is solved with a primal-dual approach [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] allowing us to compute a local minima of the energy function. The final process is summed up in Algorithm 1. We considered separate implicit gradient descent schemes over the variables u ad w, through the computation of the corresponding proximal operators. This approximation nevertheless leads to convergence in our experiments Algorithm 1 Solving problem [START_REF] Buades | Conditional image diffusion[END_REF] Initialize

w 0 = 1 N , u 0 = N i=1 w 0 i c i , z 0 = 0, ǫ → 0 and choose 0 < τ w < 1 2α and τ u , τ z > 0 such that τ u τ z < 1 8 [24] while ||u k -u k-1 || > ǫ do z k+1 = P B z k + τ z Gu k w k+1 i = P W w k i -τ w (λ u -c i 2 + α) 1 -2ατ w , u k+1 = P A u k + τ u (G * z k+1 + λ M i=1 w k i c i ) 1 + τ u λ end while
The projection operator over the convex set B reads:

P B (z) = z/ max(1, z 2 2 ) with z 2 2 = d k=1 2 l=1 z 2 k,l
. The projection over the set W is obtained via the method of [START_REF] Chen | Projection onto a simplex[END_REF] that allows to project onto a simplex. The projection on the set A is simply:

P A (u k ) = min(M k , max(m k , u k )).
Let us now give the last details necessary to implement the proposed algorithm.

In order to take into account the 2D domain of T of size width × height, the spatial operators G and G * , in Algorithm 1, are defined as follows. The vectors of ma-

trices G = [G x ; G y ] t and G * should satisfy G x u, z 1 + G y u, z 2 = u, (G * x z 1 + G * y z 2 )
. Hence, one can consider a 4-neighborhood N(p) and the finite differences Gu(i, j) = [u +

x (i, j), u + y (i, j)] t , with 1 ≤ i ≤ width and 1 ≤ j ≤ height. With this notation, the first term reads

u + x (i, j) = u(i + 1, j) -u(i, j) if i < width 0 if i = width .
We then obtain G * z = (z 1 ) - x + (z 2 ) - y , where for the first dimension

z - x (i, j) =      z(i, j) if i = 1 z(i, j) -z(i -1, j) if 1 < i < width -z(i -1, j) if i = width .

IV. EXPERIMENTAL RESULTS

Before presenting several colorized images, we provide an experimental analysis of the minimization of energy [START_REF] Luan | Natural image colorization[END_REF]. All the experiments have been run with the same parameters, i.e., τ u = 0.5, τ w = 1, τ z = 0.9/(8τ u ) in Algorithm 1 and λ = 1 and α = 0.5 in [START_REF] Luan | Natural image colorization[END_REF]. For the parameters concerning the colorization method, we have fixed the number of samples to n = 100 for each images, H = 16 for the number of histograms bins and N = 8 for the number of candidates per pixel. For the patches size, we use the following values:

Σ 1 = {3, 5}, Σ 2 = Σ 3 = {7, 9, 11}.

A. Experimental Feature Analysis

In this section, we show experimentally that the same feature is not adapted to all kind of images.

For each result, only one feature with one patch size is used. In order to observe the real influence of each feature, we do not regularize the colorization results. Hence, as there is only one candidate and no regularization, the variational model is not used. Figure 1 presents the results for the following features/distances: variance, discrete Fourier transform, Wasserstein and the SSD distances. The first feature is computed with 3 × 3 patches and the three others with 11 × 11 patches. This experiment shows the interest of combining different features. We can also observe that the candidates provided by SSD are often not well suited for colorization applications. That is the reason why we did not considered this patch metric in the rest of our experiments. By using only 100 sample points in the source image, we are able to recover an image (Figure 2(c)) that is visually close to the original one (Figure 2(a)). The corresponding energy is presented in Figure 2(d). As one can see, the energy decreases as expected but the third term in energy [START_REF] Buades | Conditional image diffusion[END_REF], that enforces the binary nature of the weights, does not reach 0 with the chosen parameters. The final distribution of the weights is shown in Figure 2(e). While at initialization the distribution is uniform (see Algorithm 1), after convergence, we get about 88% of the weights equal to 0. It corresponds to 68% of the pixels having correctly selected one candidate at convergence. For the 32% of pixels left, we can observe that, for a given pixel p, the candidates with w i (p) = 0 may corresponds to the same chrominance candidates due to the restricted initial set of color examples (n = 100). The consequence is that the different metrics may select the same sample. By removing these redundancies, we finally have 90% (Figure 2(f)) of the pixels that have correctly chosen a single color candidate. Moreover, we have observed that by increasing the α and τ w parameters, the third term tends to 0 as expected. The reconstructed image is nevertheless not visually correct with this setting. Indeed, since the initial energy is not convex, the minimization tends to reach a binary local minimum very quickly which does not always correspond to a good candidate satisfying the spatial constraints.

C. Colorization Results

This section presents several colorization results. It is important to notice that all the results presented in this section are obtained with the same set of parameters as described at the beginning of Section IV.

Figure 3 compares our approach with [START_REF] Welsh | Transferring color to greyscale images[END_REF], [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF], [START_REF] Irony | Colorization by example[END_REF] and [START_REF] Bugeau | Patch-based image colorization[END_REF] on two images. Results with [START_REF] Welsh | Transferring color to greyscale images[END_REF] have been obtained by implementing the fully-automatic method (no swatches) while the results of [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF], [START_REF] Irony | Colorization by example[END_REF] and [START_REF] Bugeau | Patch-based image colorization[END_REF] have been directly extracted from the original papers. Our method globally leads to results of the same quality as the ones of [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] or [START_REF] Irony | Colorization by example[END_REF] 1 . The main advantage of our approach is that it requires less pre-processing steps. By zooming on the images, we can observe that our results are smoother. For instance, the house is less green with our method. Furthermore, contrary to [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF], we are able to process the whole image and not only the inner part. Also we want to highlight that the result of [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] presented here has been obtained with a set of multiple reference images. Compared to our previous work in [START_REF] Bugeau | Patch-based image colorization[END_REF], our results are more shiny thanks to the post-processing step and we correctly recover the blue color in the sky of the house image thanks to the variational candidate selection. The second example in Figure 3 presents a comparison with a result extracted from [START_REF] Irony | Colorization by example[END_REF]. In their paper, the authors designed a method for image colorization relying on a predefined image segmentation. Each pixels in the target image only searches for its best match inside one segmented area of the source image. Without needing this first segmentation step, our method is able to recover an image that is visually comparable. For both examples, we clearly outperformed the results of [START_REF] Welsh | Transferring color to greyscale images[END_REF]. Figure 4 illustrates the influence of the source images. Obviously, the choice of input data has an impact on the colorization results. As illustrated in this example, depending on the initial image, the result can be different, e.g. the grass or the leopard body. To overcome this limitation, one could use many source images (as was done in [START_REF] Morimoto | Automatic colorization of grayscale images using multiple images on the web[END_REF] for image colorization and in [START_REF] Torralba | 80 million tiny images: a large database for non-parametric object and scene recognition[END_REF] for image completion). Figure 4(f) shows a colorization results with the two initial examples as input. This result clearly demonstrates the benefit of the multiple image approach. Figure 5 presents other multiple source image examples in order to colorize grayscale data from the three independent images. The first example is inspired from [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] and a comparison can be found in the associated paper. Result with [START_REF] Welsh | Transferring color to greyscale images[END_REF] Result with [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] Result with [START_REF] Bugeau | Patch-based image colorization[END_REF] Our method Source image Target image

Result with [START_REF] Welsh | Transferring color to greyscale images[END_REF] Result with [START_REF] Irony | Colorization by example[END_REF] Result with [START_REF] Bugeau | Patch-based image colorization[END_REF] Our method Fig. 3. Comparisons with other methods [START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] (that uses multiple reference images), [8] [21], and [START_REF] Irony | Colorization by example[END_REF] (see text for more details).

D. Summary of limitations

The proposed method for image colorization lead to very satisfactory results. Nevertheless, it still has several limitations. First, it requires to choose carefully the source images. Secondly, the three chosen features are not well suited to all kind of images. In particular, it seems that some effort should be given to find the good features for old medias. Also, final images are post-processed with luminance remapping which may sometimes create unwanted colors. Finally, one drawback concerns the computational time. Depending on the size of the image, the number of features used and the number of candidates per pixel, the computational time required by all the steps of the proposed colorization pipeline can go from 1 to 10 minutes. For instance, the total processing time for the 500 × 332 image of the house (Figure 3) is about 6 minutes on a standard linux computer. This point can be easily fixed with code optimization and parallel programming.

V. CONCLUSION AND FUTURE WORKS

In this work, a variational approach for image colorization has been proposed. First, several possible candidates are computed using different features and associated metrics. The final color is next obtained by solving a variational model which allows selecting automatically the best candidate while adding some regularization in the solution. Experimental results show that our approach produces visually plausible colorization results and can be competitive with respect to the complex methods proposed in the literature.

As future works, we would like to improve the colorization result to overcome the de-saturation effect of the final result. More concern should indeed be given to the reason why all colorization methods produce drab images. We will also include other complex descriptors or features to improve the characterization of complex image structures. Finally, the extension to video colorization could be considered.

Fig. 1 .

 1 Fig. 1. Influence of the features on the colorization result.

Fig. 2 .

 2 Fig. 2. Analysis of the weights and the energy (see text for details).

Fig. 4 .

 4 Fig. 4. Influence of the source image on the colorization results.

Fig. 5 .

 5 Fig. 5. Colorization with multiple source images.

Figure 6

 6 Figure 6 presents other colorization results. They demonstrate the capacity of our method to colorize different images with exactly the same parameters and without any user intervention (except for the choice of the source image). Finally, the two last examples of Figure 6 are about colorizing old postcard or photos from more recent pictures. Old media have the particularity to be dark and often contain a lot of noise. This makes the colorization process complicated as the features in old and new images do not have the same properties.

The images have been extracted from the article[START_REF] Irony | Colorization by example[END_REF] and the website[START_REF] Charpiat | Automatic image colorization via multimodal predictions[END_REF] and thus may not have the original quality.